aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--Makefile.am8
-rw-r--r--src/predict-refine.c36
-rw-r--r--tests/prediction_gradient_check.c498
3 files changed, 524 insertions, 18 deletions
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 <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/>.
+ *
+ */
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+
+#include <stdlib.h>
+#include <stdio.h>
+#include <gsl/gsl_statistics.h>
+#include <getopt.h>
+
+#include <image.h>
+#include <cell.h>
+#include <cell-utils.h>
+#include <geometry.h>
+#include <reflist.h>
+#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; i<nref; i++ ) valid[i] = 1;
+
+ scan(reflections, reflections, valid, vals, 1,
+ crystal_get_image(cr)->det);
+
+ 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;
+}
generated by cgit v1.2.3 (git 2.25.1) at 2024-06-30 16:58:04 +0000