diff options
author | Thomas White <taw@physics.org> | 2023-06-22 16:16:01 +0200 |
---|---|---|
committer | Thomas White <taw@physics.org> | 2023-07-28 13:22:05 +0200 |
commit | d874916b10e9e5a40b96ca8f878c2a5f7cc60e35 (patch) | |
tree | daf296047798bd94ad47e1b02d7acf4fd2fd3ef2 /tests/gradient_check_utils.c | |
parent | 4977cb7c384f9339037353f31abf0a6a8fd2aea7 (diff) |
Simplified gradient test
prediction_gradient_check was out of control - it had option processing,
as well as bugs of its own! This replaces it with something much
simpler.
Diffstat (limited to 'tests/gradient_check_utils.c')
-rw-r--r-- | tests/gradient_check_utils.c | 185 |
1 files changed, 185 insertions, 0 deletions
diff --git a/tests/gradient_check_utils.c b/tests/gradient_check_utils.c new file mode 100644 index 00000000..f52d43ca --- /dev/null +++ b/tests/gradient_check_utils.c @@ -0,0 +1,185 @@ +/* + * gradient_check_utils.c + * + * Check gradients for prediction refinement (common component) + * + * 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 <gsl/gsl_randist.h> +#include <getopt.h> + +#include <image.h> +#include <cell.h> +#include <cell-utils.h> +#include <geometry.h> +#include <reflist.h> + + +double **make_dev_list(struct reflpeak *rps, int n_refls, struct detgeom *det) +{ + int i; + double **vals; + + vals = malloc(3*sizeof(double *)); + vals[0] = malloc(n_refls*sizeof(double)); + vals[1] = malloc(n_refls*sizeof(double)); + vals[2] = malloc(n_refls*sizeof(double)); + + for ( i=0; i<n_refls; i++ ) { + vals[0][i] = get_exerr(rps[i].refl); + vals[1][i] = fs_dev(&rps[i], det); + vals[2][i] = ss_dev(&rps[i], det); + } + + return vals; +} + + +static UnitCell *random_rotated_cell(gsl_rng *rng) +{ + UnitCell *cell; + struct quaternion orientation; + + cell = cell_new_from_parameters(10.0e-9, 10.0e-9, 10.0e-9, + deg2rad(90.0), + deg2rad(90.0), + deg2rad(90.0)); + orientation = random_quaternion(rng); + return cell_rotate(cell, orientation); +} + + +static void rot(double *x, double *y, double ang) +{ + double nx, ny; + nx = (*x)*cos(ang) - (*y)*sin(ang); + ny = (*x)*sin(ang) + (*y)*cos(ang); + *x = nx; *y = ny; +} + + +struct reflpeak *make_test_image(int *pn_refls, struct image *image) +{ + Crystal *cr; + gsl_rng *rng; + RefList *refls; + Reflection *refl; + RefListIterator *iter; + int n_refls; + int i; + struct reflpeak *rps; + + image->detgeom = malloc(sizeof(struct detgeom)); + image->detgeom->n_panels = 1; + image->detgeom->panels = malloc(sizeof(struct detgeom_panel)); + image->detgeom->panels[0].name = "panel"; + image->detgeom->panels[0].adu_per_photon = 1.0; + image->detgeom->panels[0].max_adu = INFINITY; + image->detgeom->panels[0].fsx = 1.0; + image->detgeom->panels[0].fsy = 0.0; + image->detgeom->panels[0].fsz = 0.0; + image->detgeom->panels[0].ssx = 0.0; + image->detgeom->panels[0].ssy = 1.0; + image->detgeom->panels[0].ssz = 0.0; + rot(&image->detgeom->panels[0].fsx, + &image->detgeom->panels[0].fsy, + deg2rad(30)); + rot(&image->detgeom->panels[0].ssx, + &image->detgeom->panels[0].ssy, + deg2rad(30)); + rot(&image->detgeom->panels[0].fsx, + &image->detgeom->panels[0].fsz, + deg2rad(15)); + rot(&image->detgeom->panels[0].ssx, + &image->detgeom->panels[0].ssz, + deg2rad(15)); + image->detgeom->panels[0].cnx = -500.0; + image->detgeom->panels[0].cny = -500.0; + image->detgeom->panels[0].cnz = 1000.0; /* pixels */ + image->detgeom->panels[0].w = 1000; + image->detgeom->panels[0].h = 1000; + image->detgeom->panels[0].pixel_pitch = 75e-6; + + image->lambda = ph_en_to_lambda(eV_to_J(8000.0)); + image->div = 1e-3; + image->bw = 0.00001; + image->filename = malloc(256); + image->spectrum = spectrum_generate_gaussian(image->lambda, image->bw); + + image->crystals = NULL; + image->n_crystals = 0; + + rng = gsl_rng_alloc(gsl_rng_mt19937); + + cr = crystal_new(); + if ( cr == NULL ) { + ERROR("Failed to allocate crystal.\n"); + return NULL; + } + crystal_set_mosaicity(cr, 0.0); + crystal_set_profile_radius(cr, 0.005e9); + crystal_set_image(cr, image); + crystal_set_cell(cr, random_rotated_cell(rng)); + + refls = predict_to_res(cr, detgeom_max_resolution(image->detgeom, image->lambda)); + crystal_set_reflections(cr, refls); + n_refls = num_reflections(refls); + + /* Associate a peak with every reflection */ + rps = malloc(n_refls*sizeof(struct reflpeak)); + i = 0; + for ( refl = first_refl(refls, &iter); + refl != NULL; + refl = next_refl(refl, iter) ) + { + double fs, ss; + int pn; + struct imagefeature *pk; + + get_detector_pos(refl, &fs, &ss); + pn = get_panel_number(refl); + + pk = malloc(sizeof(struct imagefeature)); + + pk->fs = fs + gsl_ran_gaussian(rng, 5.0); + pk->ss = ss + gsl_ran_gaussian(rng, 5.0); + pk->pn = pn; + pk->intensity = 1.0; + + rps[i].peak = pk; + rps[i].refl = refl; + rps[i].Ih = 1.0; + rps[i].panel = &image->detgeom->panels[pn]; + i++; + } + + image_add_crystal(image, cr); + + gsl_rng_free(rng); + + *pn_refls = n_refls; + return rps; +} |