Refinement stuff

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

1 files changed, 120 insertions, 30 deletions

diff --git a/src/refine.c b/src/refine.c
index 7871dc1..6eff91d 100644
--- a/src/refine.c
+++ b/src/refine.c
@@ -63,11 +63,115 @@ static double refine_mean_dev(Deviation *d, int nf, SimplexVertex *s, int i) {
 		
 	}
 	
-	return fom;
+	return fom/nf;
 
 }
 
-static void refine_simplex_reflect(SimplexVertex *s, int v_worst) {
+static void refine_display_simplex(SimplexVertex *s) {
+
+	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);
+	}
+	
+}
+
+/* Expand the simplex across from vertex v_worst by factor 'fac'.
+ *	fac = -1 is a reflection
+ *	fac = +n is a 1d expansion
+ */
+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 */
+	
+	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);
+	
+}
+
+static void refine_iteration(SimplexVertex *s, Deviation *d, int nf) {
+
+	int v_worst, v_best, i;
+	double fom_worst, fom_new, fom_best;
+
+	/* Start an iteration */
+	refine_display_simplex(s);
+	
+	/* Find the least favourable vertex of the simplex */
+	v_worst = 0;
+	fom_worst = 0.0;
+	v_best = 0;
+	fom_best = 100e9;
+	for ( i=0; i<10; i++ ) {
+	
+		double fom;
+		
+		fom = refine_mean_dev(d, nf, s, i);
+		
+		printf("Mean deviation at simplex vertex %i = %f nm^-1\n", i, fom/1e9);
+		if ( fom > fom_worst ) {
+			fom_worst = fom;
+			v_worst = i;
+		}
+		if ( fom > fom_best ) {
+			fom_best = fom;
+			v_best = i;
+		}
+		
+	}
+	printf("The worst vertex is number %i\n", v_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);
+		
+	/* 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);
+		}
+		
+	} else {
+	
+		/* It's better.  Try to expand in this direction */
+		double fom_new_new;
+		SimplexVertex save;
+		
+		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 ( fom_new_new > fom_new ) {
+			/* "Got too confident" */
+			s[v_worst] = save;
+			printf("Got too confident - reverting\n");
+		} /* else yay. */
+	
+	}
+
 }
 
 void refine_do_cell(ControlContext *ctx) {
@@ -75,8 +179,7 @@ void refine_do_cell(ControlContext *ctx) {
 	SimplexVertex s[10];
 	Deviation *d;
 	double delta;
-	int i, nf, f, v_worst;
-	double fom_worst;
+	int i, nf, f, it, maxiter;
 	
 	if ( !ctx->cell_lattice ) {
 		displaywindow_error("No reciprocal unit cell has been found.", ctx->dw);
@@ -88,7 +191,7 @@ void refine_do_cell(ControlContext *ctx) {
 		return;
 	}
 	
-	/* Create the table of indicies and deviations */
+	/* Determine the size of the 'deviation table' */
 	nf = 0;
 	for ( i=0; i<ctx->images->n_images; i++ ) {
 		int j;
@@ -104,6 +207,7 @@ void refine_do_cell(ControlContext *ctx) {
 	}
 	printf("RF: There are %i partnered features in total\n", nf);
 	
+	/* Initialise the 'deviation table' */
 	d = malloc(nf*sizeof(Deviation));
 	f = 0;
 	for ( i=0; i<ctx->images->n_images; i++ ) {
@@ -159,35 +263,21 @@ void refine_do_cell(ControlContext *ctx) {
 	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;
 	memcpy(&s[1], &s[0], sizeof(SimplexVertex));	s[1].dax = delta;
-	memcpy(&s[2], &s[0], sizeof(SimplexVertex));	s[2].dbx = delta;
-	memcpy(&s[3], &s[0], sizeof(SimplexVertex));	s[3].dcx = delta;
-	memcpy(&s[4], &s[0], sizeof(SimplexVertex));	s[4].day = delta;
+	memcpy(&s[2], &s[0], sizeof(SimplexVertex));	s[2].day = delta;
+	memcpy(&s[3], &s[0], sizeof(SimplexVertex));	s[3].daz = delta;
+	memcpy(&s[4], &s[0], sizeof(SimplexVertex));	s[4].dbx = delta;
 	memcpy(&s[5], &s[0], sizeof(SimplexVertex));	s[5].dby = delta;
-	memcpy(&s[6], &s[0], sizeof(SimplexVertex));	s[6].dcy = delta;
-	memcpy(&s[7], &s[0], sizeof(SimplexVertex));	s[7].daz = delta;
-	memcpy(&s[8], &s[0], sizeof(SimplexVertex));	s[8].dbz = delta;
+	memcpy(&s[6], &s[0], sizeof(SimplexVertex));	s[6].dbz = delta;
+	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;
 	
-	/* Find the least favourable vertex of the simplex */
-	v_worst = 0;
-	fom_worst = 0.0;
-	for ( i=0; i<10; i++ ) {
-	
-		double fom;
-		
-		fom = refine_mean_dev(d, nf, s, i);
-		
-		printf("Mean deviation at simplex vertex %i = %f nm^-1\n", i, (fom/nf)/1e9);
-		if ( fom > fom_worst ) {
-			fom_worst = fom;
-			v_worst = i;
-		}
-		
+	/* Iterate */
+	maxiter = 2;
+	for ( it=0; it<maxiter; it++ ) {
+		printf("Simplex method iteration %i\n", it);
+		refine_iteration(s, d, nf);
 	}
-	printf("The worst vertex is number %i\n", v_worst);
-	
-	/* Reflect this vertex across the opposite face */
-	refine_simplex_reflect(s, v_worst);
 	
 	ctx->images->images[ctx->dw->cur_image].rflist = NULL;
 	reproject_lattice_changed(ctx);