Initial simplex refinement stuff

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

1 files changed, 86 insertions, 331 deletions

diff --git a/src/refine.c b/src/refine.c
index b18e5fb..45aec8e 100644
--- a/src/refine.c
+++ b/src/refine.c
@@ -19,394 +19,149 @@
 #include <stdio.h>
 #include <assert.h>
 #include <string.h>
-#include <gsl/gsl_matrix.h>
-#include <gsl/gsl_vector.h>
-#include <gsl/gsl_linalg.h>
 
+#include "control.h"
 #include "displaywindow.h"
-#include "gtk-valuegraph.h"
-#include "basis.h"
-#include "reflections.h"
 #include "image.h"
 #include "reproject.h"
-#include "control.h"
 #include "mapping.h"
-#include "imagedisplay.h"
-#include "utils.h"
 
-typedef enum {
-	AXIS_X = 1,
-	AXIS_Y = 2
-} Axis;
+/* A simplex is an array of ten of these */
+typedef struct {
+	double dax;	double dbx;	double dcx;
+	double day;	double dby;	double dcy;
+	double daz;	double dbz;	double dcz;
+} SimplexVertex;
+
+typedef struct {
+	signed int h;	signed int k;	signed int l;
+	double dx;	double dy;	double dz;
+} Deviation;
+
+void refine_do_sequence(ControlContext *ctx) {
+
+	
+}
 
-static void refine_fix_unconstrained(gsl_matrix *M) {
+void refine_do_cell(ControlContext *ctx) {
 
-	int row, j;
-	int modified = 0;
+	SimplexVertex s[10];
+	Deviation *d;
+	const double delta = 0.1e9;
+	int i, nf, f, v_worst;
+	double fom_worst;
 	
-	/* Find a row which is all zero */
-	for ( row=0; row<M->size1; row++ ) {
+	if ( !ctx->cell_lattice ) {
+		displaywindow_error("No reciprocal unit cell has been found.", ctx->dw);
+		return;
+	}
 	
-		int zero = 1;
+	if ( ctx->images->n_images == 0 ) {
+		displaywindow_error("There are no images to refine against.", ctx->dw);
+		return;
+	}
+	
+	/* Initialise the simplex */
+	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;
+	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[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[9], &s[0], sizeof(SimplexVertex));	s[9].dcz = delta;
+	
+	/* Create the table of indicies and deviations */
+	nf = 0;
+	for ( i=0; i<ctx->images->n_images; i++ ) {
+		int j;
 		
-		for ( j=0; j<M->size2; j++ ) {
-			if ( gsl_matrix_get(M, row, j) != 0.0 ) {
-				zero = 0;
-				break;
-			}
+		if ( !ctx->images->images[i].rflist ) {
+			ctx->images->images[i].rflist = reproject_get_reflections(&ctx->images->images[i], ctx->cell_lattice);
 		}
 		
-		if ( zero ) {
-			
-			/* Find a column which is all zero */
-			int i, col;
-			int row_altered = 0;
-			
-			for ( col=0; col<M->size2; col++ ) {
-			
-				int zero2 = 1;
-				
-				if ( row_altered ) break; /* This row is now fixed, so move on to the next */
-				
-				for ( i=0; i<M->size1; i++ ) {
-					if ( gsl_matrix_get(M, i, col) != 0.0 ) {
-						zero2 = 0;
-						break;
-					}
-				}
-				
-				if ( zero2 ) {
-					gsl_matrix_set(M, row, col, 1.0);
-					modified = 1;
-					row_altered = 1;
-				}
-			
-			}
-			
+		for ( j=0; j<ctx->images->images[i].rflist->n_features; j++ ) {
+			if ( ctx->images->images[i].rflist->features[j].partner != NULL ) nf++;
 		}
 		
 	}
+	printf("RF: There are %i partnered features in total\n", nf);
 	
-	if ( modified) printf("RF: M was modified to fix unconstrained variables.\n");
-
-}
-
-/* Use the IPR algorithm to make "cell" fit the given image */
-ImageFeature *refine_fit_image(Basis *cell, ImageRecord *image) {
-	
-	ImageFeatureList *flist;
-	int j, ns;
-	Axis axis;
-	double dax = 0.0, dbx = 0.0, dcx = 0.0;
-	double day = 0.0, dby = 0.0, dcy = 0.0;
-	double dlx = 0.0, dly = 0.0, dlz = 0.0;
-	
-	flist = image->features;
-		
-	for ( axis = AXIS_X; axis <= AXIS_Y; axis++ ) {
-	
-		gsl_permutation *perm;
-		int s;
-		double det;
-		gsl_matrix *M;
-		gsl_vector *q;
-		gsl_vector *p;
-		gsl_vector *scratch;
-		gsl_matrix *V;
-		gsl_vector *S;
+	d = malloc(nf*sizeof(Deviation));
+	f = 0;
+	for ( i=0; i<ctx->images->n_images; i++ ) {
 		
-		if ( axis == AXIS_X ) printf("RF: ------------------------------------------------------------------ Refining x coordinates -----\n");
-		if ( axis == AXIS_Y ) printf("RF: ------------------------------------------------------------------ Refining y coordinates -----\n");
+		ImageRecord *image;
+		int j;
 		
-		M = gsl_matrix_alloc(3, 3);
-		p = gsl_vector_alloc(3);
-		gsl_matrix_set_zero(M);
-		gsl_vector_set_zero(p);
+		image = &ctx->images->images[i];
 		
-		ns = 0;
-		for ( j=0; j<image->rflist->n_features; j++ ) {
+		for ( j=0; j<ctx->images->images[i].features->n_features; j++ ) {
 		
-			double val;
 			ImageFeature *rf;
-			signed int h, k, l;
-			double xy;
 			double dix, diy, dx, dy;
+			double dlx, dly, dlz;
 			double old_x, old_y;
 			
 			rf = &image->rflist->features[j];
-			
 			if ( !rf->partner ) continue;
 			
-			h = rf->reflection->h;
-			k = rf->reflection->k;
-			l = rf->reflection->l;
+			d[f].h = rf->reflection->h;
+			d[f].k = rf->reflection->k;
+			d[f].l = rf->reflection->l;
 			
 			/* Determine the difference vector */
 			dix = rf->partner->x - rf->x;
 			diy = rf->partner->y - rf->y;
-			printf("RF: Feature %3i: %3i %3i %3i dev = %+9.5f %+9.5f px ", j, h, k, l, dix, diy);
+			printf("RF: Feature %3i: %3i %3i %3i dev = %+9.5f %+9.5f px ", j, d[f].h, d[f].k, d[f].l, dix, diy);
 			
 			old_x = rf->partner->x;
 			old_y = rf->partner->y;
 			rf->partner->x = dix + rf->partner->parent->x_centre;
 			rf->partner->y = diy + rf->partner->parent->y_centre;
 			mapping_scale(rf->partner, &dx, &dy);
+			mapping_rotate(dx, dy, 0.0, &dlx, &dly, &dlz, image->omega, image->tilt);
 			rf->partner->x = old_x;
 			rf->partner->y = old_y;
-			printf("=> %+10.5f %+10.5f nm^-1\n", dx/1e9, dy/1e9);
-			
-			xy = 0;
-			switch ( axis ) {
-				case AXIS_X : xy = dx; break;
-				case AXIS_Y : xy = dy; break;
-			}
+			printf("=> %+10.5f %+10.5f %+10.5f nm^-1\n", dlx/1e9, dly/1e9, dlz/1e9);
 			
-			/* Elements of "p" */
-			val = gsl_vector_get(p, 0);	val += xy * h;	gsl_vector_set(p, 0, val);
-			val = gsl_vector_get(p, 1);	val += xy * k;	gsl_vector_set(p, 1, val);
-			val = gsl_vector_get(p, 2);	val += xy * l;	gsl_vector_set(p, 2, val);
-			gsl_matrix_get(M, 2, 2);
-			/* Elements of "M" */
-			val = gsl_matrix_get(M, 0, 0);	val += h * h;	gsl_matrix_set(M, 0, 0, val);
-			val = gsl_matrix_get(M, 0, 1);	val += k * h;	gsl_matrix_set(M, 0, 1, val);
-			val = gsl_matrix_get(M, 0, 2);	val += l * h;	gsl_matrix_set(M, 0, 2, val);
+			d[f].dx = dlx;
+			d[f].dy = dly;
+			d[f].dz = dlz;
 			
-			val = gsl_matrix_get(M, 1, 0);	val += h * k;	gsl_matrix_set(M, 1, 0, val);
-			val = gsl_matrix_get(M, 1, 1);	val += k * k;	gsl_matrix_set(M, 1, 1, val);
-			val = gsl_matrix_get(M, 1, 2);	val += l * k;	gsl_matrix_set(M, 1, 2, val);
-			
-			val = gsl_matrix_get(M, 2, 0);	val += h * l;	gsl_matrix_set(M, 2, 0, val);
-			val = gsl_matrix_get(M, 2, 1);	val += k * l;	gsl_matrix_set(M, 2, 1, val);
-			val = gsl_matrix_get(M, 2, 2);	val += l * l;	gsl_matrix_set(M, 2, 2, val);
-			
-			ns++;
+			f++;
 		
 		}
 		
-		if ( ns == 0 ) {
-			printf("RF: No partners found\n");
-			gsl_matrix_free(M);
-			gsl_vector_free(p);
-			return NULL;
-		}	
-		
-		/* Prepare for solving */
-		refine_fix_unconstrained(M);
-	//	gsl_matrix_set(M, 2, 0, 0.0);
-	//	gsl_matrix_set(M, 2, 1, 0.0);
-	//	gsl_matrix_set(M, 2, 2, 1.0);
-	//	gsl_matrix_set(M, 0, 2, 0.0);
-	//	gsl_matrix_set(M, 1, 2, 0.0);
-	//	gsl_vector_set(p, 2, 0.0);
-		matrix_vector_show(M, p, "RF: ");
-		
-		/* Calculate and display the determinant */
-		perm = gsl_permutation_alloc(M->size1);
-		gsl_linalg_LU_decomp(M, perm, &s);
-		det = gsl_linalg_LU_det(M, s);
-		printf("RF: Determinant of M = %f\n", det);
-		gsl_permutation_free(perm);		
-		
-		/* Solve the equations */
-		V = gsl_matrix_alloc(M->size2, M->size2);
-		S = gsl_vector_alloc(M->size2);
-		scratch = gsl_vector_alloc(M->size2);
-		gsl_linalg_SV_decomp(M, V, S, scratch);
-		q = gsl_vector_alloc(3);	/* This is the "answer" */
-		gsl_vector_set_zero(q);
-		gsl_linalg_SV_solve(M, V, S, p, q);
-		gsl_vector_free(scratch);
-		gsl_matrix_free(V);
-		gsl_vector_free(S);
-		gsl_matrix_free(M);
-		gsl_vector_free(p);
-		
-		switch ( axis ) {
-			case AXIS_X : {
-				dax = gsl_vector_get(q, 0);
-				dbx = gsl_vector_get(q, 1);	/* These are the deviations, in the direction "x" of the image coordinate  */
-				dcx = gsl_vector_get(q, 2);	/* system, of a,b and c */
-				break;
-			}
-			case AXIS_Y : {
-				day = gsl_vector_get(q, 0);
-				dby = gsl_vector_get(q, 1);	/* These are the deviations, in the direction "y" of the image coordinate  */
-				dcy = gsl_vector_get(q, 2);	/* system, of a,b and c */
-				break;
-			}
-		}
-		
-		gsl_vector_free(q);
-
 	}
+	assert( f == nf );
 	
-	printf("RF: ------------------------------------------------------------------ Refinement results ---------\n");
-	printf("RF: a should change by %+7.5f %+7.5f nm^-1 in the image plane\n", dax/1e9, day/1e9);
-	printf("RF: b should change by %+7.5f %+7.5f nm^-1 in the image plane\n", dbx/1e9, dby/1e9);
-	printf("RF: c should change by %+7.5f %+7.5f nm^-1 in the image plane\n", dcx/1e9, dcy/1e9);
-	
-	/* Update the cell */
-	mapping_rotate(dax, day, 0.0, &dlx, &dly, &dlz, image->omega, image->tilt);
-	printf("RF: a changed by [ %+7.5f %+7.5f %+7.5f ] nm^-1 in reciprocal space\n", dlx/1e9, dly/1e9, dlz/1e9);
-	cell->a.x += dlx;  cell->a.y += dly; cell->a.z += dlz;
-	
-	mapping_rotate(dbx, dby, 0.0, &dlx, &dly, &dlz, image->omega, image->tilt);
-	printf("RF: b changed by [ %+7.5f %+7.5f %+7.5f ] nm^-1 in reciprocal space\n", dlx/1e9, dly/1e9, dlz/1e9);
-	cell->b.x += dlx;  cell->b.y += dly;  cell->b.z += dlz;
+	/* Find the least favourable vertex of the simplex */
+	v_worst = 0;
+	fom_worst = 0;
+	for ( i=0; i<10; i++ ) {
 	
-	mapping_rotate(dcx, dcy, 0.0, &dlx, &dly, &dlz, image->omega, image->tilt);
-	printf("RF: c changed by [ %+7.5f %+7.5f %+7.5f ] nm^-1 in reciprocal space\n", dlx/1e9, dly/1e9, dlz/1e9);
-	cell->c.x += dlx;  cell->c.y += dly;  cell->c.z += dlz;
-	
-	return NULL;
-	
-}
-
-static int refine_sequence_sweep(ControlContext *ctx, double *fit, double *warp) {
-
-	int i;
-	double series_dev_max = 0;
-	double series_dev_min = +HUGE_VAL;
-	double series_dev_mean = 0;
-	int series_dev_n = 0;
-	
-	/* Ensure that ctx->cell_lattice is set up */
-	reproject_cell_to_lattice(ctx);
-	
-	for ( i=0; i<ctx->images->n_images; i++ ) {
-
-		ImageRecord *image;
-		int j, n;
-		double image_dev_mean = 0;
-		
-		image = &ctx->images->images[i];
+		double fom = 0;
+		int j;
 		
-		/* Fit this image and update ctx->cell_lattice, index the selected pattern */
-		if ( !image->rflist ) image->rflist = reproject_get_reflections(image, ctx->cell_lattice);
-		refine_fit_image(ctx->cell, image);
-		reproject_cell_to_lattice(ctx);
-		image->rflist = reproject_get_reflections(image, ctx->cell_lattice);
+		for ( j=0; j<nf; j++ ) {
 		
-		n = 0;
-		for ( j=0; j<image->rflist->n_features; j++ ) {
+			double xdev, ydev, zdev;
 			
-			double dix, diy;
+			xdev = 0;
+			ydev = 0;
+			zdev = 0;
 			
-			/* Skip if no partner */
-			if ( !image->rflist->features[j].partner ) continue;
+			fom += sqrt(xdev*xdev + ydev*ydev + zdev*zdev);
 			
-			/* Determine the difference vector */
-			dix = image->rflist->features[j].partner->x - image->rflist->features[j].x;
-			diy = image->rflist->features[j].partner->y - image->rflist->features[j].y;
-			
-			image_dev_mean += sqrt(dix*dix + diy*diy);
-			n++;
-			
-		}
-		image_dev_mean /= n;
-		
-		if ( image_dev_mean > series_dev_max ) series_dev_max = image_dev_mean;
-		if ( image_dev_mean < series_dev_min ) series_dev_min = image_dev_mean;
-		series_dev_mean += image_dev_mean;
-		series_dev_n++;
-		
-	}
-	
-	series_dev_mean /= series_dev_n;
-	*fit = series_dev_mean;
-	*warp = (series_dev_max - series_dev_min)/series_dev_min;
-	
-	return 0;
-
-}
-
-void refine_do_stack(ControlContext *ctx) {
-
-	double omega_offs;
-	GtkWidget *window_fit;
-	GtkWidget *graph_fit;
-	double *fit_vals;
-	GtkWidget *window_warp;
-	GtkWidget *graph_warp;
-	double *warp_vals;
-	size_t idx;
-	
-	fit_vals = malloc(401*sizeof(double));
-	warp_vals = malloc(401*sizeof(double));
-	idx = 0;
-
-	if ( !ctx->cell_lattice ) {
-		displaywindow_error("No reciprocal unit cell has been found.", ctx->dw);
-		return;
-	}
-	
-	if ( ctx->images->n_images == 0 ) {
-		displaywindow_error("There are no images to refine against.", ctx->dw);
-		return;
-	}
-	
-	for ( omega_offs=-2.0; omega_offs<=2.0; omega_offs+=0.01 ) {
-	
-		double fit, warp;
-		int i;
-		Basis cell_copy;
-		
-		memcpy(&cell_copy, ctx->cell, sizeof(Basis));
-		for ( i=0; i<ctx->images->n_images; i++ ) {
-			ctx->images->images[i].omega += omega_offs;
-		}
-		
-		if ( refine_sequence_sweep(ctx, &fit, &warp) ) {
-			printf("RF: Sequencer sweep failed\n");
-		}
-		printf("RF: omega_offs=%f, fit=%f, warp=%f\n", omega_offs, fit, warp);
-		fit_vals[idx] = fit;
-		warp_vals[idx++] = warp;
-		
-		for ( i=0; i<ctx->images->n_images; i++ ) {
-			ctx->images->images[i].omega -= omega_offs;
 		}
-		memcpy(ctx->cell, &cell_copy, sizeof(Basis));
 		
 	}
 	
-	displaywindow_update(ctx->dw);
-	reproject_lattice_changed(ctx);
-	
-	window_fit = gtk_window_new(GTK_WINDOW_TOPLEVEL);
-	gtk_window_set_default_size(GTK_WINDOW(window_fit), 640, 256);
-	gtk_window_set_title(GTK_WINDOW(window_fit), "Omega-Search Graph: Fit");
-	graph_fit = gtk_value_graph_new();
-	gtk_value_graph_set_data(GTK_VALUE_GRAPH(graph_fit), fit_vals, idx);
-	gtk_container_add(GTK_CONTAINER(window_fit), graph_fit);
-	gtk_widget_show_all(window_fit);
-	
-	window_warp = gtk_window_new(GTK_WINDOW_TOPLEVEL);
-	gtk_window_set_default_size(GTK_WINDOW(window_warp), 640, 256);
-	gtk_window_set_title(GTK_WINDOW(window_warp), "Omega-Search Graph: Warp");
-	graph_warp = gtk_value_graph_new();
-	gtk_value_graph_set_data(GTK_VALUE_GRAPH(graph_warp), warp_vals, idx);
-	gtk_container_add(GTK_CONTAINER(window_warp), graph_warp);
-	gtk_widget_show_all(window_warp);
-	
-}
-
-void refine_do_image(ControlContext *ctx) {
-
-	if ( !ctx->cell_lattice ) {
-		displaywindow_error("No reciprocal unit cell has been found.", ctx->dw);
-		return;
-	}
-	
-	if ( ctx->images->n_images == 0 ) {
-		displaywindow_error("There are no images to refine against.", ctx->dw);
-		return;
-	}
-		
-	ImageFeature *fitted;
-	
-	fitted = refine_fit_image(ctx->cell, &ctx->images->images[ctx->dw->cur_image]);
-	
 	ctx->images->images[ctx->dw->cur_image].rflist = NULL;
 	reproject_lattice_changed(ctx);
 	displaywindow_update(ctx->dw);