Use reference reflections for scaling

3 files changed, 104 insertions, 54 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++;
diff --git a/src/hrs-scaling.h b/src/hrs-scaling.h
index d65129b8..f65b7988 100644
--- a/src/hrs-scaling.h
+++ b/src/hrs-scaling.h
@@ -21,7 +21,8 @@
 #include "image.h"
 
 extern RefList *scale_intensities(struct image *images, int n, const char *sym,
-                                  ReflItemList *obs, char *cref);
+                                  ReflItemList *obs, char *cref,
+                                  double *reference);
 
 extern char *find_common_reflections(struct image *images, int n);
 
diff --git a/src/partialator.c b/src/partialator.c
index 64546a7b..e3eccafb 100644
--- a/src/partialator.c
+++ b/src/partialator.c
@@ -171,7 +171,7 @@ int main(int argc, char *argv[])
 	char *cref;
 	int n_usable_patterns = 0;
 	char *reference_file = NULL;
-	RefList *reference;
+	double *reference = NULL;
 
 	/* Long options */
 	const struct option longopts[] = {
@@ -230,7 +230,7 @@ int main(int argc, char *argv[])
 			break;
 
 		case 1 :
-			reference = strdup(optarg);
+			reference_file = strdup(optarg);
 			break;
 
 		case 0 :
@@ -278,9 +278,15 @@ int main(int argc, char *argv[])
 	}
 
 	if ( reference_file != NULL ) {
-		reference = read_reflections(reference_file);
+		RefList *list;
+		RefList *symmed;
+		list = read_reflections(reference_file);
 		free(reference_file);
-		if ( reference == NULL ) return 1;
+		if ( list == NULL ) return 1;
+		symmed = asymmetric_indices(list, sym);
+		reflist_free(list);
+		reference = intensities_from_list(symmed);
+		reflist_free(symmed);
 	}
 
 	n_total_patterns = count_patterns(fh);
@@ -360,7 +366,7 @@ int main(int argc, char *argv[])
 	/* Make initial estimates */
 	STATUS("Performing initial scaling.\n");
 	full = scale_intensities(images, n_usable_patterns, sym,
-	                         scalable, cref);
+	                         scalable, cref, reference);
 
 	for ( i=0; i<num_items(scalable); i++ ) {
 		Reflection *f;
@@ -425,7 +431,7 @@ int main(int argc, char *argv[])
 		/* Re-estimate all the full intensities */
 		reflist_free(full);
 		full = scale_intensities(images, n_usable_patterns,
-		                         sym, scalable, cref);
+		                         sym, scalable, cref, reference);
 
 		fclose(fhg);
 		fclose(fhp);
@@ -451,6 +457,7 @@ int main(int argc, char *argv[])
 	free_detector_geometry(det);
 	free(beam);
 	free(cref);
+	free(reference);
 	for ( i=0; i<n_usable_patterns; i++ ) {
 		cell_free(images[i].indexed_cell);
 		free(images[i].filename);