Add eFOM (needs debugging)

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

1 files changed, 83 insertions, 9 deletions

diff --git a/src/ipr.c b/src/ipr.c
index 2d78ca1..982d8c3 100644
--- a/src/ipr.c
+++ b/src/ipr.c
@@ -125,9 +125,7 @@ static void ipr_try_compact(Basis *basis) {
 
 }
 
-/* Choose an initial set of three vectors for the expansion.
- * This would probably be just as easy for the user to do... */
-static Basis *ipr_choose_initial_basis(ControlContext *ctx) {
+static Basis *ipr_choose_random_basis(ControlContext *ctx) {
 	
 	Basis *basis;
 	Reflection *reflection;
@@ -135,7 +133,6 @@ static Basis *ipr_choose_initial_basis(ControlContext *ctx) {
 	
 	basis = malloc(sizeof(Basis));
 	
-	/* Determine a vaguely sensible starting basis */
 	basis->a.x = 0;	basis->a.y = 0;	basis->a.z = 0;	basis->a.modulus = 1e30;
 	basis->b.x = 0;	basis->b.y = 0;	basis->b.z = 0;	basis->b.modulus = 1e30;
 	basis->c.x = 0;	basis->c.y = 0;	basis->c.z = 0;	basis->c.modulus = 1e30;
@@ -191,7 +188,7 @@ static Basis *ipr_choose_initial_basis(ControlContext *ctx) {
 		/* End of list? */
 		if ( !reflection ) {
 			if ( (basis->a.modulus == 1e30) || (basis->b.modulus == 1e30) || (basis->c.modulus == 1e30) ) {
-				printf("IP: Failed to find appropriate starting basis!\n");
+				//printf("IP: Failed to find appropriate starting basis!\n");
 				free(basis);
 				return NULL;
 			}
@@ -201,17 +198,94 @@ static Basis *ipr_choose_initial_basis(ControlContext *ctx) {
 	
 	ipr_try_compact(basis);
 	
-	basis->a.x = 1.842e9;  basis->a.y = 0.0;  basis->a.z = 0.0;
-	basis->b.x = 0.0;  basis->b.y = 1.842e9;  basis->b.z = 0.0;
-	basis->c.x = 0.0;  basis->c.y = 0.0;  basis->c.z = 1.842e9;
+//	basis->a.x = 1.842e9;  basis->a.y = 0.0;  basis->a.z = 0.0;
+//	basis->b.x = 0.0;  basis->b.y = 1.842e9;  basis->b.z = 0.0;
+//	basis->c.x = 0.0;  basis->c.y = 0.0;  basis->c.z = 1.842e9;
+	
+	return basis;
+
+}
+
+static double ipr_efom(ReflectionContext *rctx, Basis *basis) {
+
+	int n_indexed, n_counted;
+	Reflection *cur;
+	
+	cur = rctx->reflections;
+	n_indexed = 0;
+	n_counted = 0;
+	while ( cur ) {
+
+		if ( cur->type == REFLECTION_NORMAL ) {
+			
+			/* Can this basis approximately account for this reflection? */
+			double n;
+			double p, q, r, u, v, w, d, e, f;
+			double h, k, l;
+			
+			/* Set up the coordinate transform from hkl to xyz */
+			p = basis->a.x;  q = basis->a.y;  r = basis->a.z;
+			u = basis->b.x;  v = basis->b.y;  w = basis->b.z;
+			d = basis->c.x;  e = basis->c.y;  f = basis->c.z;
+			
+			/* Invert the matrix to get hkl from xyz */
+			n = p*(v*f - w*e) - q*(u*f - w*d) + r*(u*e - d*v);
+			h = (p*(v*f-w*e)*cur->x + -q*(u*f-w*d)*cur->y + r*(u*e-v*d)*cur->z) / n;
+			k = (-u*(q*f-r*e)*cur->x + v*(p*f-r*d)*cur->y + -w*(p*e-q*d)*cur->z) / n;
+			l = (d*(q*w-r*v)*cur->x + -e*(p*w-r*u)*cur->y + f*(p*v-q*u)*cur->z) / n;
+			
+			/* Calculate the deviations in terms of |a|, |b| and |c| */
+			h -= floor(h);  k -= floor(k);  l -= floor(l);
+			if ( h*h + k*k + l*l <= 0.1*0.1*0.1 ) n_indexed++;
+			
+			n_counted++;
+		
+		}
+		
+		cur = cur->next;
+
+	}
+	
+	return (double)n_indexed / n_counted;
 
+}
+
+static Basis *ipr_choose_initial_basis(ControlContext *ctx) {
+
+	Basis *basis;
+	Basis *best_basis;
+	double best_efom;
+	int i;
+	
+	/* Track the best basis throughout the whole procedure */
+	best_basis = NULL;  best_efom = 0;
+	
+	printf("IP: Finding initial basis using eFOM...\n");
+	for ( i=1; i<=1000; i++ ) {
+		
+		double efom;
+		
+		do { 
+			basis = ipr_choose_random_basis(ctx);
+		} while ( !basis );
+		efom = ipr_efom(ctx->reflectionctx, basis);
+		
+		if ( (efom > best_efom) || !best_basis ) {
+			if ( best_basis ) free(best_basis);
+			best_efom = efom;
+			best_basis = basis;
+			printf("IP:    %3i: eFOM=%f\n", i, efom);
+		
+		}
+		
+	}
+	
 	printf("IP: Initial basis:\n");
 	printf("IP: a = [%10.2e %10.2e %10.2e] mod=%10.2e\n", basis->a.x, basis->a.y, basis->a.z, basis->a.modulus);
 	printf("IP: b = [%10.2e %10.2e %10.2e] mod=%10.2e\n", basis->b.x, basis->b.y, basis->b.z, basis->b.modulus);
 	printf("IP: c = [%10.2e %10.2e %10.2e] mod=%10.2e\n", basis->c.x, basis->c.y, basis->c.z, basis->c.modulus);
 	
 	return basis;
-
 }
 
 /* Generate list of reflections to analyse */