Tidy up refinement unit test

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

2 files changed, 65 insertions, 23 deletions

diff --git a/src/refine.c b/src/refine.c
index e92939d..d0fc8cc 100644
--- a/src/refine.c
+++ b/src/refine.c
@@ -29,8 +29,15 @@
 #include "gtk-valuegraph.h"
 #include "utils.h"
 
+/* Divide numbers by this for display */
 #define DISPFACTOR 1.0e9
 
+/* Number of parameters */
+#define NUM_PARAMS 9
+
+/* Refine debug */
+#define REFINE_DEBUG 0
+
 /* A simplex is an array of ten of these */
 typedef struct {
 	double dax;	double dbx;	double dcx;
@@ -157,7 +164,7 @@ static void refine_simplex_transform(SimplexVertex *s, int v_worst, double fac)
 	centre.dcx = 0.0;	centre.dcy = 0.0;	centre.dcz = 0.0;
 	centre.does_nothing = 0;
 	nv = 0;
-	for ( i=0; i<5; i++ ) {
+	for ( i=0; i<=NUM_PARAMS; i++ ) {
 		if ( i != v_worst ) {
 			centre.dax += s[i].dax;  centre.day += s[i].day;  centre.daz += s[i].daz;
 			centre.dbx += s[i].dbx;  centre.dby += s[i].dby;  centre.dbz += s[i].dbz;
@@ -186,8 +193,10 @@ static void refine_simplex_transform(SimplexVertex *s, int v_worst, double fac)
 	s[v_worst].dcy = centre.dcy + fac * (s[v_worst].dcy - centre.dcy);
 	s[v_worst].dcz = centre.dcz + fac * (s[v_worst].dcz - centre.dcz);
 	
-	printf(" After transformation: ");
-	refine_display_simplex(s[v_worst]);
+	if ( REFINE_DEBUG ) {
+		printf(" After transformation: ");
+		refine_display_simplex(s[v_worst]);
+	}
 	
 }
 
@@ -206,7 +215,7 @@ static void refine_simplex_contract(SimplexVertex *s, int v, int v_best) {
 		
 }
 
-static double refine_iteration(SimplexVertex *s, Deviation *d, int nf, int debug) {
+static double refine_iteration(SimplexVertex *s, Deviation *d, int nf) {
 
 	int v_worst, v_best, v_second_worst, i;
 	double fom_worst, fom_new, fom_best, fom_second_worst;
@@ -218,14 +227,14 @@ static double refine_iteration(SimplexVertex *s, Deviation *d, int nf, int debug
 	fom_best = 100e9;
 	v_second_worst = 0;
 	fom_second_worst = 0.0;
-	if ( debug ) printf("Vertex    FoM/nm^-1    DoesNothing\n");
-	for ( i=0; i<5; i++ ) {
+	if ( REFINE_DEBUG ) printf("Vertex    FoM/nm^-1    DoesNothing\n");
+	for ( i=0; i<=NUM_PARAMS; i++ ) {
 	
 		double fom;
 		
 		fom = refine_mean_dev(d, nf, s, i);
 		
-		if ( debug ) printf("%6i     %8f     %s\n", i, fom/DISPFACTOR, s[i].does_nothing?"*":" ");
+		if ( REFINE_DEBUG ) printf("%6i     %8f     %s\n", i, fom/DISPFACTOR, s[i].does_nothing?"*":" ");
 		if ( (s[i].does_nothing == 0) && (fom > fom_worst) ) {
 			v_second_worst = v_worst;
 			fom_second_worst = fom_worst;
@@ -238,29 +247,29 @@ static double refine_iteration(SimplexVertex *s, Deviation *d, int nf, int debug
 		}
 		
 	}
-	if ( debug ) printf("The worst vertex is number %i\n", v_worst);
+	if ( REFINE_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);
-	if ( debug ) printf("New mean deviation for the worst vertex after reflection is %f nm^-1\n", fom_new/DISPFACTOR);
+	if ( REFINE_DEBUG ) printf("New mean deviation for the worst vertex after reflection is %f nm^-1\n", fom_new/DISPFACTOR);
 	if ( fom_new > fom_worst ) {
 	
 		double fom_new_new;
 		
 		/* It's worse than before.  Contract in 1D and see if that helps. */
-		if ( debug ) printf("Worse.  Trying a 1D contraction...\n");
+		if ( REFINE_DEBUG ) printf("Worse.  Trying a 1D contraction...\n");
 		refine_simplex_transform(s, v_worst, -0.5);	/* Minus puts it back on the original side of the 'good' face */
 		fom_new_new = refine_mean_dev(d, nf, s, v_worst);
-		if ( debug ) printf("Mean deviation after 1D contraction is %f nm^-1\n", fom_new_new/DISPFACTOR);
+		if ( REFINE_DEBUG ) printf("Mean deviation after 1D contraction is %f nm^-1\n", fom_new_new/DISPFACTOR);
 		if ( fom_new_new > fom_second_worst ) {
 			
 			int i;
 			
-			if ( debug ) printf("Not as good as the second worst vertex: contracting around the best vertex (%i)\n", v_best);
-			for ( i=0; i<5; i++ ) {
+			if ( REFINE_DEBUG ) printf("Not as good as the second worst vertex: contracting around the best vertex (%i)\n", v_best);
+			for ( i=0; i<=NUM_PARAMS; i++ ) {
 				if ( i != v_best ) refine_simplex_contract(s, i, v_best);
 			}
 		
@@ -272,19 +281,19 @@ static double refine_iteration(SimplexVertex *s, Deviation *d, int nf, int debug
 		double fom_new_new;
 		SimplexVertex save;
 		
-		if ( debug ) printf("This is better.  Trying to expand...\n");
+		if ( REFINE_DEBUG ) printf("This is better.  Trying to expand...\n");
 		
 		save = s[v_worst];
 		refine_simplex_transform(s, v_worst, 2.0);	/* +ve means stay on this side of the 'good' face */
 		/* Better? */
 		fom_new_new = refine_mean_dev(d, nf, s, v_worst);
-		if ( debug ) printf("Mean deviation after expansion is %f nm^-1\n", fom_new_new/DISPFACTOR);
+		if ( REFINE_DEBUG ) printf("Mean deviation after expansion is %f nm^-1\n", fom_new_new/DISPFACTOR);
 		if ( fom_new_new > fom_new ) {
 			/* "Got too confident" */
 			s[v_worst] = save;
-			if ( debug ) printf("Got too confident - reverting\n");
+			if ( REFINE_DEBUG ) printf("Got too confident - reverting\n");
 		} else {
-			if ( debug ) printf("Better still.  Great.\n");
+			if ( REFINE_DEBUG ) printf("Better still.  Great.\n");
 		}
 	
 	} else {
@@ -294,6 +303,19 @@ static double refine_iteration(SimplexVertex *s, Deviation *d, int nf, int debug
 		
 	}
 	
+	/* Check convergence and return */
+	fom_worst = 0;
+	fom_best = 100e9;
+	for ( i=0; i<=NUM_PARAMS; i++ ) {
+		double fom;
+		fom = refine_mean_dev(d, nf, s, i);
+		if ( (s[i].does_nothing == 0) && (fom > fom_worst) ) {
+			fom_worst = fom;
+		}
+		if ( fom < fom_best ) {
+			fom_best = fom;
+		}
+	}
 	return fom_worst - fom_best;
 
 }
@@ -306,7 +328,6 @@ double refine_do_cell(ControlContext *ctx) {
 	int i, nf, f, it, maxiter;
 	const double tol = 0.0001e9;	/* Stopping condition */
 	//const double tol = 0.001;	/* For testing */
-	int debug = 1;
 	
 	if ( !ctx->cell_lattice ) {
 		displaywindow_error("No reciprocal unit cell has been found.", ctx->dw);
@@ -332,7 +353,7 @@ double refine_do_cell(ControlContext *ctx) {
 		}
 		
 	}
-	if ( debug ) printf("RF: There are %i partnered features in total\n", nf);
+	if ( REFINE_DEBUG ) printf("RF: There are %i partnered features in total\n", nf);
 	
 	/* Initialise the 'deviation table' */
 	d = malloc(nf*sizeof(Deviation));
@@ -410,10 +431,10 @@ double refine_do_cell(ControlContext *ctx) {
 		//	refine_display_simplex(s[i]);
 		//}
 		
-		if ( debug ) printf("------------------- Simplex method iteration %i -------------------\n", it);
-		conv = refine_iteration(s, d, nf, debug);
+		if ( REFINE_DEBUG ) printf("------------------- Simplex method iteration %i -------------------\n", it);
+		conv = refine_iteration(s, d, nf);
 		if ( conv < tol ) {
-			if ( debug ) printf("RF: Converged after %i iterations (%f nm^-1)\n", it, conv/DISPFACTOR);
+			if ( REFINE_DEBUG ) printf("RF: Converged after %i iterations (%f nm^-1)\n", it, conv/DISPFACTOR);
 			break;
 		}
 		
diff --git a/src/refinetest.c b/src/refinetest.c
index f24862c..10e0b5f 100644
--- a/src/refinetest.c
+++ b/src/refinetest.c
@@ -31,6 +31,7 @@ int main(int argc, char *argv[]) {
 	ControlContext *ctx;
 	ReflectionList *reflections_real;
 	Basis *cell_real;
+	int fail;
 	
 	ctx = control_ctx_new();
 	ctx->omega = deg2rad(0.0);
@@ -54,7 +55,7 @@ int main(int argc, char *argv[]) {
 	/* 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);
+//	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));
@@ -70,8 +71,28 @@ int main(int argc, char *argv[]) {
 	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);
 	
+	fail = 0;
+	if ( fabs(ctx->cell->a.x - cell_real->a.x) > 0.01e9 ) {
+		fprintf(stderr, "refinetest: ax not inferred correctly\n");
+		fail = 1;
+	}
+	if ( fabs(ctx->cell->a.y - cell_real->a.y) > 0.01e9 ) {
+		fprintf(stderr, "refinetest: ay not inferred correctly\n");
+		fail = 1;
+	}
+	if ( fabs(ctx->cell->b.x - cell_real->b.x) > 0.01e9 ) {
+		fprintf(stderr, "refinetest: bx not inferred correctly\n");
+		fail = 1;
+	}
+	if ( fabs(ctx->cell->b.y - cell_real->b.y) > 0.01e9 ) {
+		fprintf(stderr, "refinetest: by not inferred correctly\n");
+		fail = 1;
+	}
+	
 	free(ctx);
 	
+	if ( fail ) return 1;
+	
 	return 0;
 	
 }