Make refinement work?

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

1 files changed, 34 insertions, 18 deletions

diff --git a/src/refine.c b/src/refine.c
index 6eff91d..4b6793e 100644
--- a/src/refine.c
+++ b/src/refine.c
@@ -91,22 +91,22 @@ static void refine_simplex_transform(SimplexVertex *s, int v_worst, double fac)
 	v_face = 0;
 	if ( v_worst == 0 ) v_face = 1;
 
-	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);
+	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);
 	
 }
 
 static void refine_iteration(SimplexVertex *s, Deviation *d, int nf) {
 
-	int v_worst, v_best, i;
-	double fom_worst, fom_new, fom_best;
+	int v_worst, v_best, v_second_worst, i;
+	double fom_worst, fom_new, fom_best, fom_second_worst;
 
 	/* Start an iteration */
 	refine_display_simplex(s);
@@ -116,6 +116,8 @@ static void refine_iteration(SimplexVertex *s, Deviation *d, int nf) {
 	fom_worst = 0.0;
 	v_best = 0;
 	fom_best = 100e9;
+	v_second_worst = 0;
+	fom_second_worst = 0.0;
 	for ( i=0; i<10; i++ ) {
 	
 		double fom;
@@ -124,6 +126,8 @@ static void refine_iteration(SimplexVertex *s, Deviation *d, int nf) {
 		
 		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;
 			fom_worst = fom;
 			v_worst = i;
 		}
@@ -134,21 +138,33 @@ static void refine_iteration(SimplexVertex *s, Deviation *d, int nf) {
 		
 	}
 	printf("The worst vertex is number %i\n", v_worst);
+	printf("The second worst vertex is number %i\n", v_second_worst);
 	printf("The best vertex is number %i\n", v_best);
 	
 	/* Reflect this vertex across the opposite face */
-	refine_simplex_transform(s, v_worst, -1);
+	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 ( fom_new > fom_worst ) {
 	
-		/* It's worse than before.  Contract around the best vertex */
-		int i;
-		printf("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);
+		double fom_new_new;
+		
+		/* It's worse than before.  Contract in 1D and see if that helps. */
+		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 ( 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);
+			for ( i=0; i<10; i++ ) {
+				if ( i != v_best ) refine_simplex_transform(s, i, 0.5);
+			}
+		
 		}
 		
 	} else {
@@ -273,7 +289,7 @@ void refine_do_cell(ControlContext *ctx) {
 	memcpy(&s[9], &s[0], sizeof(SimplexVertex));	s[9].dcz = delta;
 	
 	/* Iterate */
-	maxiter = 2;
+	maxiter = 500;
 	for ( it=0; it<maxiter; it++ ) {
 		printf("Simplex method iteration %i\n", it);
 		refine_iteration(s, d, nf);