Sequencer etc

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

2 files changed, 54 insertions, 24 deletions

diff --git a/src/refine.c b/src/refine.c
index 4b6793e..4d2a12e 100644
--- a/src/refine.c
+++ b/src/refine.c
@@ -25,6 +25,7 @@
 #include "image.h"
 #include "reproject.h"
 #include "mapping.h"
+#include "refine.h"
 
 /* A simplex is an array of ten of these */
 typedef struct {
@@ -40,7 +41,26 @@ typedef struct {
 
 void refine_do_sequence(ControlContext *ctx) {
 
-	
+	double omega_offs;
+
+	for ( omega_offs=-2.0; omega_offs<=2.0; omega_offs+=0.01 ) {
+
+		double fit;
+		int i;
+		
+		for ( i=0; i<ctx->images->n_images; i++ ) {
+			ctx->images->images[i].omega += omega_offs;
+		}
+		
+		fit = refine_do_cell(ctx);
+		printf("RF:                           omega_offs=%f, fit=%f nm^-1\n", omega_offs, fit/1e9);
+		
+		for ( i=0; i<ctx->images->n_images; i++ ) {
+			ctx->images->images[i].omega -= omega_offs;
+		}
+
+	}
+
 }
 
 static double refine_mean_dev(Deviation *d, int nf, SimplexVertex *s, int i) {
@@ -67,6 +87,7 @@ static double refine_mean_dev(Deviation *d, int nf, SimplexVertex *s, int i) {
 
 }
 
+#if 0
 static void refine_display_simplex(SimplexVertex *s) {
 
 	int i;
@@ -79,6 +100,7 @@ static void refine_display_simplex(SimplexVertex *s) {
 	}
 	
 }
+#endif
 
 /* Expand the simplex across from vertex v_worst by factor 'fac'.
  *	fac = -1 is a reflection
@@ -103,14 +125,13 @@ static void refine_simplex_transform(SimplexVertex *s, int v_worst, double fac)
 	
 }
 
-static void refine_iteration(SimplexVertex *s, Deviation *d, int nf) {
+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;
 
-	/* Start an iteration */
-	refine_display_simplex(s);
-	
+//	printf("---\n");
+
 	/* Find the least favourable vertex of the simplex */
 	v_worst = 0;
 	fom_worst = 0.0;
@@ -124,43 +145,40 @@ static void refine_iteration(SimplexVertex *s, Deviation *d, int nf) {
 		
 		fom = refine_mean_dev(d, nf, s, i);
 		
-		printf("Mean deviation at simplex vertex %i = %f nm^-1\n", i, fom/1e9);
+		//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;
 		}
-		if ( fom > fom_best ) {
+		if ( fom < fom_best ) {
 			fom_best = fom;
 			v_best = i;
 		}
 		
 	}
-	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.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);
+//	printf("New mean deviation for the worst vertex after reflection is %f nm^-1\n", fom_new/1e9);
 	if ( fom_new > fom_worst ) {
 	
 		double fom_new_new;
 		
 		/* It's worse than before.  Contract in 1D and see if that helps. */
-		printf("Worse.  Trying a 1D contraction...\n");
+//		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);
+//		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);
+//			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);
 			}
@@ -173,38 +191,41 @@ static void refine_iteration(SimplexVertex *s, Deviation *d, int nf) {
 		double fom_new_new;
 		SimplexVertex save;
 		
-		printf("This is better.  Trying to expand...\n");
+//		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);
+//		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");
+//			printf("Got too confident - reverting\n");
 		} /* else yay. */
 	
 	}
+	
+	return fom_worst - fom_best;
 
 }
 
-void refine_do_cell(ControlContext *ctx) {
+double refine_do_cell(ControlContext *ctx) {
 
 	SimplexVertex s[10];
 	Deviation *d;
 	double delta;
 	int i, nf, f, it, maxiter;
+	const double tol = 0.001e9;	/* Stopping condition */
 	
 	if ( !ctx->cell_lattice ) {
 		displaywindow_error("No reciprocal unit cell has been found.", ctx->dw);
-		return;
+		return -1;
 	}
 	
 	if ( ctx->images->n_images == 0 ) {
 		displaywindow_error("There are no images to refine against.", ctx->dw);
-		return;
+		return -1;
 	}
 	
 	/* Determine the size of the 'deviation table' */
@@ -291,13 +312,22 @@ void refine_do_cell(ControlContext *ctx) {
 	/* Iterate */
 	maxiter = 500;
 	for ( it=0; it<maxiter; it++ ) {
-		printf("Simplex method iteration %i\n", it);
-		refine_iteration(s, d, nf);
+	
+		double conv;
+		
+//		printf("Simplex method iteration %i\n", it);
+		conv = refine_iteration(s, d, nf);
+		if ( conv < tol ) {
+			printf("Converged after %i iterations (%f)\n", it, conv);
+			break;
+		}
 	}
 	
 	ctx->images->images[ctx->dw->cur_image].rflist = NULL;
 	reproject_lattice_changed(ctx);
 	displaywindow_update(ctx->dw);
+	
+	return refine_mean_dev(d, nf, s, 0);
 		
 }
 
diff --git a/src/refine.h b/src/refine.h
index d3ca6d9..38dee22 100644
--- a/src/refine.h
+++ b/src/refine.h
@@ -19,7 +19,7 @@
 #include "control.h"
 
 extern void refine_do_sequence(ControlContext *ctx);
-extern void refine_do_cell(ControlContext *ctx);
+extern double refine_do_cell(ControlContext *ctx);
 
 #endif	/* REFINE_H */