From 8a8e940575562086ab166e88a1f1927f420c564d Mon Sep 17 00:00:00 2001 From: Thomas White Date: Mon, 16 Mar 2015 17:16:10 +0100 Subject: Add gradient test program and fix gradients --- Makefile.am | 8 +- src/predict-refine.c | 36 +-- tests/prediction_gradient_check.c | 498 ++++++++++++++++++++++++++++++++++++++ 3 files changed, 524 insertions(+), 18 deletions(-) create mode 100644 tests/prediction_gradient_check.c diff --git a/Makefile.am b/Makefile.am index 73d9a89e..b340dc86 100644 --- a/Makefile.am +++ b/Makefile.am @@ -11,7 +11,8 @@ noinst_PROGRAMS = tests/list_check tests/integration_check \ tests/pr_p_gradient_check tests/symmetry_check \ tests/centering_check tests/transformation_check \ tests/cell_check tests/ring_check \ - tests/prof2d_check tests/ambi_check + tests/prof2d_check tests/ambi_check \ + tests/prediction_gradient_check MERGE_CHECKS = tests/first_merge_check tests/second_merge_check \ tests/third_merge_check tests/fourth_merge_check @@ -24,7 +25,8 @@ PARTIAL_CHECKS = tests/partialator_merge_check_1 \ TESTS = tests/list_check $(MERGE_CHECKS) $(PARTIAL_CHECKS) \ tests/integration_check \ tests/symmetry_check tests/centering_check tests/transformation_check \ - tests/cell_check tests/ring_check tests/prof2d_check tests/ambi_check + tests/cell_check tests/ring_check tests/prof2d_check tests/ambi_check \ + tests/prediction_gradient_check EXTRA_DIST += $(MERGE_CHECKS) $(PARTIAL_CHECKS) EXTRA_DIST += relnotes-0.6.0 @@ -122,6 +124,8 @@ tests_ambi_check_SOURCES = tests/ambi_check.c tests_pr_p_gradient_check_SOURCES = tests/pr_p_gradient_check.c \ src/post-refinement.c +tests_prediction_gradient_check_SOURCES = tests/prediction_gradient_check.c + tests_centering_check_SOURCES = tests/centering_check.c tests_transformation_check_SOURCES = tests/transformation_check.c diff --git a/src/predict-refine.c b/src/predict-refine.c index 5d0c37db..dee1ad55 100644 --- a/src/predict-refine.c +++ b/src/predict-refine.c @@ -120,8 +120,8 @@ static void twod_mapping(double fs, double ss, double *px, double *py, xs = fs*p->fsx + ss*p->ssx; ys = fs*p->fsy + ss*p->ssy; - *px = xs + p->cnx; - *py = ys + p->cny; + *px = (xs + p->cnx) / p->res; + *py = (ys + p->cny) / p->res; } @@ -171,11 +171,8 @@ static double r_gradient(UnitCell *cell, int k, Reflection *refl, azi = atan2(yl, xl); - /* For many gradients, just multiply the above number by the gradient - * of excitation error wrt whatever. */ switch ( k ) { - /* Cell parameters and orientation */ case REF_ASX : return - hs * sin(phi) * cos(azi); @@ -210,46 +207,53 @@ static double r_gradient(UnitCell *cell, int k, Reflection *refl, } -static double pos_gradient(int param, struct reflpeak *rp, struct detector *det) +/* Returns d(xh-xpk)/dP + d(yh-ypk)/dP, where P = any parameter */ +static double pos_gradient(int param, struct reflpeak *rp, struct detector *det, + double lambda, UnitCell *cell) { signed int h, k, l; double xpk, ypk, xh, yh; double fsh, ssh; + double tt, clen, azi, azf; twod_mapping(rp->peak->fs, rp->peak->ss, &xpk, &ypk, det); get_detector_pos(rp->refl, &fsh, &ssh); twod_mapping(fsh, ssh, &xh, &yh, det); get_indices(rp->refl, &h, &k, &l); + tt = asin(lambda * resolution(cell, h, k, l)); + clen = find_panel(det, fsh, ssh)->clen; + azi = atan2(yh, xh); + azf = 2.0*(cos(azi) + sin(azi)); /* FIXME: Why factor of 2? */ + switch ( param ) { - /* Cell parameters and orientation */ case REF_ASX : - return 2.0 * h * (xh-xpk); + return h * lambda * clen / cos(tt); case REF_BSX : - return 2.0 * k * (xh-xpk); + return k * lambda * clen / cos(tt); case REF_CSX : - return 2.0 * l * (xh-xpk); + return l * lambda * clen / cos(tt); case REF_ASY : - return 2.0 * h * (yh-ypk); + return h * lambda * clen / cos(tt); case REF_BSY : - return 2.0 * k * (yh-ypk); + return k * lambda * clen / cos(tt); case REF_CSY : - return 2.0 * l * (yh-ypk); + return l * lambda * clen / cos(tt); case REF_ASZ : - return 0.0; + return -h * lambda * clen * azf * sin(tt) / (cos(tt)*cos(tt)); case REF_BSZ : - return 0.0; + return -k * lambda * clen * azf * sin(tt) / (cos(tt)*cos(tt)); case REF_CSZ : - return 0.0; + return -l * lambda * clen * azf * sin(tt) / (cos(tt)*cos(tt)); } diff --git a/tests/prediction_gradient_check.c b/tests/prediction_gradient_check.c new file mode 100644 index 00000000..09b4a868 --- /dev/null +++ b/tests/prediction_gradient_check.c @@ -0,0 +1,498 @@ +/* + * prediction_gradient_check.c + * + * Check partiality gradients for prediction refinement + * + * Copyright © 2012-2015 Deutsches Elektronen-Synchrotron DESY, + * a research centre of the Helmholtz Association. + * + * Authors: + * 2012-2015 Thomas White + * + * 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 . + * + */ + +#ifdef HAVE_CONFIG_H +#include +#endif + + +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include "../src/predict-refine.c" + + +int checkrxy; + + +static void scan(RefList *reflections, RefList *compare, + int *valid, long double *vals[3], int idx, + struct detector *det) +{ + int i; + Reflection *refl; + RefListIterator *iter; + + i = 0; + for ( refl = first_refl(reflections, &iter); + refl != NULL; + refl = next_refl(refl, iter) ) + { + signed int h, k, l; + Reflection *refl2; + double rlow, rhigh, p; + double fs, ss, xh, yh; + + get_indices(refl, &h, &k, &l); + refl2 = find_refl(compare, h, k, l); + if ( refl2 == NULL ) { + valid[i] = 0; + i++; + continue; + } + + get_partial(refl2, &rlow, &rhigh, &p); + get_detector_pos(refl2, &fs, &ss); + twod_mapping(fs, ss, &xh, &yh, det); + + switch ( checkrxy ) { + + case 0 : + vals[idx][i] = (rlow + rhigh)/2.0; + break; + + case 1 : + vals[idx][i] = xh + yh; + break; + } + + i++; + } +} + + +static UnitCell *new_shifted_cell(UnitCell *input, int k, double shift) +{ + UnitCell *cell; + double asx, asy, asz; + double bsx, bsy, bsz; + double csx, csy, csz; + + cell = cell_new(); + cell_get_reciprocal(input, &asx, &asy, &asz, &bsx, &bsy, &bsz, + &csx, &csy, &csz); + switch ( k ) + { + case REF_ASX : asx += shift; break; + case REF_ASY : asy += shift; break; + case REF_ASZ : asz += shift; break; + case REF_BSX : bsx += shift; break; + case REF_BSY : bsy += shift; break; + case REF_BSZ : bsz += shift; break; + case REF_CSX : csx += shift; break; + case REF_CSY : csy += shift; break; + case REF_CSZ : csz += shift; break; + } + cell_set_reciprocal(cell, asx, asy, asz, bsx, bsy, bsz, csx, csy, csz); + + return cell; +} + + +static Crystal *new_shifted_crystal(Crystal *cr, int refine, double incr_val) +{ + Crystal *cr_new; + UnitCell *cell; + + cr_new = crystal_copy(cr); + if ( cr_new == NULL ) { + ERROR("Failed to allocate crystal.\n"); + return NULL; + } + + crystal_set_image(cr_new, crystal_get_image(cr)); + + switch ( refine ) { + + case REF_ASX : + case REF_ASY : + case REF_ASZ : + case REF_BSX : + case REF_BSY : + case REF_BSZ : + case REF_CSX : + case REF_CSY : + case REF_CSZ : + cell = new_shifted_cell(crystal_get_cell(cr), refine, + incr_val); + crystal_set_cell(cr_new, cell); + break; + + default: + ERROR("Can't shift %i\n", refine); + break; + + } + + return cr_new; +} + + +static void calc_either_side(Crystal *cr, double incr_val, + int *valid, long double *vals[3], int refine) +{ + RefList *compare; + struct image *image = crystal_get_image(cr); + Crystal *cr_new; + + cr_new = new_shifted_crystal(cr, refine, -incr_val); + compare = find_intersections(image, cr_new, PMODEL_SCSPHERE); + scan(crystal_get_reflections(cr), compare, valid, vals, 0, + crystal_get_image(cr)->det); + cell_free(crystal_get_cell(cr_new)); + crystal_free(cr_new); + reflist_free(compare); + + cr_new = new_shifted_crystal(cr, refine, +incr_val); + compare = find_intersections(image, cr_new, PMODEL_SCSPHERE); + scan(crystal_get_reflections(cr), compare, valid, vals, 2, + crystal_get_image(cr)->det); + cell_free(crystal_get_cell(cr_new)); + crystal_free(cr_new); + reflist_free(compare); +} + + +static double test_gradients(Crystal *cr, double incr_val, int refine, + const char *str, const char *file, + int quiet, int plot) +{ + Reflection *refl; + RefListIterator *iter; + long double *vals[3]; + int i; + int *valid; + int nref; + int n_good, n_invalid, n_small, n_nan, n_bad; + RefList *reflections; + FILE *fh = NULL; + int ntot = 0; + double total = 0.0; + char tmp[32]; + double *vec1; + double *vec2; + int n_line; + double cc; + + reflections = find_intersections(crystal_get_image(cr), cr, + PMODEL_SCSPHERE); + crystal_set_reflections(cr, reflections); + + nref = num_reflections(reflections); + if ( nref < 10 ) { + ERROR("Too few reflections found. Failing test by default.\n"); + return 0.0; + } + + vals[0] = malloc(nref*sizeof(long double)); + vals[1] = malloc(nref*sizeof(long double)); + vals[2] = malloc(nref*sizeof(long double)); + if ( (vals[0] == NULL) || (vals[1] == NULL) || (vals[2] == NULL) ) { + ERROR("Couldn't allocate memory.\n"); + return 0.0; + } + + valid = malloc(nref*sizeof(int)); + if ( valid == NULL ) { + ERROR("Couldn't allocate memory.\n"); + return 0.0; + } + for ( i=0; idet); + + calc_either_side(cr, incr_val, valid, vals, refine); + + if ( plot ) { + snprintf(tmp, 32, "gradient-test-%s.dat", file); + fh = fopen(tmp, "w"); + } + + vec1 = malloc(nref*sizeof(double)); + vec2 = malloc(nref*sizeof(double)); + if ( (vec1 == NULL) || (vec2 == NULL) ) { + ERROR("Couldn't allocate memory.\n"); + return 0.0; + } + + n_invalid = 0; n_good = 0; + n_nan = 0; n_small = 0; n_bad = 0; n_line = 0; + i = 0; + for ( refl = first_refl(reflections, &iter); + refl != NULL; + refl = next_refl(refl, iter) ) + { + long double grad1, grad2, grad; + double cgrad; + signed int h, k, l; + struct reflpeak rp; + + get_indices(refl, &h, &k, &l); + + if ( !valid[i] ) { + n_invalid++; + i++; + } else { + + double r1, r2, p; + + grad1 = (vals[1][i] - vals[0][i]) / incr_val; + grad2 = (vals[2][i] - vals[1][i]) / incr_val; + grad = (grad1 + grad2) / 2.0; + i++; + + if ( checkrxy == 0 ) { + + cgrad = r_gradient(crystal_get_cell(cr), refine, + refl, crystal_get_image(cr)); + + } else { + + struct imagefeature pk; + + pk.fs = 0.0; + pk.ss = 0.0; + rp.refl = refl; + rp.peak = &pk; + + cgrad = pos_gradient(refine, &rp, + crystal_get_image(cr)->det, + crystal_get_image(cr)->lambda, + crystal_get_cell(cr)); + } + + get_partial(refl, &r1, &r2, &p); + + if ( isnan(cgrad) ) { + n_nan++; + continue; + } + + if ( plot ) { + fprintf(fh, "%e %Le\n", cgrad, grad); + } + + vec1[n_line] = cgrad; + vec2[n_line] = grad; + n_line++; + + if ( (fabs(cgrad) < 5e-12) && (fabs(grad) < 5e-12) ) { + n_small++; + continue; + } + + total += fabs(cgrad - grad); + ntot++; + + if ( !within_tolerance(grad, cgrad, 5.0) + || !within_tolerance(cgrad, grad, 5.0) ) + { + + if ( !quiet ) { + STATUS("!- %s %3i %3i %3i" + " %10.2Le %10.2e ratio = %5.2Lf" + " %10.2e %10.2e\n", + str, h, k, l, grad, cgrad, + cgrad/grad, r1, r2); + } + n_bad++; + + } else { + + STATUS("OK %s %3i %3i %3i" + " %10.2Le %10.2e ratio = %5.2Lf" + " %10.2e %10.2e\n", + str, h, k, l, grad, cgrad, cgrad/grad, + r1, r2); + + n_good++; + + } + + } + + } + + STATUS("%3s: %3i within 5%%, %3i outside, %3i nan, %3i invalid, " + "%3i small. ", str, n_good, n_bad, n_nan, n_invalid, n_small); + + if ( plot ) { + fclose(fh); + } + + cc = gsl_stats_correlation(vec1, 1, vec2, 1, n_line); + STATUS("CC = %+f\n", cc); + return cc; +} + + +int main(int argc, char *argv[]) +{ + struct image image; + const double incr_frac = 1.0/100000.0; + double incr_val; + double ax, ay, az; + double bx, by, bz; + double cx, cy, cz; + UnitCell *cell; + Crystal *cr; + struct quaternion orientation; + int fail = 0; + int quiet = 0; + int plot = 0; + int c; + gsl_rng *rng; + UnitCell *rot; + double val; + + const struct option longopts[] = { + {"quiet", 0, &quiet, 1}, + {"plot", 0, &plot, 1}, + {0, 0, NULL, 0} + }; + + while ((c = getopt_long(argc, argv, "", longopts, NULL)) != -1) { + switch (c) { + + case 0 : + break; + + case '?' : + break; + + default : + ERROR("Unhandled option '%c'\n", c); + break; + + } + + } + + image.width = 1024; + image.height = 1024; + image.det = simple_geometry(&image); + image.det->panels[0].res = 13333.3; + image.det->panels[0].clen = 80e-3; + image.det->panels[0].coffset = 0.0; + + image.lambda = ph_en_to_lambda(eV_to_J(8000.0)); + image.div = 1e-3; + image.bw = 0.01; + image.filename = malloc(256); + + cr = crystal_new(); + if ( cr == NULL ) { + ERROR("Failed to allocate crystal.\n"); + return 1; + } + crystal_set_mosaicity(cr, 0.0); + crystal_set_profile_radius(cr, 0.005e9); + crystal_set_image(cr, &image); + + cell = cell_new_from_parameters(10.0e-9, 10.0e-9, 10.0e-9, + deg2rad(90.0), + deg2rad(90.0), + deg2rad(90.0)); + + rng = gsl_rng_alloc(gsl_rng_mt19937); + + for ( checkrxy=1; checkrxy<2; checkrxy++ ) { + + + switch ( checkrxy ) { + case 0 : + STATUS("Excitation error:\n"); + break; + case 1: + STATUS("x+y coordinate:\n"); + break; + default: + STATUS("WTF??\n"); + break; + } + + orientation = random_quaternion(rng); + rot = cell_rotate(cell, orientation); + crystal_set_cell(cr, rot); + + cell_get_reciprocal(rot, &ax, &ay, &az, + &bx, &by, &bz, &cx, &cy, &cz); + + incr_val = incr_frac * ax; + val = test_gradients(cr, incr_val, REF_ASX, "ax*", "ax", + quiet, plot); + if ( val < 0.99 ) fail = 1; + incr_val = incr_frac * bx; + val = test_gradients(cr, incr_val, REF_BSX, "bx*", "bx", + quiet, plot); + if ( val < 0.99 ) fail = 1; + incr_val = incr_frac * cx; + val = test_gradients(cr, incr_val, REF_CSX, "cx*", "cx", + quiet, plot); + if ( val < 0.99 ) fail = 1; + + incr_val = incr_frac * ay; + val = test_gradients(cr, incr_val, REF_ASY, "ay*", "ay", + quiet, plot); + if ( val < 0.99 ) fail = 1; + incr_val = incr_frac * by; + val = test_gradients(cr, incr_val, REF_BSY, "by*", "by", + quiet, plot); + if ( val < 0.99 ) fail = 1; + incr_val = incr_frac * cy; + val = test_gradients(cr, incr_val, REF_CSY, "cy*", "cy", + quiet, plot); + if ( val < 0.99 ) fail = 1; + + incr_val = incr_frac * az; + val = test_gradients(cr, incr_val, REF_ASZ, "az*", "az", + quiet, plot); + if ( val < 0.99 ) fail = 1; + incr_val = incr_frac * bz; + val = test_gradients(cr, incr_val, REF_BSZ, "bz*", "bz", + quiet, plot); + if ( val < 0.99 ) fail = 1; + incr_val = incr_frac * cz; + val = test_gradients(cr, incr_val, REF_CSZ, "cz*", "cz", + quiet, plot); + if ( val < 0.99 ) fail = 1; + + } + + gsl_rng_free(rng); + + return fail; +} -- cgit v1.2.3