Get rid of "image.raw_reflections"

9 files changed, 130 insertions, 90 deletions

diff --git a/src/geometry.c b/src/geometry.c
index db261e6b..d5845e5d 100644
--- a/src/geometry.c
+++ b/src/geometry.c
@@ -24,6 +24,7 @@
 #include "beam-parameters.h"
 #include "reflist.h"
 #include "reflist-utils.h"
+#include "symmetry.h"
 
 
 static signed int locate_peak(double x, double y, double z, double k,
@@ -293,10 +294,8 @@ double integrate_all(struct image *image, RefList *reflections)
 
 
 /* Decide which reflections can be scaled */
-static int select_scalable_reflections(RefList *list)
+static void select_scalable_reflections(RefList *list, ReflItemList *sc_l)
 {
-	int n_scalable = 0;
-
 	Reflection *refl;
 	RefListIterator *iter;
 
@@ -312,16 +311,27 @@ static int select_scalable_reflections(RefList *list)
 		if ( v < 0.1 ) scalable = 0;
 
 		set_scalable(refl, scalable);
-		if ( scalable ) n_scalable++;
+		if ( scalable && (sc_l != NULL) ) {
 
-	}
+			signed int h, k, l;
+
+			get_indices(refl, &h, &k, &l);  /* Should already be
+			                                 * asymmetric */
+			if ( !find_item(sc_l, h, k, l) ) {
+				add_item(sc_l, h, k, l);
+			}
+
+		}
 
-	return n_scalable;
+	}
 }
 
 
-/* Calculate partialities and apply them to the image's raw_reflections */
+/* Calculate partialities and apply them to the image's raw_reflections,
+ * returning a ReflItemList of the currentl scalable (asymmetric) reflections.
+ */
 void update_partialities(struct image *image, const char *sym,
+                         ReflItemList *scalable,
                          int *n_expected, int *n_found, int *n_notfound)
 {
 	Reflection *refl;
@@ -335,56 +345,54 @@ void update_partialities(struct image *image, const char *sym,
 	      refl = next_refl(refl, iter) )
 	{
 
-		Reflection *peak_in_pattern;
+		Reflection *p_peak;
 		double r1, r2, p, x, y;
 		signed int h, k, l;
+		signed int ha, ka, la;
 		int clamp1, clamp2;
+		int found = 0;
 
 		/* Get predicted indices and location */
 		get_indices(refl, &h, &k, &l);
 		get_detector_pos(refl, &x, &y);
 		if ( n_expected != NULL ) (*n_expected)++;
 
+		/* Get the asymmetric indices, with which the reflection in the
+		 * image's list will be indexed. */
+		get_asymm(h, k, l, &ha, &ka, &la, sym);
+
 		/* 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)++;
+		p_peak = find_refl(image->reflections, ha, ka, la);
+		do {
+
+			signed int hs, ks, ls;
+
+			if ( p_peak != NULL ) {
+				get_symmetric_indices(p_peak, &hs, &ks, &ls);
+				if ( (hs==h) && (ks==k) && (ls==l) ) found = 1;
+			}
+
+			if ( !found && (p_peak != NULL ) ) {
+				p_peak = next_found_refl(p_peak);
+			}
+
+		} while ( !found && (p_peak != NULL) );
+		if ( !found ) {
+			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);
+		set_partial(p_peak, r1, r2, p, clamp1, clamp2);
 
 		/* Transfer detector location */
 		get_detector_pos(refl, &x, &y);
-		set_detector_pos(peak_in_pattern, 0.0, x, y);
+		set_detector_pos(p_peak, 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);
-	select_scalable_reflections(image->reflections);
-
-	/* Need these lists to work fast */
-	optimise_reflist(image->reflections);
+	select_scalable_reflections(image->reflections, scalable);
 }
diff --git a/src/geometry.h b/src/geometry.h
index f8c007fa..171368a7 100644
--- a/src/geometry.h
+++ b/src/geometry.h
@@ -24,8 +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);
+extern void update_partialities(struct image *image, const char *sym,
+                                ReflItemList *scalable,
+                                int *n_expected, int *n_found, int *n_notfound);
+
 #endif	/* GEOMETRY_H */
diff --git a/src/image.h b/src/image.h
index 22e0ddae..d6e22aee 100644
--- a/src/image.h
+++ b/src/image.h
@@ -87,7 +87,6 @@ typedef struct _imagefeaturelist ImageFeatureList;
  *    int                     height;
  *
  *    RefList                 *reflections;
- *    RefList                 *raw_reflections;
  *
  *    ImageFeatureList        *features;
  * };
@@ -151,9 +150,6 @@ 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 9f931496..dc29d0d6 100644
--- a/src/partialator.c
+++ b/src/partialator.c
@@ -155,7 +155,7 @@ int main(int argc, char *argv[])
 	FILE *fh;
 	int nthreads = 1;
 	struct detector *det;
-	ReflItemList *obs;
+	ReflItemList *scalable;
 	int i;
 	int n_total_patterns;
 	struct image *images;
@@ -282,9 +282,11 @@ int main(int argc, char *argv[])
 
 	/* Fill in what we know about the images so far */
 	rewind(fh);
-	obs = new_items();
+	scalable = new_items();
 	for ( i=0; i<n_total_patterns; i++ ) {
 
+		RefList *as;
+
 		images[n_usable_patterns].det = NULL;
 
 		if ( read_chunk(fh, &images[n_usable_patterns]) != 0 ) {
@@ -316,13 +318,14 @@ int main(int argc, char *argv[])
 		images[n_usable_patterns].beam = NULL;
 
 		/* This is the raw list of reflections */
-		images[n_usable_patterns].raw_reflections =
-		                          images[n_usable_patterns].reflections;
-		images[n_usable_patterns].reflections = NULL;
+		as = asymmetric_indices(images[n_usable_patterns].reflections,
+		                        sym);
+		optimise_reflist(as);
+		reflist_free(images[n_usable_patterns].reflections);
+		images[n_usable_patterns].reflections = as;
 
-		update_partialities_and_asymm(&images[n_usable_patterns], sym,
-		                              obs, &n_expected, &n_found,
-		                              &n_notfound);
+		update_partialities(&images[n_usable_patterns], sym, scalable,
+		                    &n_expected, &n_found, &n_notfound);
 
 		progress_bar(i, n_total_patterns-1, "Loading pattern data");
 		n_usable_patterns++;
@@ -331,14 +334,15 @@ int main(int argc, char *argv[])
 	fclose(fh);
 	STATUS("Found %5.2f%% of the expected peaks (missed %i of %i).\n",
 	       100.0 * (double)n_found / n_expected, n_notfound, n_expected);
-	STATUS("Mean measurements per unique reflection: %5.2f\n",
-	       (double)n_found / num_items(obs));
+	STATUS("Mean measurements per scalable unique reflection: %5.2f\n",
+	       (double)n_found / num_items(scalable));
 
 	cref = find_common_reflections(images, n_usable_patterns);
 
 	/* Make initial estimates */
 	STATUS("Performing initial scaling.\n");
-	full = scale_intensities(images, n_usable_patterns, sym, obs, cref);
+	full = scale_intensities(images, n_usable_patterns, sym,
+	                         scalable, cref);
 
 	/* Iterate */
 	for ( i=0; i<n_iter; i++ ) {
@@ -364,13 +368,13 @@ int main(int argc, char *argv[])
 		}
 
 		/* Refine the geometry of all patterns to get the best fit */
-		refine_all(images, n_total_patterns, det, sym, obs, full,
+		refine_all(images, n_total_patterns, det, sym, scalable, full,
 		           nthreads, fhg, fhp);
 
 		/* Re-estimate all the full intensities */
 		reflist_free(full);
 		full = scale_intensities(images, n_usable_patterns,
-		                         sym, obs, cref);
+		                         sym, scalable, cref);
 
 		fclose(fhg);
 		fclose(fhp);
@@ -388,10 +392,9 @@ int main(int argc, char *argv[])
 	/* Clean up */
 	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);
+	delete_items(scalable);
 	free(sym);
 	free(outfile);
 	free_detector_geometry(det);
diff --git a/src/post-refinement.c b/src/post-refinement.c
index d804dbc2..1fd4d13b 100644
--- a/src/post-refinement.c
+++ b/src/post-refinement.c
@@ -95,21 +95,21 @@ static double gradient(struct image *image, int k, Reflection *refl, double r)
 	double bsx, bsy, bsz;
 	double csx, csy, csz;
 	double xl, yl, zl;
-	signed int hi, ki, li;
+	signed int hs, ks, ls;
 	double r1, r2, p;
 	int clamp_low, clamp_high;
 
-	get_indices(refl, &hi, &ki, &li);
+	get_symmetric_indices(refl, &hs, &ks, &ls);
 
 	cell_get_reciprocal(image->indexed_cell, &asx, &asy, &asz,
 	                                         &bsx, &bsy, &bsz,
 	                                         &csx, &csy, &csz);
-	xl = hi*asx + ki*bsx + li*csx;
-	yl = hi*asy + ki*bsy + li*csy;
-	zl = hi*asz + ki*bsz + li*csz;
+	xl = hs*asx + ks*bsx + ls*csx;
+	yl = hs*asy + ks*bsy + ls*csy;
+	zl = hs*asz + ks*bsz + ls*csz;
 
-	ds = 2.0 * resolution(image->indexed_cell, hi, ki, li);
-	tt = angle_between(0.0, 0.0, 1.0,  xl, yl, zl+k);
+	ds = 2.0 * resolution(image->indexed_cell, hs, ks, ls);
+	tt = angle_between(0.0, 0.0, 1.0, xl, yl, zl+1.0/image->lambda);
 	azi = angle_between(1.0, 0.0, 0.0, xl, yl, 0.0);
 
 	get_partial(refl, &r1, &r2, &p, &clamp_low, &clamp_high);
@@ -142,23 +142,23 @@ static double gradient(struct image *image, int k, Reflection *refl, double r)
 
 	/* Cell parameters and orientation */
 	case REF_ASX :
-		return hi * sin(tt) * cos(azi) * g;
+		return hs * sin(tt) * cos(azi) * g;
 	case REF_BSX :
-		return ki * sin(tt) * g;
+		return ks * sin(tt) * g;
 	case REF_CSX :
-		return li * sin(tt) * g;
+		return ls * sin(tt) * g;
 	case REF_ASY :
-		return hi * sin(tt) * g;
+		return hs * sin(tt) * g;
 	case REF_BSY :
-		return ki * sin(tt) * g;
+		return ks * sin(tt) * g;
 	case REF_CSY :
-		return li * sin(tt) * g;
+		return ls * sin(tt) * g;
 	case REF_ASZ :
-		return hi * cos(tt) * g;
+		return hs * cos(tt) * g;
 	case REF_BSZ :
-		return ki * cos(tt) * g;
+		return ks * cos(tt) * g;
 	case REF_CSZ :
-		return li * cos(tt) * g;
+		return ls * cos(tt) * g;
 
 	}
 
@@ -285,8 +285,7 @@ static double pr_iterate(struct image *image, const RefList *full,
 		/* Calculate all gradients for this reflection */
 		for ( k=0; k<NUM_PARAMS; k++ ) {
 			double gr;
-			gr = gradient(image, k, refl,
-			              image->profile_radius);
+			gr = gradient(image, k, refl, image->profile_radius);
 			gradients[k] = gr;
 		}
 
@@ -403,8 +402,7 @@ static void plot_curve(struct image *image, const RefList *full,
 		ax = origval + (double)i*shval;
 		cell_set_reciprocal(cell, ax, ay, az, bx, by, bz, cx, cy, cz);
 
-		update_partialities_and_asymm(image, sym,
-		                              NULL, NULL, NULL, NULL);
+		update_partialities(image, sym, NULL, NULL, NULL, NULL);
 
 		dev = mean_partial_dev(image, full, sym);
 		STATUS("%i %e %e\n", i, ax, dev);
@@ -431,8 +429,7 @@ void pr_refine(struct image *image, const RefList *full, const char *sym)
 		double dev;
 
 		max_shift = pr_iterate(image, full, sym);
-		update_partialities_and_asymm(image, sym,
-		                              NULL, NULL, NULL, NULL);
+		update_partialities(image, sym, NULL, NULL, NULL, NULL);
 
 		dev = mean_partial_dev(image, full, sym);
 		STATUS("PR Iteration %2i: max shift = %5.2f dev = %5.2f\n",
diff --git a/src/reflist-utils.c b/src/reflist-utils.c
index 0a5c4863..016fa317 100644
--- a/src/reflist-utils.c
+++ b/src/reflist-utils.c
@@ -362,7 +362,7 @@ RefList *read_reflections(const char *filename)
 }
 
 
-RefList *asymmetric_indices(RefList *in, const char *sym, ReflItemList *obs)
+RefList *asymmetric_indices(RefList *in, const char *sym)
 {
 	Reflection *refl;
 	RefListIterator *iter;
@@ -386,12 +386,7 @@ RefList *asymmetric_indices(RefList *in, const char *sym, ReflItemList *obs)
 		assert(cr != NULL);
 
 		copy_data(cr, refl);
-
-		if ( obs != NULL ) {
-			if ( !find_item(obs, ha, ka, la) ) {
-				add_item(obs, ha, ka, la);
-			}
-		}
+		set_symmetric_indices(cr, h, k, l);
 
 	}
 
diff --git a/src/reflist-utils.h b/src/reflist-utils.h
index 2efa3a76..775e375a 100644
--- a/src/reflist-utils.h
+++ b/src/reflist-utils.h
@@ -41,7 +41,6 @@ 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);
+extern RefList *asymmetric_indices(RefList *in, const char *sym);
 
 #endif	/* REFLIST_UTILS_H */
diff --git a/src/reflist.c b/src/reflist.c
index e420b736..d23b952c 100644
--- a/src/reflist.c
+++ b/src/reflist.c
@@ -45,6 +45,11 @@
 
 struct _refldata {
 
+	/* Symmetric indices (i.e. the "real" indices) */
+	signed int hs;
+	signed int ks;
+	signed int ls;
+
 	/* Partiality and related geometrical stuff */
 	double r1;  /* First excitation error */
 	double r2;  /* Second excitation error */
@@ -288,6 +293,29 @@ void get_indices(const Reflection *refl,
 
 
 /**
+ * get_symmetric_indices:
+ * @refl: A %Reflection
+ * @h: Location at which to store the 'h' index of the reflection
+ * @k: Location at which to store the 'k' index of the reflection
+ * @l: Location at which to store the 'l' index of the reflection
+ *
+ * This function gives the symmetric indices, that is, the "real" indices before
+ * squashing down to the asymmetric reciprocal unit.  This may be useful if the
+ * list is indexed according to the asymmetric indices, but you still need
+ * access to the symmetric version.  This happens during post-refinement.
+ *
+ **/
+void get_symmetric_indices(const Reflection *refl,
+                                  signed int *hs, signed int *ks,
+                                  signed int *ls)
+{
+	*hs = refl->data.hs;
+	*ks = refl->data.ks;
+	*ls = refl->data.ls;
+}
+
+
+/**
  * get_partiality:
  * @refl: A %Reflection
  *
@@ -489,6 +517,15 @@ void set_ph(Reflection *refl, double phase)
 }
 
 
+void set_symmetric_indices(Reflection *refl,
+                           signed int hs, signed int ks, signed int ls)
+{
+	refl->data.hs = hs;
+	refl->data.ks = ks;
+	refl->data.ls = ls;
+}
+
+
 /********************************* Insertion **********************************/
 
 static void insert_node(Reflection *head, Reflection *new)
diff --git a/src/reflist.h b/src/reflist.h
index f380ba85..580ea33c 100644
--- a/src/reflist.h
+++ b/src/reflist.h
@@ -47,6 +47,9 @@ extern void get_detector_pos(const Reflection *refl, double *fs, double *ss);
 extern double get_partiality(const Reflection *refl);
 extern void get_indices(const Reflection *refl,
                         signed int *h, signed int *k, signed int *l);
+extern void get_symmetric_indices(const Reflection *refl,
+                                  signed int *hs, signed int *ks,
+                                  signed int *ls);
 extern double get_intensity(const Reflection *refl);
 extern void get_partial(const Reflection *refl, double *r1, double *r2,
                         double *p, int *clamp_low, int *clamp_high);
@@ -68,6 +71,8 @@ extern void set_redundancy(Reflection *refl, int red);
 extern void set_sum_squared_dev(Reflection *refl, double dev);
 extern void set_esd_intensity(Reflection *refl, double esd);
 extern void set_ph(Reflection *refl, double phase);
+extern void set_symmetric_indices(Reflection *refl,
+                                  signed int hs, signed int ks, signed int ls);
 
 /* Insertion */
 extern Reflection *add_refl(RefList *list,