Refinement progress

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

1 files changed, 80 insertions, 32 deletions

diff --git a/src/refine.c b/src/refine.c
index 7bb2585..ac2ba38 100644
--- a/src/refine.c
+++ b/src/refine.c
@@ -34,6 +34,7 @@ typedef struct {
 	double dax;	double dbx;	double dcx;
 	double day;	double dby;	double dcy;
 	double daz;	double dbz;	double dcz;
+	int does_nothing;	/* If non-zero, this vertex doesn't appear to affect the FoM */
 } SimplexVertex;
 
 typedef struct {
@@ -127,20 +128,14 @@ static double refine_mean_dev(Deviation *d, int nf, SimplexVertex *s, int i) {
 
 }
 
-#if 0
-static void refine_display_simplex(SimplexVertex *s) {
+static void refine_display_simplex(SimplexVertex sv) {
 
-	int i;
-	
-	for ( i=0; i<10; i++ ) {
-		printf("Vertex %i: %8f %8f %8f %8f %8f %8f %8f %8f %8f\n",
-				i, s[i].dax/1e9, s[i].day/1e9, s[i].daz/1e9,
-				   s[i].dbx/1e9, s[i].dby/1e9, s[i].dbz/1e9,
-				   s[i].dcx/1e9, s[i].dcy/1e9, s[i].dcz/1e9);
-	}
+	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);
 	
 }
-#endif
 
 /* Expand the simplex across from vertex v_worst by factor 'fac'.
  *	fac = -1 is a reflection
@@ -148,21 +143,62 @@ static void refine_display_simplex(SimplexVertex *s) {
  */
 static void refine_simplex_transform(SimplexVertex *s, int v_worst, double fac) {
 
-	int v_face;	/* A simplex vertex which is on the face being reflected across */
+	SimplexVertex centre;
+	int i, nv;
+	
+	/* Average the coordinates of the non-worst vertices to
+	 *  get the centre of the opposite face. */
+	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;
+	nv = 0;
+	for ( i=0; i<10; 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;
+			nv++;
+		}
+	}
+	centre.dax /= nv;	centre.day /= nv;	centre.daz /= nv;
+	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]);
 	
-	v_face = 0;
-	if ( v_worst == 0 ) v_face = 1;
+	/* Do the transformation */
+	s[v_worst].dax = centre.dax + fac * (s[v_worst].dax - centre.dax);
+	s[v_worst].day = centre.day + fac * (s[v_worst].day - centre.day);
+	s[v_worst].daz = centre.daz + fac * (s[v_worst].daz - centre.daz);
+	s[v_worst].dbx = centre.dbx + fac * (s[v_worst].dbx - centre.dbx);
+	s[v_worst].dby = centre.dby + fac * (s[v_worst].dby - centre.dby);
+	s[v_worst].dbz = centre.dbz + fac * (s[v_worst].dbz - centre.dbz);
+	s[v_worst].dcx = centre.dcx + fac * (s[v_worst].dcx - centre.dcx);
+	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]);
+	
+}
 
-	s[v_worst].dax = s[v_face].dax + fac * (s[v_worst].dax - s[v_face].dax);
-	s[v_worst].day = s[v_face].day + fac * (s[v_worst].day - s[v_face].dax);
-	s[v_worst].daz = s[v_face].daz + fac * (s[v_worst].daz - s[v_face].dax);
-	s[v_worst].dbx = s[v_face].dbx + fac * (s[v_worst].dbx - s[v_face].dbx);
-	s[v_worst].dby = s[v_face].dby + fac * (s[v_worst].dby - s[v_face].dbx);
-	s[v_worst].dbz = s[v_face].dbz + fac * (s[v_worst].dbz - s[v_face].dbx);
-	s[v_worst].dcx = s[v_face].dcx + fac * (s[v_worst].dcx - s[v_face].dcx);
-	s[v_worst].dcy = s[v_face].dcy + fac * (s[v_worst].dcy - s[v_face].dcx);
-	s[v_worst].dcz = s[v_face].dcz + fac * (s[v_worst].dcz - s[v_face].dcx);
+/* Contract vertex v of simplex s towards vertex v_best */
+static void refine_simplex_contract(SimplexVertex *s, int v, int v_best) {
 	
+	s[v].dax = s[v_best].dax + 0.5 * (s[v].dax - s[v_best].dax);
+	s[v].day = s[v_best].day + 0.5 * (s[v].day - s[v_best].day);
+	s[v].daz = s[v_best].daz + 0.5 * (s[v].daz - s[v_best].daz);
+	s[v].dbx = s[v_best].dbx + 0.5 * (s[v].dbx - s[v_best].dbx);
+	s[v].dby = s[v_best].dby + 0.5 * (s[v].dby - s[v_best].dby);
+	s[v].dbz = s[v_best].dbz + 0.5 * (s[v].dbz - s[v_best].dbz);
+	s[v].dcx = s[v_best].dcx + 0.5 * (s[v].dcx - s[v_best].dcx);
+	s[v].dcy = s[v_best].dcy + 0.5 * (s[v].dcy - s[v_best].dcy);
+	s[v].dcz = s[v_best].dcz + 0.5 * (s[v].dcz - s[v_best].dcz);
+		
 }
 
 static double refine_iteration(SimplexVertex *s, Deviation *d, int nf, int debug) {
@@ -177,14 +213,15 @@ 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<10; i++ ) {
 	
 		double fom;
 		
 		fom = refine_mean_dev(d, nf, s, i);
 		
-		if ( debug ) printf("Mean deviation at simplex vertex %i = %f nm^-1\n", i, fom/1e9);
-		if ( fom > fom_worst ) {
+		if ( debug ) printf("%6i     %8f     %s\n", i, fom/1e9, s[i].does_nothing?"*":" ");
+		if ( (s[i].does_nothing == 0) && (fom > fom_worst) ) {
 			v_second_worst = v_worst;
 			fom_second_worst = fom_worst;
 			fom_worst = fom;
@@ -210,7 +247,7 @@ static double refine_iteration(SimplexVertex *s, Deviation *d, int nf, int debug
 		
 		/* It's worse than before.  Contract in 1D and see if that helps. */
 		if ( debug ) printf("Worse.  Trying a 1D contraction...\n");
-		refine_simplex_transform(s, v_worst, 0.5);
+		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 ( fom_new_new > fom_second_worst ) {
@@ -219,14 +256,13 @@ static double refine_iteration(SimplexVertex *s, Deviation *d, int nf, int debug
 			
 			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);
+				if ( i != v_best ) refine_simplex_contract(s, i, v_best);
 			}
 		
 		}
 		
-	} else {
+	} else if ( fom_new < fom_worst ) {
 	
-		#if 0
 		/* It's better.  Try to expand in this direction */
 		double fom_new_new;
 		SimplexVertex save;
@@ -234,7 +270,7 @@ static double refine_iteration(SimplexVertex *s, Deviation *d, int nf, int debug
 		if ( debug ) printf("This is better.  Trying to expand...\n");
 		
 		save = s[v_worst];
-		refine_simplex_transform(s, v_worst, 2);
+		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);
@@ -242,9 +278,15 @@ static double refine_iteration(SimplexVertex *s, Deviation *d, int nf, int debug
 			/* "Got too confident" */
 			s[v_worst] = save;
 			if ( debug ) printf("Got too confident - reverting\n");
-		} /* else yay. */
-		#endif
+		} else {
+			if ( debug ) printf("Better still.  Great.\n");
+		}
+	
+	} else {
 	
+		printf("No change!\n");
+		s[v_worst].does_nothing = 1;
+		
 	}
 	
 	return fom_worst - fom_best;
@@ -341,6 +383,7 @@ double refine_do_cell(ControlContext *ctx) {
 	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;
+	s[0].does_nothing = 0;
 	memcpy(&s[1], &s[0], sizeof(SimplexVertex));	s[1].dax = delta;
 	memcpy(&s[2], &s[0], sizeof(SimplexVertex));	s[2].day = delta;
 	memcpy(&s[3], &s[0], sizeof(SimplexVertex));	s[3].daz = delta;
@@ -350,6 +393,11 @@ double refine_do_cell(ControlContext *ctx) {
 	memcpy(&s[7], &s[0], sizeof(SimplexVertex));	s[7].dcx = delta;
 	memcpy(&s[8], &s[0], sizeof(SimplexVertex));	s[8].dcy = delta;
 	memcpy(&s[9], &s[0], sizeof(SimplexVertex));	s[9].dcz = delta;
+	s[3].does_nothing = 1;
+	s[6].does_nothing = 1;
+	s[7].does_nothing = 1;
+	s[8].does_nothing = 1;
+	s[9].does_nothing = 1;
 	
 	/* Iterate */
 	maxiter = 500;