First round of scaling optimisations

2 files changed, 39 insertions, 11 deletions

diff --git a/src/hrs-scaling.c b/src/hrs-scaling.c
index 0b0b700a..46cb9cc6 100644
--- a/src/hrs-scaling.c
+++ b/src/hrs-scaling.c
@@ -91,15 +91,12 @@ static double s_vha(signed int hat, signed int kat, signed int lat,
 		for ( hi=0; hi<image->n_cpeaks; hi++ ) {
 
 			double ic, sigi;
-			signed int ha, ka, la;
 
 			if ( !spots[hi].scalable ) continue;
 
-			get_asymm(spots[hi].h, spots[hi].k, spots[hi].l,
-				  &ha, &ka, &la, sym);
-			if ( ha != hat ) continue;
-			if ( ka != kat ) continue;
-			if ( la != lat ) continue;
+			if ( spots[hi].h != hat ) continue;
+			if ( spots[hi].k != kat ) continue;
+			if ( spots[hi].l != lat ) continue;
 
 			ic = spots[hi].intensity / spots[hi].p;
 			sigi = sqrt(fabs(ic));
@@ -153,15 +150,32 @@ static double iterate_scale(struct image *images, int n,
 	int a;
 	double max_shift;
 	int n_ref;
+	double *uh_arr;
+	double *vh_arr;
+	int h;  /* Reflection index */
 
 	M = gsl_matrix_calloc(n, n);
 	v = gsl_vector_calloc(n);
 	n_ref = num_items(obs);
 
+	uh_arr = new_list_intensity();
+	vh_arr = new_list_intensity();
+	for ( h=0; h<n_ref; h++ ) {
+
+		double uh, vh;
+		struct refl_item *it = get_item(obs, h);
+
+		uh = s_uh(images, n, it->h, it->k, it->l, sym);
+		vh = s_vh(images, n, it->h, it->k, it->l, sym);
+
+		set_intensity(uh_arr, it->h, it->k, it->l, uh);
+		set_intensity(vh_arr, it->h, it->k, it->l, vh);
+
+	}
+
 	for ( a=0; a<n; a++ ) {  /* "Equation number": one equation per frame */
 
 		int b;  /* Frame (scale factor) number */
-		int h;  /* Reflection index */
 		double vc_tot = 0.0;
 		struct image *image_a = &images[a];
 
@@ -176,8 +190,8 @@ static double iterate_scale(struct image *images, int n,
 			/* Determine the "solution" vector component */
 			vha = s_vha(h, k, l, images, n, sym, a);
 			uha = s_uha(h, k, l, images, n, sym, a);
-			uh = s_uh(images, n, h, k, l, sym);
-			vh = s_vh(images, n, h, k, l, sym);
+			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 */
 			rha = vha - image_a->osf * uha * Ih;
diff --git a/src/partialator.c b/src/partialator.c
index 420638bf..d1ad3768 100644
--- a/src/partialator.c
+++ b/src/partialator.c
@@ -163,7 +163,19 @@ static void refine_all(struct image *images, int n_total_patterns,
 }
 
 
-static void integrate_image(struct image *image, ReflItemList *obs)
+static void uniquify(struct cpeak *spot, const char *sym)
+{
+	signed int ha, ka, la;
+
+	get_asymm(spot->h, spot->k, spot->l, &ha, &ka, &la, sym);
+	spot->h = ha;
+	spot->k = ka;
+	spot->l = la;
+}
+
+
+static void integrate_image(struct image *image, ReflItemList *obs,
+                            const char *sym)
 {
 	struct cpeak *spots;
 	int j, n;
@@ -196,6 +208,8 @@ static void integrate_image(struct image *image, ReflItemList *obs)
 		float i_partial;
 		float xc, yc;
 
+		uniquify(&spots[j], sym);
+
 		h = spots[j].h;
 		k = spots[j].k;
 		l = spots[j].l;
@@ -449,7 +463,7 @@ int main(int argc, char *argv[])
 
 		/* Get reflections from this image.
 		 * FIXME: Use the ones from the stream */
-		integrate_image(&images[i], obs);
+		integrate_image(&images[i], obs, sym);
 
 		progress_bar(i, n_total_patterns-1, "Loading pattern data");