Refinement test framework

git-svn-id: svn://cook.msm.cam.ac.uk:745/diff-tomo/dtr@271 bf6ca9ba-c028-0410-8290-897cf20841d1

6 files changed, 117 insertions, 20 deletions

diff --git a/Makefile.am b/Makefile.am
index 9ac3ebb..3f6d2bf 100644
--- a/Makefile.am
+++ b/Makefile.am
@@ -3,3 +3,5 @@ EXTRA_DIST = configure src/displaywindow.h src/trackball.h src/reflections.h src
 		src/itrans-zaefferer.h src/itrans-stat.h src/mapping.h src/reproject.h src/prealign.h \
 		src/dirax.h src/image.h src/refine.h src/gtk-valuegraph.h src/intensities.h src/glbits.h src/gtk-valuegraph.h
 SUBDIRS = src data
+TESTS = src/refinetest
+
diff --git a/src/Makefile.am b/src/Makefile.am
index e9aae1c..ea94888 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -1,10 +1,13 @@
-bin_PROGRAMS = dtr
+bin_PROGRAMS = dtr refinetest
 
 dtr_SOURCES = main.c displaywindow.c trackball.c reflections.c readpng.c mrc.c imagedisplay.c utils.c itrans.c qdrp.c cache.c \
 		itrans-threshold.c itrans-zaefferer.c itrans-stat.c control.c mapping.c reproject.c prealign.c basis.c \
 		dirax.c image.c refine.c gtk-valuegraph.c intensities.c glbits.c
 dtr_LDADD = @LIBS@ @GTK_LIBS@ -lm @GTKGLEXT_LIBS@ -lgsl -lgslcblas -lutil
 
+refinetest_SOURCES = refinetest.c reflections.c basis.c utils.c reproject.c image.c refine.c mapping.c control.c
+refinetest_LDADD = @LIBS@ @GTK_LIBS@ -lm -lgsl -lgslcblas -lutil
+
 AM_CFLAGS = -Wall -g @CFLAGS@ @GTK_CFLAGS@ @GTKGLEXT_CFLAGS@
 AM_CPPFLAGS = -DDATADIR=\""$(datadir)"\"
 
diff --git a/src/refine.c b/src/refine.c
index b56cc05..7bb2585 100644
--- a/src/refine.c
+++ b/src/refine.c
@@ -93,11 +93,11 @@ void refine_do_sequence(ControlContext *ctx) {
 	gtk_widget_show_all(window_fit);
 	
 	/* Perform final refinement */
-	refine_do_cell(ctx);
 	printf("Best omega offset = %f deg (%f nm^-1)\n", rad2deg(omega_offs_best), fit_best/1e9);
 	for ( j=0; j<ctx->images->n_images; j++ ) {
 		ctx->images->images[j].omega += omega_offs_best;
 	}
+	refine_do_cell(ctx);
 	reproject_lattice_changed(ctx);
 	mapping_adjust_axis(ctx, omega_offs_best);
 
@@ -165,13 +165,11 @@ static void refine_simplex_transform(SimplexVertex *s, int v_worst, double fac)
 	
 }
 
-static double refine_iteration(SimplexVertex *s, Deviation *d, int nf) {
+static double refine_iteration(SimplexVertex *s, Deviation *d, int nf, int debug) {
 
 	int v_worst, v_best, v_second_worst, i;
 	double fom_worst, fom_new, fom_best, fom_second_worst;
-
-//	printf("---\n");
-
+	
 	/* Find the least favourable vertex of the simplex */
 	v_worst = 0;
 	fom_worst = 0.0;
@@ -185,7 +183,7 @@ static double refine_iteration(SimplexVertex *s, Deviation *d, int nf) {
 		
 		fom = refine_mean_dev(d, nf, s, i);
 		
-		//printf("Mean deviation at simplex vertex %i = %f nm^-1\n", i, fom/1e9);
+		if ( debug ) printf("Mean deviation at simplex vertex %i = %f nm^-1\n", i, fom/1e9);
 		if ( fom > fom_worst ) {
 			v_second_worst = v_worst;
 			fom_second_worst = fom_worst;
@@ -198,27 +196,28 @@ static double refine_iteration(SimplexVertex *s, Deviation *d, int nf) {
 		}
 		
 	}
+	if ( debug ) printf("The worst vertex is number %i\n", v_worst);
 	
 	/* Reflect this vertex across the opposite face */
 	refine_simplex_transform(s, v_worst, -1.0);
 	
 	/* Is the worst vertex any better? */
 	fom_new = refine_mean_dev(d, nf, s, v_worst);
-//	printf("New mean deviation for the worst vertex after reflection is %f nm^-1\n", fom_new/1e9);
+	if ( debug ) printf("New mean deviation for the worst vertex after reflection is %f nm^-1\n", fom_new/1e9);
 	if ( fom_new > fom_worst ) {
 	
 		double fom_new_new;
 		
 		/* It's worse than before.  Contract in 1D and see if that helps. */
-//		printf("Worse.  Trying a 1D contraction...\n");
+		if ( debug ) printf("Worse.  Trying a 1D contraction...\n");
 		refine_simplex_transform(s, v_worst, 0.5);
 		fom_new_new = refine_mean_dev(d, nf, s, v_worst);
-//		printf("Mean deviation after 1D contraction is %f nm^-1\n", fom_new_new/1e9);
+		if ( debug ) printf("Mean deviation after 1D contraction is %f nm^-1\n", fom_new_new/1e9);
 		if ( fom_new_new > fom_second_worst ) {
 			
 			int i;
 			
-//			printf("Not as good as the second worst vertex: contracting around the best vertex (%i)\n", v_best);
+			if ( debug ) printf("Not as good as the second worst vertex: contracting around the best vertex (%i)\n", v_best);
 			for ( i=0; i<10; i++ ) {
 				if ( i != v_best ) refine_simplex_transform(s, i, 0.5);
 			}
@@ -227,22 +226,24 @@ static double refine_iteration(SimplexVertex *s, Deviation *d, int nf) {
 		
 	} else {
 	
+		#if 0
 		/* It's better.  Try to expand in this direction */
 		double fom_new_new;
 		SimplexVertex save;
 		
-//		printf("This is better.  Trying to expand...\n");
+		if ( debug ) printf("This is better.  Trying to expand...\n");
 		
 		save = s[v_worst];
 		refine_simplex_transform(s, v_worst, 2);
 		/* Better? */
 		fom_new_new = refine_mean_dev(d, nf, s, v_worst);
-//		printf("Mean deviation after expansion is %f nm^-1\n", fom_new_new/1e9);
+		if ( debug ) printf("Mean deviation after expansion is %f nm^-1\n", fom_new_new/1e9);
 		if ( fom_new_new > fom_new ) {
 			/* "Got too confident" */
 			s[v_worst] = save;
-//			printf("Got too confident - reverting\n");
+			if ( debug ) printf("Got too confident - reverting\n");
 		} /* else yay. */
+		#endif
 	
 	}
 	
@@ -257,6 +258,7 @@ double refine_do_cell(ControlContext *ctx) {
 	double delta;
 	int i, nf, f, it, maxiter;
 	const double tol = 0.001e9;	/* Stopping condition */
+	int debug = 1;
 	
 	if ( !ctx->cell_lattice ) {
 		displaywindow_error("No reciprocal unit cell has been found.", ctx->dw);
@@ -282,7 +284,7 @@ double refine_do_cell(ControlContext *ctx) {
 		}
 		
 	}
-	printf("RF: There are %i partnered features in total\n", nf);
+	if ( debug ) printf("RF: There are %i partnered features in total\n", nf);
 	
 	/* Initialise the 'deviation table' */
 	d = malloc(nf*sizeof(Deviation));
@@ -355,13 +357,15 @@ double refine_do_cell(ControlContext *ctx) {
 	
 		double conv;
 		
-//		printf("Simplex method iteration %i\n", it);
-		conv = refine_iteration(s, d, nf);
+		if ( debug ) printf("------------------- Simplex method iteration %i -------------------\n", it);
+		conv = refine_iteration(s, d, nf, debug);
 		if ( conv < tol ) {
-			printf("Converged after %i iterations (%f nm^-1)\n", it, conv/1e9);
+			if ( debug ) printf("RF: Converged after %i iterations (%f nm^-1)\n", it, conv/1e9);
 			break;
 		}
+		
 	}
+	if ( it == maxiter ) printf("RF: Did not converge.\n");
 	
 	/* Apply the final values to the cell */
 	ctx->cell->a.x += s[0].dax;	ctx->cell->b.x += s[0].dbx;	ctx->cell->c.x += s[0].dcx;
diff --git a/src/refinetest.c b/src/refinetest.c
new file mode 100644
index 0000000..f24862c
--- /dev/null
+++ b/src/refinetest.c
@@ -0,0 +1,88 @@
+/*
+ * refinectx.c
+ *
+ * Unit test for refinement procedure
+ *
+ * (c) 2007-2008 Thomas White <taw27@cam.ac.uk>
+ *
+ *  dtr - Diffraction Tomography Reconstruction
+ *
+ */
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include <stdlib.h>
+#include <math.h>
+
+#include "basis.h"
+#include "reflections.h"
+#include "image.h"
+#include "reproject.h"
+#include "refine.h"
+#include "utils.h"
+#include "mapping.h"
+#include "control.h"
+#include "displaywindow.h"
+
+int main(int argc, char *argv[]) {
+
+	ControlContext *ctx;
+	ReflectionList *reflections_real;
+	Basis *cell_real;
+	
+	ctx = control_ctx_new();
+	ctx->omega = deg2rad(0.0);
+	ctx->lambda = lambda(300.0e3);	/* 300 keV */
+	ctx->fmode = FORMULATION_PIXELSIZE;
+	ctx->x_centre = 256;
+	ctx->y_centre = 256;
+	ctx->pixel_size = 5e7;
+	image_add(ctx->images, NULL, 512, 512, deg2rad(0.0), ctx);
+	ctx->images->images[0].features = image_feature_list_new();
+	
+	/* Fudge to avoid horrifying pointer-related death */
+	ctx->dw = malloc(sizeof(DisplayWindow));
+	ctx->dw->cur_image = 0;
+	
+	/* The "true" cell */
+	cell_real = malloc(sizeof(Basis));
+	cell_real->a.x = 5.0e9;  cell_real->a.y = 0.0e9;  cell_real->a.z = 0.0e9;
+	cell_real->b.x = 0.0e9;  cell_real->b.y = 5.0e9;  cell_real->b.z = 0.0e9;
+	cell_real->c.x = 0.0e9;  cell_real->c.y = 0.0e9;  cell_real->c.z = 5.0e9;
+	/* The "real" reflections */
+	reflections_real = reflection_list_from_cell(cell_real);
+	ctx->images->images[0].features = reproject_get_reflections(&ctx->images->images[0], reflections_real);
+	printf("RT: %i test features generated.\n", ctx->images->images[0].features->n_features);
+	
+	/* The "model" cell to be refined */
+	ctx->cell = malloc(sizeof(Basis));
+	ctx->cell->a.x = 5.2e9;  ctx->cell->a.y = 0.1e9;  ctx->cell->a.z = 0.0e9;
+	ctx->cell->b.x = 0.2e9;  ctx->cell->b.y = 4.8e9;  ctx->cell->b.z = 0.0e9;
+	ctx->cell->c.x = 0.0e9;  ctx->cell->c.y = 0.0e9;  ctx->cell->c.z = 5.0e9;
+	ctx->cell_lattice = reflection_list_from_cell(ctx->cell);
+	ctx->images->images[0].rflist = reproject_get_reflections(&ctx->images->images[0], ctx->cell_lattice);
+	reproject_partner_features(ctx->images->images[0].rflist, &ctx->images->images[0]);
+	
+	refine_do_cell(ctx);
+	image_feature_list_free(ctx->images->images[0].rflist);
+	reflection_list_from_new_cell(ctx->cell_lattice, ctx->cell);
+	ctx->images->images[0].rflist = reproject_get_reflections(&ctx->images->images[0], ctx->cell_lattice);
+	
+	free(ctx);
+	
+	return 0;
+	
+}
+
+/* Dummy function stubs */
+#include "gtk-valuegraph.h"
+void displaywindow_update_imagestack(DisplayWindow *dw) { };
+void displaywindow_enable_cell_functions(DisplayWindow *dw, gboolean g) { };
+void displaywindow_update(DisplayWindow *dw) { };
+void displaywindow_error(const char *msg, DisplayWindow *dw) { };
+guint gtk_value_graph_get_type() { return 0; };
+GtkWidget *gtk_value_graph_new() { return NULL; };
+void gtk_value_graph_set_data(GtkValueGraph *vg, double *data, unsigned int n) { };
+
diff --git a/src/reproject.c b/src/reproject.c
index 4036863..e1ce13e 100644
--- a/src/reproject.c
+++ b/src/reproject.c
@@ -3,7 +3,7 @@
  *
  * Synthesize diffraction patterns
  *
- * (c) 2007 Thomas White <taw27@cam.ac.uk>
+ * (c) 2007-2008 Thomas White <taw27@cam.ac.uk>
  *
  *  dtr - Diffraction Tomography Reconstruction
  *
diff --git a/src/reproject.h b/src/reproject.h
index 26d28bb..93f010a 100644
--- a/src/reproject.h
+++ b/src/reproject.h
@@ -3,7 +3,7 @@
  *
  * Synthesize diffraction patterns
  *
- * (c) 2007 Thomas White <taw27@cam.ac.uk>
+ * (c) 2007-2008 Thomas White <taw27@cam.ac.uk>
  *
  *  dtr - Diffraction Tomography Reconstruction
  *