diff options
| -rw-r--r-- | Makefile.am | 11 | ||||
| -rw-r--r-- | Makefile.in | 35 | ||||
| -rw-r--r-- | src/cell.c | 8 | ||||
| -rw-r--r-- | src/geometry.c | 1 | ||||
| -rw-r--r-- | src/post-refinement.c | 19 | ||||
| -rw-r--r-- | src/post-refinement.h | 20 | ||||
| -rw-r--r-- | src/utils.h | 7 | ||||
| -rw-r--r-- | tests/.gitignore | 1 | ||||
| -rw-r--r-- | tests/pr_gradient_check.c | 270 |
9 files changed, 340 insertions, 32 deletions
diff --git a/Makefile.am b/Makefile.am index 5b71af1e..2c4672ce 100644 --- a/Makefile.am +++ b/Makefile.am @@ -7,11 +7,12 @@ bin_PROGRAMS = src/pattern_sim src/process_hkl src/get_hkl src/indexamajig \ src/calibrate_detector src/partialator \ src/check_hkl src/sum_stack src/partial_sim -noinst_PROGRAMS = tests/list_check tests/integration_check +noinst_PROGRAMS = tests/list_check tests/integration_check \ + tests/pr_gradient_check TESTS = tests/list_check tests/first_merge_check tests/second_merge_check \ tests/third_merge_check tests/fourth_merge_check \ - tests/integration_check + tests/integration_check tests/pr_gradient_check EXTRA_DIST += tests/first_merge_check tests/second_merge_check \ tests/third_merge_check tests/fourth_merge_check \ @@ -134,6 +135,12 @@ tests_integration_check_SOURCES = tests/integration_check.c src/peaks.c \ src/thread-pool.c src/utils.c src/image.c \ src/hdf5-file.c +tests_pr_gradient_check_SOURCES = tests/pr_gradient_check.c src/detector.c \ + src/cell.c src/geometry.c src/reflist.c \ + src/thread-pool.c src/utils.c src/peaks.c \ + src/symmetry.c src/image.c src/hdf5-file.c \ + src/post-refinement.c + INCLUDES = "-I$(top_srcdir)/data" EXTRA_DIST += src/cell.h src/hdf5-file.h src/image.h src/utils.h \ diff --git a/Makefile.in b/Makefile.in index 5a7d3fd4..1413d922 100644 --- a/Makefile.in +++ b/Makefile.in @@ -43,11 +43,12 @@ bin_PROGRAMS = src/pattern_sim$(EXEEXT) src/process_hkl$(EXEEXT) \ src/sum_stack$(EXEEXT) src/partial_sim$(EXEEXT) \ $(am__EXEEXT_1) $(am__EXEEXT_2) noinst_PROGRAMS = tests/list_check$(EXEEXT) \ - tests/integration_check$(EXEEXT) $(am__EXEEXT_3) + tests/integration_check$(EXEEXT) \ + tests/pr_gradient_check$(EXEEXT) $(am__EXEEXT_3) TESTS = tests/list_check$(EXEEXT) tests/first_merge_check \ tests/second_merge_check tests/third_merge_check \ tests/fourth_merge_check tests/integration_check$(EXEEXT) \ - $(am__EXEEXT_3) + tests/pr_gradient_check$(EXEEXT) $(am__EXEEXT_3) @BUILD_HDFSEE_TRUE@am__append_1 = src/hdfsee @BUILD_CUBEIT_TRUE@am__append_2 = src/cubeit @HAVE_OPENCL_TRUE@am__append_3 = src/diffraction-gpu.c src/cl-utils.c @@ -275,6 +276,17 @@ am_tests_list_check_OBJECTS = tests/list_check.$(OBJEXT) \ tests_list_check_OBJECTS = $(am_tests_list_check_OBJECTS) tests_list_check_LDADD = $(LDADD) tests_list_check_DEPENDENCIES = $(top_builddir)/lib/libgnu.a +am_tests_pr_gradient_check_OBJECTS = \ + tests/pr_gradient_check.$(OBJEXT) src/detector.$(OBJEXT) \ + src/cell.$(OBJEXT) src/geometry.$(OBJEXT) \ + src/reflist.$(OBJEXT) src/thread-pool.$(OBJEXT) \ + src/utils.$(OBJEXT) src/peaks.$(OBJEXT) src/symmetry.$(OBJEXT) \ + src/image.$(OBJEXT) src/hdf5-file.$(OBJEXT) \ + src/post-refinement.$(OBJEXT) +tests_pr_gradient_check_OBJECTS = \ + $(am_tests_pr_gradient_check_OBJECTS) +tests_pr_gradient_check_LDADD = $(LDADD) +tests_pr_gradient_check_DEPENDENCIES = $(top_builddir)/lib/libgnu.a DEFAULT_INCLUDES = -I.@am__isrc@ depcomp = $(SHELL) $(top_srcdir)/build-aux/depcomp am__depfiles_maybe = depfiles @@ -303,7 +315,8 @@ SOURCES = $(src_calibrate_detector_SOURCES) $(src_check_hkl_SOURCES) \ $(src_powder_plot_SOURCES) $(src_process_hkl_SOURCES) \ $(src_render_hkl_SOURCES) $(src_sum_stack_SOURCES) \ $(tests_gpu_sim_check_SOURCES) \ - $(tests_integration_check_SOURCES) $(tests_list_check_SOURCES) + $(tests_integration_check_SOURCES) $(tests_list_check_SOURCES) \ + $(tests_pr_gradient_check_SOURCES) DIST_SOURCES = $(src_calibrate_detector_SOURCES) \ $(src_check_hkl_SOURCES) $(src_compare_hkl_SOURCES) \ $(am__src_cubeit_SOURCES_DIST) $(src_get_hkl_SOURCES) \ @@ -313,7 +326,8 @@ DIST_SOURCES = $(src_calibrate_detector_SOURCES) \ $(src_powder_plot_SOURCES) $(src_process_hkl_SOURCES) \ $(src_render_hkl_SOURCES) $(src_sum_stack_SOURCES) \ $(am__tests_gpu_sim_check_SOURCES_DIST) \ - $(tests_integration_check_SOURCES) $(tests_list_check_SOURCES) + $(tests_integration_check_SOURCES) $(tests_list_check_SOURCES) \ + $(tests_pr_gradient_check_SOURCES) RECURSIVE_TARGETS = all-recursive check-recursive dvi-recursive \ html-recursive info-recursive install-data-recursive \ install-dvi-recursive install-exec-recursive \ @@ -711,6 +725,12 @@ tests_integration_check_SOURCES = tests/integration_check.c src/peaks.c \ src/thread-pool.c src/utils.c src/image.c \ src/hdf5-file.c +tests_pr_gradient_check_SOURCES = tests/pr_gradient_check.c src/detector.c \ + src/cell.c src/geometry.c src/reflist.c \ + src/thread-pool.c src/utils.c src/peaks.c \ + src/symmetry.c src/image.c src/hdf5-file.c \ + src/post-refinement.c + INCLUDES = "-I$(top_srcdir)/data" crystfeldir = $(datadir)/crystfel crystfel_DATA = data/diffraction.cl data/defs.h data/hdfsee.ui @@ -976,6 +996,11 @@ tests/list_check.$(OBJEXT): tests/$(am__dirstamp) \ tests/list_check$(EXEEXT): $(tests_list_check_OBJECTS) $(tests_list_check_DEPENDENCIES) tests/$(am__dirstamp) @rm -f tests/list_check$(EXEEXT) $(AM_V_CCLD)$(LINK) $(tests_list_check_OBJECTS) $(tests_list_check_LDADD) $(LIBS) +tests/pr_gradient_check.$(OBJEXT): tests/$(am__dirstamp) \ + tests/$(DEPDIR)/$(am__dirstamp) +tests/pr_gradient_check$(EXEEXT): $(tests_pr_gradient_check_OBJECTS) $(tests_pr_gradient_check_DEPENDENCIES) tests/$(am__dirstamp) + @rm -f tests/pr_gradient_check$(EXEEXT) + $(AM_V_CCLD)$(LINK) $(tests_pr_gradient_check_OBJECTS) $(tests_pr_gradient_check_LDADD) $(LIBS) mostlyclean-compile: -rm -f *.$(OBJEXT) @@ -1022,6 +1047,7 @@ mostlyclean-compile: -rm -f tests/gpu_sim_check.$(OBJEXT) -rm -f tests/integration_check.$(OBJEXT) -rm -f tests/list_check.$(OBJEXT) + -rm -f tests/pr_gradient_check.$(OBJEXT) distclean-compile: -rm -f *.tab.c @@ -1069,6 +1095,7 @@ distclean-compile: @AMDEP_TRUE@@am__include@ @am__quote@tests/$(DEPDIR)/gpu_sim_check.Po@am__quote@ @AMDEP_TRUE@@am__include@ @am__quote@tests/$(DEPDIR)/integration_check.Po@am__quote@ @AMDEP_TRUE@@am__include@ @am__quote@tests/$(DEPDIR)/list_check.Po@am__quote@ +@AMDEP_TRUE@@am__include@ @am__quote@tests/$(DEPDIR)/pr_gradient_check.Po@am__quote@ .c.o: @am__fastdepCC_TRUE@ $(AM_V_CC)depbase=`echo $@ | sed 's|[^/]*$$|$(DEPDIR)/&|;s|\.o$$||'`;\ @@ -614,14 +614,6 @@ void cell_print(UnitCell *cell) #define MAX_CAND (1024) -static int within_tolerance(double a, double b, double percent) -{ - double tol = a * (percent/100.0); - if ( fabs(b-a) < tol ) return 1; - return 0; -} - - static int right_handed(struct rvec a, struct rvec b, struct rvec c) { struct rvec aCb; diff --git a/src/geometry.c b/src/geometry.c index c5e85b12..de2ee1b4 100644 --- a/src/geometry.c +++ b/src/geometry.c @@ -217,6 +217,7 @@ static int check_reflection(struct image *image, double mres, int output, refl = add_refl(reflections, h, k, l); set_detector_pos(refl, 0.0, xda, yda); set_partial(refl, rlow, rhigh, part, clamp_low, clamp_high); + set_symmetric_indices(refl, h, k, l); set_redundancy(refl, 1); if ( output ) { diff --git a/src/post-refinement.c b/src/post-refinement.c index 4ebaffd9..285fe037 100644 --- a/src/post-refinement.c +++ b/src/post-refinement.c @@ -32,23 +32,6 @@ #define MAX_CYCLES (5) -/* Refineable parameters */ -enum { - REF_ASX, - REF_BSX, - REF_CSX, - REF_ASY, - REF_BSY, - REF_CSY, - REF_ASZ, - REF_BSZ, - REF_CSZ, - NUM_PARAMS, - REF_DIV, - REF_R, -}; - - /* Returns dp/dr at "r" */ static double partiality_gradient(double r, double profile_radius) { @@ -86,7 +69,7 @@ static double partiality_rgradient(double r, double profile_radius) /* Return the gradient of parameter 'k' given the current status of 'image'. */ -static double gradient(struct image *image, int k, Reflection *refl, double r) +double gradient(struct image *image, int k, Reflection *refl, double r) { double ds, tt, azix, aziy; double nom, den; diff --git a/src/post-refinement.h b/src/post-refinement.h index 3c8d8587..5b71cbd5 100644 --- a/src/post-refinement.h +++ b/src/post-refinement.h @@ -24,8 +24,28 @@ #include "utils.h" +/* Refineable parameters */ +enum { + REF_ASX, + REF_BSX, + REF_CSX, + REF_ASY, + REF_BSY, + REF_CSY, + REF_ASZ, + REF_BSZ, + REF_CSZ, + NUM_PARAMS, + REF_DIV, + REF_R, +}; + + extern void pr_refine(struct image *image, const RefList *full, const char *sym); +/* Exported so it can be poked by tests/pr_gradient_check */ +extern double gradient(struct image *image, int k, Reflection *refl, double r); + #endif /* POST_REFINEMENT_H */ diff --git a/src/utils.h b/src/utils.h index 4433023f..a11cd59e 100644 --- a/src/utils.h +++ b/src/utils.h @@ -140,6 +140,13 @@ static inline double angle_between(double x1, double y1, double z1, return acos(cosine); } +static inline int within_tolerance(double a, double b, double percent) +{ + double tol = a * (percent/100.0); + if ( fabs(b-a) < fabs(tol) ) return 1; + return 0; +} + /* ----------------------------- Useful macros ------------------------------ */ diff --git a/tests/.gitignore b/tests/.gitignore index 938b5a4a..898fcd89 100644 --- a/tests/.gitignore +++ b/tests/.gitignore @@ -3,4 +3,5 @@ list_check gpu_sim_check integration_check +pr_gradient_check .dirstamp diff --git a/tests/pr_gradient_check.c b/tests/pr_gradient_check.c new file mode 100644 index 00000000..617de4df --- /dev/null +++ b/tests/pr_gradient_check.c @@ -0,0 +1,270 @@ +/* + * pr_gradient_check.c + * + * Check gradients for post refinement + * + * (c) 2011 Thomas White <taw@physics.org> + * + * Part of CrystFEL - crystallography with a FEL + * + */ + + +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif + + +#include <stdlib.h> +#include <stdio.h> + +#include "../src/image.h" +#include "../src/cell.h" +#include "../src/geometry.h" +#include "../src/reflist.h" +#include "../src/post-refinement.h" + + +static void scan_partialities(RefList *reflections, RefList *compare, + int *valid, double *vals[3], int idx) +{ + 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 r1, r2, p; + int clamp_low, clamp_high; + + get_indices(refl, &h, &k, &l); + refl2 = find_refl(compare, h, k, l); + if ( refl2 == NULL ) { + valid[i] = 0; + i++; + continue; + } + + get_partial(refl2, &r1, &r2, &p, &clamp_low, &clamp_high); + if ( clamp_low && clamp_high ) { + if ( !within_tolerance(p, 1.0, 0.001) ) { + + signed int h, k, l; + + get_indices(refl, &h, &k, &l); + + ERROR("%3i %3i %3i - double clamped but" + " partiality not close to 1.0 (%5.2f)\n", + h, k, l, p); + + } + valid[i] = 0; + } + + vals[idx][i] = p; + 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 void calc_either_side(struct image *image, double incr_val, + int *valid, double *vals[3], int refine) +{ + RefList *compare; + UnitCell *cell; + + cell = new_shifted_cell(image->indexed_cell, refine, -incr_val); + compare = find_intersections(image, cell, 0); + scan_partialities(image->reflections, compare, valid, vals, 0); + cell_free(cell); + + cell = new_shifted_cell(image->indexed_cell, refine, +incr_val); + compare = find_intersections(image, cell, 0); + scan_partialities(image->reflections, compare, valid, vals, 2); + cell_free(cell); +} + + +static int test_gradients(struct image *image, double incr_val, int refine, + const char *str) +{ + Reflection *refl; + RefListIterator *iter; + double *vals[3]; + int i; + int *valid; + int nref; + int n_acc, n_valid; + + image->reflections = find_intersections(image, image->indexed_cell, 0); + + nref = num_reflections(image->reflections); + if ( nref < 10 ) { + ERROR("Too few reflections found. Failing test by default.\n"); + return -1; + } + + vals[0] = malloc(nref*sizeof(double)); + vals[1] = malloc(nref*sizeof(double)); + vals[2] = malloc(nref*sizeof(double)); + if ( (vals[0] == NULL) || (vals[1] == NULL) || (vals[2] == NULL) ) { + ERROR("Couldn't allocate memory.\n"); + return -1; + } + + valid = malloc(nref*sizeof(int)); + if ( valid == NULL ) { + ERROR("Couldn't allocate memory.\n"); + return -1; + } + for ( i=0; i<nref; i++ ) valid[i] = 1; + + scan_partialities(image->reflections, image->reflections, + valid, vals, 1); + + calc_either_side(image, incr_val, valid, vals, refine); + + n_valid = nref; n_acc = 0; + i = 0; + for ( refl = first_refl(image->reflections, &iter); + refl != NULL; + refl = next_refl(refl, iter) ) + { + + double grad1, grad2, grad; + double cgrad; + signed int h, k, l; + + get_indices(refl, &h, &k, &l); + + if ( !valid[i] ) { + n_valid--; + } else { + + grad1 = (vals[1][i] - vals[0][i]) / incr_val; + grad2 = (vals[2][i] - vals[1][i]) / incr_val; + grad = (grad1 + grad2) / 2.0; + + cgrad = gradient(image, refine, refl, + image->profile_radius); + + if ( !within_tolerance(grad, cgrad, 10.0) ) { + STATUS("!- %s %5.2f %5.2f %5.2f" + " -> %10.2e %10.2e -> %10.2e %10.2e\n", + str, vals[0][i], vals[1][i], vals[2][i], + grad1, grad2, grad, cgrad); + } else { + //STATUS("OK %s %5.2f %5.2f %5.2f" + // " -> %10.2e %10.2e" + // " -> %10.2e %10.2e\n", + // str, vals[0][i], vals[1][i], vals[2][i], + // grad1, grad2, grad, cgrad); + n_acc++; + } + + } + + i++; + + } + + STATUS("%s: %i out of %i valid gradients were accurate.\n", + str, n_acc, n_valid); + + if ( n_acc != n_valid ) return 1; + + return 0; +} + +int main(int argc, char *argv[]) +{ + struct image image; + const double incr_frac = 1.0/1000.0; + double incr_val; + double ax, ay, az; + double bx, by, bz; + double cx, cy, cz; + UnitCell *cell; + struct quaternion orientation; + + image.width = 1024; + image.height = 1024; + image.det = simple_geometry(&image); + + image.lambda = ph_en_to_lambda(eV_to_J(2000.0)); + image.div = 0.009; + image.bw = 0.01; + image.m = 0.0; + image.profile_radius = 0.005e9; + image.i0_available = 0; + image.filename = malloc(256); + + 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(); + image.indexed_cell = cell_rotate(cell, orientation); + + cell_get_reciprocal(image.indexed_cell, + &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz); + + incr_val = incr_frac * ax; + test_gradients(&image, incr_val, REF_ASX, "ax*"); + incr_val = incr_frac * ay; + test_gradients(&image, incr_val, REF_ASY, "ay*"); + incr_val = incr_frac * az; + test_gradients(&image, incr_val, REF_ASZ, "az*"); + + incr_val = incr_frac * bx; + test_gradients(&image, incr_val, REF_BSX, "bx*"); + incr_val = incr_frac * by; + test_gradients(&image, incr_val, REF_BSY, "by*"); + incr_val = incr_frac * bz; + test_gradients(&image, incr_val, REF_BSZ, "bz*"); + + incr_val = incr_frac * cx; + test_gradients(&image, incr_val, REF_CSX, "cx*"); + incr_val = incr_frac * cy; + test_gradients(&image, incr_val, REF_CSY, "cy*"); + incr_val = incr_frac * cz; + test_gradients(&image, incr_val, REF_CSZ, "cz*"); + + return 0; +} |