PR gradient check

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$$||'`;\
diff --git a/src/cell.c b/src/cell.c
index ff9ff3d9..ee00248d 100644
--- a/src/cell.c
+++ b/src/cell.c
@@ -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;
+}