Use reference reflections for scaling

1 files changed, 89 insertions, 47 deletions

diff --git a/src/hrs-scaling.c b/src/hrs-scaling.c
index bde7ad98..9b5005d8 100644
--- a/src/hrs-scaling.c
+++ b/src/hrs-scaling.c
@@ -135,8 +135,62 @@ static double s_vh(struct image *images, int n,
 }
 
 
+static gsl_vector *solve_by_eigenvalue_filtration(gsl_vector *v, gsl_matrix *M)
+{
+	gsl_eigen_symmv_workspace *work;
+	gsl_vector *e_val;
+	gsl_matrix *e_vec;
+	int val, n, frame;
+	gsl_vector *shifts;
+
+	n = v->size;
+	if ( v->size != M->size1 ) return NULL;
+	if ( v->size != M->size2 ) return NULL;
+
+	/* Diagonalise */
+	work = gsl_eigen_symmv_alloc(n);
+	e_val = gsl_vector_alloc(n);
+	e_vec = gsl_matrix_alloc(n, n);
+	val = gsl_eigen_symmv(M, e_val, e_vec, work);
+	if ( val ) {
+		ERROR("Couldn't diagonalise matrix.\n");
+		return NULL;
+	}
+	gsl_eigen_symmv_free(work);
+	val = gsl_eigen_symmv_sort(e_val, e_vec, GSL_EIGEN_SORT_ABS_DESC);
+
+	/* Rotate the "solution vector" */
+	gsl_vector *rprime;
+	rprime = gsl_vector_alloc(n);
+	val = gsl_blas_dgemv(CblasTrans, 1.0, e_vec, v, 0.0, rprime);
+
+	/* Solve (now easy) */
+	gsl_vector *sprime;
+	sprime = gsl_vector_alloc(n);
+	for ( frame=0; frame<n-1; frame++ ) {
+		double num, den;
+		num = gsl_vector_get(rprime, frame);
+		den = gsl_vector_get(e_val, frame);
+		gsl_vector_set(sprime, frame, num/den);
+	}
+	gsl_vector_set(sprime, n-1, 0.0);  /* Set last shift to zero */
+
+	/* Rotate back */
+	shifts = gsl_vector_alloc(n);
+	val = gsl_blas_dgemv(CblasNoTrans, 1.0, e_vec, sprime, 0.0, shifts);
+
+	gsl_vector_free(sprime);
+	gsl_vector_free(rprime);
+	gsl_matrix_free(e_vec);
+	gsl_vector_free(e_val);
+
+	return shifts;
+}
+
+
 static double iterate_scale(struct image *images, int n,
-                            ReflItemList *obs, const char *sym, char *cref)
+                            ReflItemList *obs, const char *sym, char *cref,
+                            double *reference)
 {
 	gsl_matrix *M;
 	gsl_vector *v;
@@ -197,16 +251,26 @@ static double iterate_scale(struct image *images, int n,
 
 			int b;  /* Frame (scale factor) number */
 			struct image *image_a = &images[a];
-			double vc, Ih, uh, vh, rha, vha, uha;
+			double vc, Ih, uh, rha, vha, uha;
 			double vval;
 
 			/* Determine the "solution" vector component */
 			uha = uha_arr[a];
 			vha = vha_arr[a];
 			uh = lookup_intensity(uh_arr, h, k, l);
-			vh = lookup_intensity(vh_arr, h, k, l);
-			Ih = vh / uh;
-			if ( isnan(Ih) ) Ih = 0.0;  /* 0 / 0 = 0, not NaN */
+
+
+			if ( !reference ) {
+				double vh;
+				vh = lookup_intensity(vh_arr, h, k, l);
+				Ih = vh / uh;
+				/* 0 / 0 = 0, not NaN */
+				if ( isnan(Ih) ) Ih = 0.0;
+			} else {
+				/* Look up by asymmetric indices */
+				Ih = lookup_intensity(reference, h, k, l);
+			}
+
 			rha = vha - image_a->osf * uha * Ih;
 			vc = Ih * rha;
 
@@ -220,6 +284,9 @@ static double iterate_scale(struct image *images, int n,
 				/* Matrix is symmetric */
 				if ( b > a ) continue;
 
+				/* Off-diagonals zero if reference available */
+				if ( (reference != NULL) && (a!=b) ) continue;
+
 				/* Zero if no common reflections */
 				if ( cref[a+n*b] != 0 ) continue;
 
@@ -252,43 +319,24 @@ static double iterate_scale(struct image *images, int n,
 	free(uha_arr);
 	free(vha_arr);
 
-	/* Fox and Holmes method */
-	gsl_eigen_symmv_workspace *work;
-	gsl_vector *e_val;
-	gsl_matrix *e_vec;
-	int val;
-
-	/* Diagonalise */
-	work = gsl_eigen_symmv_alloc(n);
-	e_val = gsl_vector_alloc(n);
-	e_vec = gsl_matrix_alloc(n, n);
-	val = gsl_eigen_symmv(M, e_val, e_vec, work);
-	if ( val ) {
-		ERROR("Couldn't diagonalise matrix.\n");
-		return 0.0;
+	show_matrix_eqn(M, v, n);
+
+	if ( reference == NULL ) {
+		shifts = solve_by_eigenvalue_filtration(v, M);
+	} else {
+		shifts = gsl_vector_alloc(n);
+		for ( frame=0; frame<n-1; frame++ ) {
+			double num, den;
+			num = gsl_vector_get(v, frame);
+			den = gsl_matrix_get(M, frame, frame);
+			gsl_vector_set(shifts, frame, num/den);
+		}
 	}
-	gsl_eigen_symmv_free(work);
-	val = gsl_eigen_symmv_sort(e_val, e_vec, GSL_EIGEN_SORT_ABS_DESC);
 
-	/* Rotate the "solution vector" */
-	gsl_vector *rprime;
-	rprime = gsl_vector_alloc(n);
-	val = gsl_blas_dgemv(CblasTrans, 1.0, e_vec, v, 0.0, rprime);
-
-	/* Solve (now easy) */
-	gsl_vector *sprime;
-	sprime = gsl_vector_alloc(n);
-	for ( frame=0; frame<n-1; frame++ ) {
-		double num, den;
-		num = gsl_vector_get(rprime, frame);
-		den = gsl_vector_get(e_val, frame);
-		gsl_vector_set(sprime, frame, num/den);
-	}
-	gsl_vector_set(sprime, n-1, 0.0);  /* Set last shift to zero */
+	gsl_matrix_free(M);
+	gsl_vector_free(v);
 
-	/* Rotate back */
-	shifts = gsl_vector_alloc(n);
-	val = gsl_blas_dgemv(CblasNoTrans, 1.0, e_vec, sprime, 0.0, shifts);
+	if ( shifts == NULL ) return 0.0;
 
 	/* Apply shifts */
 	max_shift = 0.0;
@@ -305,12 +353,6 @@ static double iterate_scale(struct image *images, int n,
 	}
 
 	gsl_vector_free(shifts);
-	gsl_vector_free(sprime);
-	gsl_vector_free(rprime);
-	gsl_matrix_free(e_vec);
-	gsl_vector_free(e_val);
-	gsl_matrix_free(M);
-	gsl_vector_free(v);
 
 	return max_shift;
 }
@@ -394,7 +436,7 @@ static void normalise_osfs(struct image *images, int n)
 
 /* Scale the stack of images */
 RefList *scale_intensities(struct image *images, int n, const char *sym,
-                           ReflItemList *obs, char *cref)
+                           ReflItemList *obs, char *cref, double *reference)
 {
 	int m;
 	RefList *full;
@@ -408,7 +450,7 @@ RefList *scale_intensities(struct image *images, int n, const char *sym,
 	i = 0;
 	do {
 
-		max_shift = iterate_scale(images, n, obs, sym, cref);
+		max_shift = iterate_scale(images, n, obs, sym, cref, reference);
 		STATUS("Scaling iteration %2i: max shift = %5.2f\n",
 		       i, max_shift);
 		i++;