Set matrix elements to zero if no common reflections (NB broken)

3 files changed, 54 insertions, 7 deletions

diff --git a/src/hrs-scaling.c b/src/hrs-scaling.c
index 99490251..e2adaaf6 100644
--- a/src/hrs-scaling.c
+++ b/src/hrs-scaling.c
@@ -35,6 +35,45 @@
 #define MAX_CYCLES (30)
 
 
+char *find_common_reflections(struct image *images, int n)
+{
+	int i;
+	char *cref;
+
+	cref = calloc(n*n, sizeof(char));
+
+	for ( i=0; i<n; i++ ) {
+
+		Reflection *refl;
+		RefListIterator *iter;
+
+		for ( refl = first_refl(images[i].reflections, &iter);
+		      refl != NULL;
+		      refl = next_refl(refl, iter) ) {
+
+			signed int h, k, l;
+			int j;
+
+			get_indices(refl, &h, &k, &l);
+
+			for ( j=0; j<i; j++ ) {
+				Reflection *r2;
+				r2 = find_refl(images[j].reflections, h, k, l);
+				if ( r2 != NULL ) {
+					cref[i+n*j] = 1;
+					cref[j+n*i] = 1;
+					break;
+				}
+			}
+
+		}
+
+	}
+
+	return cref;
+}
+
+
 static void s_uhavha(signed int hat, signed int kat, signed int lat,
                      struct image *image, double *uha, double *vha)
 {
@@ -97,7 +136,7 @@ static double s_vh(struct image *images, int n,
 
 
 static double iterate_scale(struct image *images, int n,
-                            ReflItemList *obs, const char *sym)
+                            ReflItemList *obs, const char *sym, char *cref)
 {
 	gsl_matrix *M;
 	gsl_vector *v;
@@ -181,6 +220,9 @@ static double iterate_scale(struct image *images, int n,
 				/* Matrix is symmetric */
 				if ( b > a ) continue;
 
+				/* Zero if no common reflections */
+				if ( cref[a+n*b] != 0 ) continue;
+
 				uhb = uha_arr[b];
 				vhb = vha_arr[b];
 				rhb = vhb - image_b->osf * uhb * Ih;
@@ -336,7 +378,7 @@ static void normalise_osfs(struct image *images, int n)
 
 /* Scale the stack of images */
 double *scale_intensities(struct image *images, int n, const char *sym,
-                          ReflItemList *obs)
+                          ReflItemList *obs, char *cref)
 {
 	int m;
 	double *I_full;
@@ -350,7 +392,7 @@ double *scale_intensities(struct image *images, int n, const char *sym,
 	i = 0;
 	do {
 
-		max_shift = iterate_scale(images, n, obs, sym);
+		max_shift = iterate_scale(images, n, obs, sym, cref);
 		STATUS("Iteration %2i: max shift = %5.2f\n", i, max_shift);
 		i++;
 		normalise_osfs(images, n);
diff --git a/src/hrs-scaling.h b/src/hrs-scaling.h
index ae384aa8..04a321e3 100644
--- a/src/hrs-scaling.h
+++ b/src/hrs-scaling.h
@@ -20,8 +20,9 @@
 
 #include "image.h"
 
-extern double *scale_intensities(struct image *image, int n, const char *sym,
-                                 ReflItemList *obs);
+extern double *scale_intensities(struct image *images, int n, const char *sym,
+                                 ReflItemList *obs, char *cref);
 
+extern char *find_common_reflections(struct image *images, int n);
 
 #endif	/* HRS_SCALING_H */
diff --git a/src/partialator.c b/src/partialator.c
index 309343cb..e99b2291 100644
--- a/src/partialator.c
+++ b/src/partialator.c
@@ -161,6 +161,7 @@ int main(int argc, char *argv[])
 	int n_found = 0;
 	int n_expected = 0;
 	int n_notfound = 0;
+	char *cref;
 
 	/* Long options */
 	const struct option longopts[] = {
@@ -379,12 +380,13 @@ int main(int argc, char *argv[])
 	STATUS("Mean measurements per unique reflection: %5.2f\n",
 	       (double)n_found / num_items(obs));
 
+	cref = find_common_reflections(images, n_total_patterns);
 	cts = new_list_count();
 
 	/* Make initial estimates */
 	STATUS("Performing initial scaling.\n");
 	select_scalable_reflections(images, n_total_patterns);
-	I_full = scale_intensities(images, n_total_patterns, sym, obs);
+	I_full = scale_intensities(images, n_total_patterns, sym, obs, cref);
 
 	/* Iterate */
 	for ( i=0; i<n_iter; i++ ) {
@@ -416,7 +418,8 @@ int main(int argc, char *argv[])
 		/* Re-estimate all the full intensities */
 		free(I_full);
 		select_scalable_reflections(images, n_total_patterns);
-		I_full = scale_intensities(images, n_total_patterns, sym, obs);
+		I_full = scale_intensities(images, n_total_patterns,
+		                           sym, obs, cref);
 
 		fclose(fhg);
 		fclose(fhp);
@@ -443,6 +446,7 @@ int main(int argc, char *argv[])
 	free(det);
 	free(beam);
 	free(cts);
+	free(cref);
 	for ( i=0; i<n_total_patterns; i++ ) {
 		cell_free(images[i].indexed_cell);
 		free(images[i].filename);