Work on post refinement

7 files changed, 156 insertions, 85 deletions

diff --git a/src/geometry.c b/src/geometry.c
index 749c194e..3ff39f31 100644
--- a/src/geometry.c
+++ b/src/geometry.c
@@ -23,9 +23,7 @@
 #include "peaks.h"
 #include "beam-parameters.h"
 #include "reflist.h"
-
-
-#define MAX_CPEAKS (256 * 256)
+#include "reflist-utils.h"
 
 
 static signed int locate_peak(double x, double y, double z, double k,
@@ -292,3 +290,72 @@ double integrate_all(struct image *image, RefList *reflections)
 
 	return itot;
 }
+
+
+/* Calculate partialities and apply them to the image's raw_reflections */
+void update_partialities(struct image *image, const char *sym,
+                         int *n_expected, int *n_found, int *n_notfound)
+{
+	Reflection *refl;
+	RefListIterator *iter;
+	RefList *predicted;
+
+	predicted = find_intersections(image, image->indexed_cell, 0);
+
+	for ( refl = first_refl(predicted, &iter);
+	      refl != NULL;
+	      refl = next_refl(refl, iter) )
+	{
+
+		Reflection *peak_in_pattern;
+		double r1, r2, p, x, y;
+		signed int h, k, l;
+		int clamp1, clamp2, scalable;
+
+		/* Get predicted indices and location */
+		get_indices(refl, &h, &k, &l);
+		get_detector_pos(refl, &x, &y);
+		if ( n_expected != NULL ) (*n_expected)++;
+
+		/* Look for this reflection in the pattern */
+		peak_in_pattern = find_refl(image->raw_reflections, h, k, l);
+		if ( peak_in_pattern == NULL ) {
+			if ( n_notfound != NULL ) (*n_notfound)++;
+			continue;
+		}
+		if ( n_found != NULL ) (*n_found)++;
+
+		/* Transfer partiality stuff */
+		get_partial(refl, &r1, &r2, &p, &clamp1, &clamp2);
+		set_partial(peak_in_pattern, r1, r2, p, clamp1, clamp2);
+
+		/* Transfer detector location */
+		get_detector_pos(refl, &x, &y);
+		set_detector_pos(peak_in_pattern, 0.0, x, y);
+
+	}
+
+	reflist_free(predicted);
+}
+
+
+void update_partialities_and_asymm(struct image *image, const char *sym,
+                                   ReflItemList *obs,
+                                   int *n_expected, int *n_found,
+                                   int *n_notfound)
+{
+	/* Get rid of the old list, about to be replaced */
+	reflist_free(image->reflections);
+	image->reflections = NULL;
+
+	/* Fill in partialities */
+	update_partialities(image, sym, n_expected, n_found, n_notfound);
+
+	/* Rewrite the reflections with the asymmetric indices
+	 * to get the list used for scaling and post refinement */
+	image->reflections = asymmetric_indices(image->raw_reflections,
+	                                        sym, obs);
+
+	/* Need these lists to work fast */
+	optimise_reflist(image->reflections);
+}
diff --git a/src/geometry.h b/src/geometry.h
index a0761456..f8c007fa 100644
--- a/src/geometry.h
+++ b/src/geometry.h
@@ -24,5 +24,8 @@ extern RefList *find_intersections(struct image *image, UnitCell *cell,
 
 extern double integrate_all(struct image *image, RefList *reflections);
 
-
+extern void update_partialities_and_asymm(struct image *image, const char *sym,
+                                          ReflItemList *obs,
+                                          int *n_expected, int *n_found,
+                                          int *n_notfound);
 #endif	/* GEOMETRY_H */
diff --git a/src/hrs-scaling.c b/src/hrs-scaling.c
index 03458670..4b9b932d 100644
--- a/src/hrs-scaling.c
+++ b/src/hrs-scaling.c
@@ -302,6 +302,7 @@ 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);
diff --git a/src/image.h b/src/image.h
index d6e22aee..22e0ddae 100644
--- a/src/image.h
+++ b/src/image.h
@@ -87,6 +87,7 @@ typedef struct _imagefeaturelist ImageFeatureList;
  *    int                     height;
  *
  *    RefList                 *reflections;
+ *    RefList                 *raw_reflections;
  *
  *    ImageFeatureList        *features;
  * };
@@ -150,6 +151,9 @@ struct image {
 	/* Integrated (or about-to-be-integrated) reflections */
 	RefList                 *reflections;
 
+	/* Raw version of "reflections", e.g. in point group 1 */
+	RefList                 *raw_reflections;
+
 	/* Detected peaks */
 	ImageFeatureList        *features;
 
diff --git a/src/partialator.c b/src/partialator.c
index 33dcad68..54022e9e 100644
--- a/src/partialator.c
+++ b/src/partialator.c
@@ -197,6 +197,7 @@ int main(int argc, char *argv[])
 	int n_expected = 0;
 	int n_notfound = 0;
 	char *cref;
+	int n_usable_patterns = 0;
 
 	/* Long options */
 	const struct option longopts[] = {
@@ -315,93 +316,47 @@ int main(int argc, char *argv[])
 	obs = new_items();
 	for ( i=0; i<n_total_patterns; i++ ) {
 
-		RefList *predicted;
-		RefList *measured;
-		Reflection *refl;
-		RefListIterator *iter;
+		images[n_usable_patterns].det = NULL;
 
-		if ( read_chunk(fh, &images[i]) != 0 ) {
+		if ( read_chunk(fh, &images[n_usable_patterns]) != 0 ) {
 			/* Should not happen, because we counted the patterns
 			 * earlier. */
 			ERROR("Failed to read chunk from the input stream.\n");
 			return 1;
 		}
 
-		/* FIXME: This leaves gaps in the array */
-		if ( images[i].indexed_cell == NULL ) continue;
-
 		/* Won't be needing this, if it exists */
-		image_feature_list_free(images[i].features);
-		images[i].features = NULL;
-
-		images[i].div = beam->divergence;
-		images[i].bw = beam->bandwidth;
-		images[i].det = det;
-		images[i].width = det->max_fs;
-		images[i].height = det->max_ss;
-		images[i].osf = 1.0;
-		images[i].profile_radius = 0.005e9;
-
-		/* Muppet proofing */
-		images[i].data = NULL;
-		images[i].flags = NULL;
-		images[i].beam = NULL;
-
-		/* Calculate initial partialities and fill in intensities from
-		 * the stream */
-		predicted = find_intersections(&images[i],
-		                               images[i].indexed_cell, 0);
-
-		/* We start again with a new reflection list, this time with
-		 * the asymmetric indices */
-		measured = images[i].reflections;
-		images[i].reflections = reflist_new();
-
-		for ( refl = first_refl(predicted, &iter);
-		      refl != NULL;
-		      refl = next_refl(refl, iter) ) {
+		image_feature_list_free(images[n_usable_patterns].features);
+		images[n_usable_patterns].features = NULL;
 
-			Reflection *peak_in_pattern;
-			Reflection *new;
-			signed int h, k, l, ha, ka, la;
-			double r1, r2, p, x, y;
-			int clamp1, clamp2;
-
-			/* Get predicted indices and location */
-			get_indices(refl, &h, &k, &l);
-			get_detector_pos(refl, &x, &y);
-			n_expected++;
-
-			/* Look for this reflection in the pattern */
-			peak_in_pattern = find_refl(measured, h, k, l);
-			if ( peak_in_pattern == NULL ) {
-				n_notfound++;
-				continue;
-			}
-			n_found++;
+		/* "n_usable_patterns" will not be incremented in this case */
+		if ( images[n_usable_patterns].indexed_cell == NULL ) continue;
 
-			/* Put it into the asymmetric cell */
-			get_asymm(h, k, l, &ha, &ka, &la, sym);
-			if ( find_item(obs, ha, ka, la) == 0 ) {
-				add_item(obs, ha, ka, la);
-			}
+		/* Fill in initial estimates of stuff */
+		images[n_usable_patterns].div = beam->divergence;
+		images[n_usable_patterns].bw = beam->bandwidth;
+		images[n_usable_patterns].det = det;
+		images[n_usable_patterns].width = det->max_fs;
+		images[n_usable_patterns].height = det->max_ss;
+		images[n_usable_patterns].osf = 1.0;
+		images[n_usable_patterns].profile_radius = 0.005e9;
 
-			/* Create new reflection and copy data across */
-			new = add_refl(images[i].reflections, ha, ka, la);
-			get_partial(refl, &r1, &r2, &p, &clamp1, &clamp2);
-			get_detector_pos(refl, &x, &y);
-			set_int(new, get_intensity(peak_in_pattern));
-			set_partial(new, r1, r2, p, clamp1, clamp2);
-			set_detector_pos(new, 0.0, x, y);
+		/* Muppet proofing */
+		images[n_usable_patterns].data = NULL;
+		images[n_usable_patterns].flags = NULL;
+		images[n_usable_patterns].beam = NULL;
 
-		}
-		reflist_free(measured);
-		reflist_free(predicted);
+		/* This is the raw list of reflections */
+		images[n_usable_patterns].raw_reflections =
+		                          images[n_usable_patterns].reflections;
+		images[n_usable_patterns].reflections = NULL;
 
-		/* Do magic on the reflection list to make things go faster */
-		optimise_reflist(images[i].reflections);
+		update_partialities_and_asymm(&images[n_usable_patterns], sym,
+		                              obs, &n_expected, &n_found,
+		                              &n_notfound);
 
 		progress_bar(i, n_total_patterns-1, "Loading pattern data");
+		n_usable_patterns++;
 
 	}
 	fclose(fh);
@@ -410,12 +365,12 @@ 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);
+	cref = find_common_reflections(images, n_usable_patterns);
 
 	/* Make initial estimates */
 	STATUS("Performing initial scaling.\n");
 	select_scalable_reflections(images, n_total_patterns);
-	full = scale_intensities(images, n_total_patterns, sym, obs, cref);
+	full = scale_intensities(images, n_usable_patterns, sym, obs, cref);
 
 	/* Iterate */
 	for ( i=0; i<n_iter; i++ ) {
@@ -446,8 +401,8 @@ int main(int argc, char *argv[])
 
 		/* Re-estimate all the full intensities */
 		reflist_free(full);
-		select_scalable_reflections(images, n_total_patterns);
-		full = scale_intensities(images, n_total_patterns,
+		select_scalable_reflections(images, n_usable_patterns);
+		full = scale_intensities(images, n_usable_patterns,
 		                         sym, obs, cref);
 
 		fclose(fhg);
@@ -456,7 +411,7 @@ int main(int argc, char *argv[])
 	}
 
 	STATUS("Final scale factors:\n");
-	for ( i=0; i<n_total_patterns; i++ ) {
+	for ( i=0; i<n_usable_patterns; i++ ) {
 		STATUS("%4i : %5.2f\n", i, images[i].osf);
 	}
 
@@ -464,18 +419,18 @@ int main(int argc, char *argv[])
 	write_reflist(outfile, full, images[0].indexed_cell);
 
 	/* Clean up */
-	for ( i=0; i<n_total_patterns; i++ ) {
+	for ( i=0; i<n_usable_patterns; i++ ) {
 		reflist_free(images[i].reflections);
+		reflist_free(images[i].raw_reflections);
 	}
 	reflist_free(full);
 	delete_items(obs);
 	free(sym);
 	free(outfile);
-	free(det->panels);
-	free(det);
+	free_detector_geometry(det);
 	free(beam);
 	free(cref);
-	for ( i=0; i<n_total_patterns; i++ ) {
+	for ( i=0; i<n_usable_patterns; i++ ) {
 		cell_free(images[i].indexed_cell);
 		free(images[i].filename);
 	}
diff --git a/src/reflist-utils.c b/src/reflist-utils.c
index ea590aaf..0a5c4863 100644
--- a/src/reflist-utils.c
+++ b/src/reflist-utils.c
@@ -11,6 +11,7 @@
 
 
 #include <stdio.h>
+#include <assert.h>
 
 
 #include "reflist.h"
@@ -359,3 +360,40 @@ RefList *read_reflections(const char *filename)
 
 	return out;
 }
+
+
+RefList *asymmetric_indices(RefList *in, const char *sym, ReflItemList *obs)
+{
+	Reflection *refl;
+	RefListIterator *iter;
+	RefList *new;
+
+	new = reflist_new();
+
+	for ( refl = first_refl(in, &iter);
+	      refl != NULL;
+	      refl = next_refl(refl, iter) ) {
+
+		signed int h, k, l;
+		signed int ha, ka, la;
+		Reflection *cr;
+
+		get_indices(refl, &h, &k, &l);
+
+		get_asymm(h, k, l, &ha, &ka, &la, sym);
+
+		cr = add_refl(new, ha, ka, la);
+		assert(cr != NULL);
+
+		copy_data(cr, refl);
+
+		if ( obs != NULL ) {
+			if ( !find_item(obs, ha, ka, la) ) {
+				add_item(obs, ha, ka, la);
+			}
+		}
+
+	}
+
+	return new;
+}
diff --git a/src/reflist-utils.h b/src/reflist-utils.h
index b7ecd3ea..2efa3a76 100644
--- a/src/reflist-utils.h
+++ b/src/reflist-utils.h
@@ -41,4 +41,7 @@ extern int find_equiv_in_list(RefList *list, signed int h, signed int k,
                               signed int l, const char *sym, signed int *hu,
                               signed int *ku, signed int *lu);
 
+extern RefList *asymmetric_indices(RefList *in, const char *sym,
+                                   ReflItemList *obs);
+
 #endif	/* REFLIST_UTILS_H */