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

1 files changed, 63 insertions, 2 deletions

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);