diff options
Diffstat (limited to 'tests/gradient_check.c')
-rw-r--r-- | tests/gradient_check.c | 187 |
1 files changed, 187 insertions, 0 deletions
diff --git a/tests/gradient_check.c b/tests/gradient_check.c new file mode 100644 index 00000000..001bff7a --- /dev/null +++ b/tests/gradient_check.c @@ -0,0 +1,187 @@ +/* + * gradient_panel_x.c + * + * Check gradients for prediction refinement + * + * Copyright © 2012-2023 Deutsches Elektronen-Synchrotron DESY, + * a research centre of the Helmholtz Association. + * + * Authors: + * 2012-2023 Thomas White <taw@physics.org> + * + * This file is part of CrystFEL. + * + * CrystFEL is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * CrystFEL is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with CrystFEL. If not, see <http://www.gnu.org/licenses/>. + * + */ + +#include <stdlib.h> +#include <stdio.h> +#include <getopt.h> + +#include <image.h> +#include <geometry.h> +#include <predict-refine.h> + +#include "gradient_check_utils.h" + + +int main(int argc, char *argv[]) +{ + struct image image; + struct reflpeak *rps; + int n_refls; + double **before; + double **after; + int i; + int n_wrong_r = 0; + int n_wrong_fs = 0; + int n_wrong_ss = 0; + int n_wrong_obsr = 0; + int fail = 0; + double step; + gsl_matrix **panel_matrices; + int didsomething = 0; + const double cx = 0.03; /* Detector is a 7.5 cm side length square */ + const double cy = 0.02; + const double cz = 0.01; + + rps = make_test_image(&n_refls, &image); + panel_matrices = make_panel_minvs(image.detgeom); + + before = make_dev_list(rps, n_refls, image.detgeom); + + #ifdef TRANSLATE_PANEL + struct detgeom_panel *p = &image.detgeom->panels[0]; + step = 0.01e-3; /* metres */ + image.detgeom->panels[0].THING_TO_MOVE += step/p->pixel_pitch; + didsomething = 1; + #endif + + #ifdef ROTATE_PANEL_X + struct detgeom_panel *p = &image.detgeom->panels[0]; + step = deg2rad(0.01); + rotate2d(&p->cny, &p->cnz, cy/p->pixel_pitch, cz/p->pixel_pitch, step); + rotate2d(&p->fsy, &p->fsz, 0, 0, step); + rotate2d(&p->ssy, &p->ssz, 0, 0, step); + didsomething = 1; + #endif + + #ifdef ROTATE_PANEL_Y + struct detgeom_panel *p = &image.detgeom->panels[0]; + step = deg2rad(0.01); + rotate2d(&p->cnz, &p->cnx, cz/p->pixel_pitch, cx/p->pixel_pitch, step); + rotate2d(&p->fsz, &p->fsx, 0, 0, step); + rotate2d(&p->ssz, &p->ssx, 0, 0, step); + didsomething = 1; + #endif + + #ifdef ROTATE_PANEL_Z + struct detgeom_panel *p = &image.detgeom->panels[0]; + step = deg2rad(0.01); + rotate2d(&p->cnx, &p->cny, cx/p->pixel_pitch, cy/p->pixel_pitch, step); + rotate2d(&p->fsx, &p->fsy, 0, 0, step); + rotate2d(&p->ssx, &p->ssy, 0, 0, step); + didsomething = 1; + #endif + + #ifdef CHANGE_CELL + double asx, asy, asz, bsx, bsy, bsz, csx, csy, csz; + UnitCell *cell = crystal_get_cell(image.crystals[0]); + step = 0.5e5; + cell_get_reciprocal(cell, &asx, &asy, &asz, + &bsx, &bsy, &bsz, + &csx, &csy, &csz); + THING_TO_MOVE += step; + cell_set_reciprocal(cell, asx, asy, asz, + bsx, bsy, bsz, + csx, csy, csz); + didsomething = 1; + #endif + + if ( !didsomething ) { + fprintf(stderr, "Nothing changed. Check the build system.\n"); + return 1; + } + + update_predictions(image.crystals[0]); + after = make_dev_list(rps, n_refls, image.detgeom); + + for ( i=0; i<n_refls; i++ ) { + + float calc[3]; + double obs[3]; + + calc[0] = r_gradient(TEST_GPARAM, rps[i].refl, + crystal_get_cell(image.crystals[0]), + image.lambda); + + fs_ss_gradient(TEST_GPARAM, rps[i].refl, + crystal_get_cell(image.crystals[0]), + &image.detgeom->panels[rps[i].peak->pn], + panel_matrices[rps[i].peak->pn], cx, cy, cz, + &calc[1], &calc[2]); + + obs[0] = (after[0][i] - before[0][i]) / step; + obs[1] = (after[1][i] - before[1][i]) / step; + obs[2] = (after[2][i] - before[2][i]) / step; + + #ifdef TRANSLATE_PANEL + if ( fabs(calc[0]) > 1e-12 ) n_wrong_r++; /* Should be zero */ + if ( fabs(obs[0]) > 1e-12 ) n_wrong_obsr++; /* Should also be zero */ + if ( fabs(obs[1] - calc[1]) > 10.0 ) n_wrong_fs++; + if ( fabs(obs[2] - calc[2]) > 10.0 ) n_wrong_ss++; + #endif + + #if defined(ROTATE_PANEL_X) || defined(ROTATE_PANEL_Y) || defined(ROTATE_PANEL_Z) + if ( fabs(calc[0]) > 1e-12 ) n_wrong_r++; /* Should be zero */ + if ( fabs(obs[0]) > 1e-12 ) n_wrong_obsr++; /* Should also be zero */ + if ( fabs(obs[1] - calc[1]) > 1.0 ) n_wrong_fs++; /* Units are pixels/rad */ + if ( fabs(obs[2] - calc[2]) > 1.0 ) n_wrong_ss++; /* (numbers are big) */ + #endif + + #ifdef CHANGE_CELL + if ( fabs(obs[0] - calc[0]) > 1e-2 ) n_wrong_r++; + if ( fabs(obs[1] - calc[1]) > 1e-8 ) n_wrong_fs++; + if ( fabs(obs[2] - calc[2]) > 1e-8 ) n_wrong_ss++; + #endif + + } + + if ( n_wrong_r > 0 ) { + fprintf(stderr, "%i out of %i R gradients were wrong.\n", + n_wrong_r, n_refls); + fail = 1; + } + + if ( n_wrong_fs > 0 ) { + fprintf(stderr, "%i out of %i fs gradients were wrong.\n", + n_wrong_fs, n_refls); + fail = 1; + } + + if ( n_wrong_ss > 0 ) { + fprintf(stderr, "%i out of %i ss gradients were wrong.\n", + n_wrong_ss, n_refls); + fail = 1; + } + + if ( n_wrong_obsr > 0 ) { + fprintf(stderr, "%i out of %i observed R gradients were not zero as expected\n", + n_wrong_obsr, n_refls); + fail = 1; + } + + return fail; +} |