Test minimisation works (yay!)

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

1 files changed, 33 insertions, 25 deletions

diff --git a/src/refine.c b/src/refine.c
index 45402d1..fc983ab 100644
--- a/src/refine.c
+++ b/src/refine.c
@@ -29,6 +29,8 @@
 #include "gtk-valuegraph.h"
 #include "utils.h"
 
+#define DISPFACTOR 1.0
+
 /* A simplex is an array of ten of these */
 typedef struct {
 	double dax;	double dbx;	double dcx;
@@ -71,7 +73,7 @@ void refine_do_sequence(ControlContext *ctx) {
 		reproject_lattice_changed(ctx);
 		
 		fit = refine_do_cell(ctx);
-		printf("RF:                           omega_offs=%f deg, fit=%f nm^-1\n", rad2deg(omega_offs), fit/1e9);
+		printf("RF:                           omega_offs=%f deg, fit=%f nm^-1\n", rad2deg(omega_offs), fit/DISPFACTOR);
 		fit_vals[idx++] = fit;
 		if ( fit < fit_best ) {
 			fit_best = fit;
@@ -94,7 +96,7 @@ void refine_do_sequence(ControlContext *ctx) {
 	gtk_widget_show_all(window_fit);
 	
 	/* Perform final refinement */
-	printf("Best omega offset = %f deg (%f nm^-1)\n", rad2deg(omega_offs_best), fit_best/1e9);
+	printf("Best omega offset = %f deg (%f nm^-1)\n", rad2deg(omega_offs_best), fit_best/DISPFACTOR);
 	for ( j=0; j<ctx->images->n_images; j++ ) {
 		ctx->images->images[j].omega += omega_offs_best;
 	}
@@ -124,18 +126,18 @@ static double refine_mean_dev(Deviation *d, int nf, SimplexVertex *s, int i) {
 		
 	}
 	
-//	return fabs(s[i].dax-10.0) + fabs(s[i].day-20.0) + fabs(s[i].dbx-30.0) + fabs(s[i].dby-40.0);
+	return fabs(s[i].dax-10.0) + fabs(s[i].day-20.0) + fabs(s[i].dbx-30.0) + fabs(s[i].dby-40.0);
 	
-	return fom/nf;
+//	return fom/nf;
 
 }
 
 static void refine_display_simplex(SimplexVertex sv) {
 
-	printf("%8f %8f %8f %8f %8f %8f %8f %8f %8f\n",
-				sv.dax/1e9, sv.day/1e9, sv.daz/1e9,
-				sv.dbx/1e9, sv.dby/1e9, sv.dbz/1e9,
-				sv.dcx/1e9, sv.dcy/1e9, sv.dcz/1e9);
+	printf("%9.6f %9.6f %9.6f %9.6f %9.6f %9.6f %9.6f %9.6f %9.6f\n",
+				sv.dax/DISPFACTOR, sv.day/DISPFACTOR, sv.daz/DISPFACTOR,
+				sv.dbx/DISPFACTOR, sv.dby/DISPFACTOR, sv.dbz/DISPFACTOR,
+				sv.dcx/DISPFACTOR, sv.dcy/DISPFACTOR, sv.dcz/DISPFACTOR);
 	
 }
 
@@ -153,12 +155,13 @@ static void refine_simplex_transform(SimplexVertex *s, int v_worst, double fac)
 	centre.dax = 0.0;	centre.day = 0.0;	centre.daz = 0.0;
 	centre.dbx = 0.0;	centre.dby = 0.0;	centre.dbz = 0.0;
 	centre.dcx = 0.0;	centre.dcy = 0.0;	centre.dcz = 0.0;
+	centre.does_nothing = 0;
 	nv = 0;
 	for ( i=0; i<5; i++ ) {
 		if ( i != v_worst ) {
-			centre.dax += s[i].dax;  centre.dax += s[i].day;  centre.dax += s[i].daz;
-			centre.dax += s[i].dbx;  centre.dax += s[i].dby;  centre.dax += s[i].dbz;
-			centre.dax += s[i].dcx;  centre.dax += s[i].dcy;  centre.dax += s[i].dcz;
+			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;
+			centre.dcx += s[i].dcx;  centre.dcy += s[i].dcy;  centre.dcz += s[i].dcz;
 			nv++;
 		}
 	}
@@ -166,11 +169,11 @@ static void refine_simplex_transform(SimplexVertex *s, int v_worst, double fac)
 	centre.dbx /= nv;	centre.dby /= nv;	centre.dbz /= nv;
 	centre.dcx /= nv;	centre.dcy /= nv;	centre.dcz /= nv;
 	
-//	printf("             Midpoint: ");
-//	refine_display_simplex(centre);
+	printf("Before transformation: ");
+	refine_display_simplex(s[v_worst]);
 	
-//	printf("Before transformation: ");
-//	refine_display_simplex(s[v_worst]);
+	printf("             Midpoint: ");
+	refine_display_simplex(centre);
 	
 	/* Do the transformation */
 	s[v_worst].dax = centre.dax + fac * (s[v_worst].dax - centre.dax);
@@ -183,8 +186,8 @@ 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]);
+	printf(" After transformation: ");
+	refine_display_simplex(s[v_worst]);
 	
 }
 
@@ -222,7 +225,7 @@ static double refine_iteration(SimplexVertex *s, Deviation *d, int nf, int debug
 		
 		fom = refine_mean_dev(d, nf, s, i);
 		
-		if ( debug ) printf("%6i     %8f     %s\n", i, fom/1e9, s[i].does_nothing?"*":" ");
+		if ( 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;
@@ -242,7 +245,7 @@ static double refine_iteration(SimplexVertex *s, Deviation *d, int nf, int debug
 	
 	/* 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/1e9);
+	if ( 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;
@@ -251,7 +254,7 @@ static double refine_iteration(SimplexVertex *s, Deviation *d, int nf, int debug
 		if ( 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/1e9);
+		if ( 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;
@@ -275,7 +278,7 @@ static double refine_iteration(SimplexVertex *s, Deviation *d, int nf, int debug
 		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/1e9);
+		if ( 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;
@@ -301,7 +304,8 @@ double refine_do_cell(ControlContext *ctx) {
 	Deviation *d;
 	double delta;
 	int i, nf, f, it, maxiter;
-	const double tol = 0.001e9;	/* Stopping condition */
+	//const double tol = 0.001e9;	/* Stopping condition */
+	const double tol = 0.001;	/* For testing */
 	int debug = 1;
 	
 	if ( !ctx->cell_lattice ) {
@@ -367,7 +371,7 @@ double refine_do_cell(ControlContext *ctx) {
 			mapping_rotate(dx, dy, 0.0, &dlx, &dly, &dlz, image->omega, image->tilt);
 			rf->partner->x = old_x;
 			rf->partner->y = old_y;
-			//printf("=> %+10.5f %+10.5f %+10.5f nm^-1\n", dlx/1e9, dly/1e9, dlz/1e9);
+			//printf("=> %+10.5f %+10.5f %+10.5f nm^-1\n", dlx/DISPFACTOR, dly/DISPFACTOR, dlz/DISPFACTOR);
 			
 			d[f].dx = dlx;
 			d[f].dy = dly;
@@ -381,7 +385,7 @@ double refine_do_cell(ControlContext *ctx) {
 	assert( f == nf );
 	
 	/* Initialise the simplex */
-	delta = 10.0;//0.01e9;
+	delta = 5.0;//0.01e9;
 	s[0].dax = 0.0;	s[0].dbx = 0.0;	s[0].dcx = 0.0;
 	s[0].day = 0.0;	s[0].dby = 0.0;	s[0].dcy = 0.0;
 	s[0].daz = 0.0;	s[0].dbz = 0.0;	s[0].dcz = 0.0;
@@ -402,10 +406,14 @@ double refine_do_cell(ControlContext *ctx) {
 	
 		double conv;
 		
+		for ( i=0; i<5; i++ ) {
+			refine_display_simplex(s[i]);
+		}
+		
 		if ( debug ) printf("------------------- Simplex method iteration %i -------------------\n", it);
 		conv = refine_iteration(s, d, nf, debug);
 		if ( conv < tol ) {
-			if ( debug ) printf("RF: 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/DISPFACTOR);
 			break;
 		}