Automatically fix unconstrained variables in refinement

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

2 files changed, 65 insertions, 4 deletions

diff --git a/src/iprtest.c b/src/iprtest.c
index 321c095..cb2fb10 100644
--- a/src/iprtest.c
+++ b/src/iprtest.c
@@ -117,7 +117,7 @@ int main(int argc, char *argv[]) {
 	
 	ctx->image = malloc(sizeof(ImageRecord));
 	ctx->image->image = NULL;
-	ctx->image->tilt = 0.0;
+	ctx->image->tilt = 20.0;
 	ctx->image->omega = 0.0;
 	ctx->image->slop = 0.0;
 	ctx->image->fmode = FORMULATION_PIXELSIZE;
@@ -141,7 +141,7 @@ int main(int argc, char *argv[]) {
 	ctx->image->features = reproject_get_reflections(ctx->image, ctx->reflections);
 	/* The "model" cell */
 	ctx->cell->a.x = 5.2e9;  ctx->cell->a.y = 0.1e9;  ctx->cell->a.z = 0.0;
-	ctx->cell->b.x = 0.1e9;  ctx->cell->b.y = 4.1e9;  ctx->cell->b.z = 0.0;
+	ctx->cell->b.x = 0.1e9;  ctx->cell->b.y = 4.5e9;  ctx->cell->b.z = 0.0;
 	ctx->cell->c.x = 0.0;  ctx->cell->c.y = 0.0;  ctx->cell->c.z = 5.0e9;
 	ctx->gen = reflection_list_from_cell(ctx->cell);
 	ctx->image->rflist = reproject_get_reflections(ctx->image, ctx->gen);
diff --git a/src/refine.c b/src/refine.c
index 7d77ef6..14681c9 100644
--- a/src/refine.c
+++ b/src/refine.c
@@ -39,6 +39,59 @@ typedef enum {
 	AXIS_Y = 2
 } Axis;
 
+static void refine_fix_unconstrained(gsl_matrix *M) {
+
+	int row, j;
+	int modified = 0;
+	
+	/* Find a row which is all zero */
+	for ( row=0; row<M->size1; row++ ) {
+	
+		int zero = 1;
+		
+		for ( j=0; j<M->size2; j++ ) {
+			if ( gsl_matrix_get(M, row, j) != 0.0 ) {
+				zero = 0;
+				break;
+			}
+		}
+		
+		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;
+				
+				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;
+					continue;
+				}
+			
+			}
+			
+		}
+		
+	}
+	
+	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) {
 	
@@ -89,6 +142,15 @@ ImageFeature *refine_fit_image(Basis *cell, ImageRecord *image) {
 			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);
 			
+			//if ( (h<0) || (k<0) || (l<0) ) {
+			//	h = -h;
+			//	k = -k;
+			//	l = -l;
+			//	dix = -dix;
+			//	diy = -diy;
+			//	printf("\n          Fudge: %3i %3i %3i dev = %+9.5f %+9.5f px ", h, k, l, dix, diy);
+			//}
+			
 			old_x = rf->partner->x;
 			old_y = rf->partner->y;
 			rf->partner->x = dix + rf->partner->parent->x_centre;
@@ -134,8 +196,7 @@ ImageFeature *refine_fit_image(Basis *cell, ImageRecord *image) {
 		}	
 		
 		/* Do the fitting */
-		//gsl_matrix_set(M, 1, 1, 1.0);
-		gsl_matrix_set(M, 2, 2, 1.0);
+		refine_fix_unconstrained(M);
 		matrix_vector_show(M, p);
 		perm = gsl_permutation_alloc(M->size1);
 		gsl_linalg_LU_decomp(M, perm, &s);