partialator: Save reflection list

7 files changed, 66 insertions, 57 deletions

diff --git a/src/hrs-scaling.c b/src/hrs-scaling.c
index e2adaaf6..03458670 100644
--- a/src/hrs-scaling.c
+++ b/src/hrs-scaling.c
@@ -311,19 +311,20 @@ static double iterate_scale(struct image *images, int n,
 }
 
 
-static double *lsq_intensities(struct image *images, int n,
-                               ReflItemList *obs, const char *sym)
+static RefList *lsq_intensities(struct image *images, int n,
+                                ReflItemList *obs, const char *sym)
 {
-	double *I_full;
+	RefList *full;
 	int i;
 
-	I_full = new_list_intensity();
+	full = reflist_new();
 	for ( i=0; i<num_items(obs); i++ ) {
 
 		struct refl_item *it = get_item(obs, i);
 		double num = 0.0;
 		double den = 0.0;
 		int m;
+		Reflection *new;
 
 		/* For each frame */
 		for ( m=0; m<n; m++ ) {
@@ -351,11 +352,12 @@ static double *lsq_intensities(struct image *images, int n,
 
 		}
 
-		set_intensity(I_full, it->h, it->k, it->l, num/den);
+		new = add_refl(full, it->h, it->k, it->l);
+		set_int(new, num/den);
 
 	}
 
-	return I_full;
+	return full;
 }
 
 
@@ -377,11 +379,11 @@ 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, char *cref)
+RefList *scale_intensities(struct image *images, int n, const char *sym,
+                           ReflItemList *obs, char *cref)
 {
 	int m;
-	double *I_full;
+	RefList *full;
 	int i;
 	double max_shift;
 
@@ -399,6 +401,6 @@ double *scale_intensities(struct image *images, int n, const char *sym,
 
 	} while ( (max_shift > 0.01) && (i < MAX_CYCLES) );
 
-	I_full = lsq_intensities(images, n, obs, sym);
-	return I_full;
+	full = lsq_intensities(images, n, obs, sym);
+	return full;
 }
diff --git a/src/hrs-scaling.h b/src/hrs-scaling.h
index 04a321e3..d65129b8 100644
--- a/src/hrs-scaling.h
+++ b/src/hrs-scaling.h
@@ -20,8 +20,8 @@
 
 #include "image.h"
 
-extern double *scale_intensities(struct image *images, int n, const char *sym,
-                                 ReflItemList *obs, char *cref);
+extern RefList *scale_intensities(struct image *images, int n, const char *sym,
+                                  ReflItemList *obs, char *cref);
 
 extern char *find_common_reflections(struct image *images, int n);
 
diff --git a/src/partialator.c b/src/partialator.c
index 553a216a..700159ee 100644
--- a/src/partialator.c
+++ b/src/partialator.c
@@ -64,7 +64,7 @@ struct refine_args
 {
 	const char *sym;
 	ReflItemList *obs;
-	double *i_full;
+	RefList *full;
 	struct image *image;
 	FILE *graph;
 	FILE *pgraph;
@@ -77,13 +77,13 @@ static void refine_image(int mytask, void *tasks)
 	struct refine_args *pargs = &all_args[mytask];
 	struct image *image = pargs->image;
 
-	pr_refine(image, pargs->i_full, pargs->sym);
+	pr_refine(image, pargs->full, pargs->sym);
 }
 
 
 static void refine_all(struct image *images, int n_total_patterns,
                        struct detector *det, const char *sym,
-                       ReflItemList *obs, double *i_full, int nthreads,
+                       ReflItemList *obs, RefList *full, int nthreads,
                        FILE *graph, FILE *pgraph)
 {
 	struct refine_args *tasks;
@@ -94,7 +94,7 @@ static void refine_all(struct image *images, int n_total_patterns,
 
 		tasks[i].sym = sym;
 		tasks[i].obs = obs;
-		tasks[i].i_full = i_full;
+		tasks[i].full = full;
 		tasks[i].image = &images[i];
 		tasks[i].graph = graph;
 		tasks[i].pgraph = pgraph;
@@ -157,7 +157,7 @@ int main(int argc, char *argv[])
 	struct image *images;
 	int n_iter = 10;
 	struct beam_params *beam = NULL;
-	double *I_full;
+	RefList *full;
 	int n_found = 0;
 	int n_expected = 0;
 	int n_notfound = 0;
@@ -286,6 +286,9 @@ int main(int argc, char *argv[])
 			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;
@@ -372,7 +375,7 @@ int main(int argc, char *argv[])
 	/* 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, cref);
+	full = scale_intensities(images, n_total_patterns, sym, obs, cref);
 
 	/* Iterate */
 	for ( i=0; i<n_iter; i++ ) {
@@ -398,14 +401,14 @@ 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, I_full,
+		refine_all(images, n_total_patterns, det, sym, obs, full,
 		           nthreads, fhg, fhp);
 
 		/* Re-estimate all the full intensities */
-		free(I_full);
+		reflist_free(full);
 		select_scalable_reflections(images, n_total_patterns);
-		I_full = scale_intensities(images, n_total_patterns,
-		                           sym, obs, cref);
+		full = scale_intensities(images, n_total_patterns,
+		                         sym, obs, cref);
 
 		fclose(fhg);
 		fclose(fhp);
@@ -418,13 +421,13 @@ int main(int argc, char *argv[])
 	}
 
 	/* Output results */
-	//write_reflist(outfile, obs, I_full, NULL, NULL, cts, NULL);
+	write_reflist(outfile, full, images[0].indexed_cell);
 
 	/* Clean up */
 	for ( i=0; i<n_total_patterns; i++ ) {
 		reflist_free(images[i].reflections);
 	}
-	free(I_full);
+	reflist_free(full);
 	delete_items(obs);
 	free(sym);
 	free(outfile);
diff --git a/src/post-refinement.c b/src/post-refinement.c
index 0bd26f76..d5f9d3f8 100644
--- a/src/post-refinement.c
+++ b/src/post-refinement.c
@@ -248,8 +248,8 @@ static void apply_shift(struct image *image, int k, double shift)
 }
 
 
-/* Perform one cycle of post refinement on 'image' against 'i_full' */
-static double pr_iterate(struct image *image, const double *i_full,
+/* Perform one cycle of post refinement on 'image' against 'full' */
+static double pr_iterate(struct image *image, const RefList *full,
                          const char *sym)
 {
 	gsl_matrix *M;
@@ -277,6 +277,7 @@ static double pr_iterate(struct image *image, const double *i_full,
 		double I_partial;
 		int k;
 		double p;
+		Reflection *match;
 
 		get_indices(refl, &hind, &kind, &lind);
 
@@ -286,7 +287,9 @@ static double pr_iterate(struct image *image, const double *i_full,
 		I_partial = get_intensity(refl);
 
 		get_asymm(hind, kind, lind, &ha, &ka, &la, sym);
-		I_full = lookup_intensity(i_full, ha, ka, la);
+		match = find_refl(full, ha, ka, la);
+		assert(match != NULL);
+		I_full = get_intensity(match);
 		p = get_partiality(refl);
 		delta_I = I_partial - (p * I_full / image->osf);
 
@@ -338,14 +341,14 @@ static double pr_iterate(struct image *image, const double *i_full,
 }
 
 
-void pr_refine(struct image *image, const double *i_full, const char *sym)
+void pr_refine(struct image *image, const RefList *full, const char *sym)
 {
 	double max_shift;
 	int i;
 
 	i = 0;
 	do {
-		max_shift = pr_iterate(image, i_full, sym);
+		max_shift = pr_iterate(image, full, sym);
 		STATUS("Iteration %2i: max shift = %5.2f\n", i, max_shift);
 		i++;
 	} while ( (max_shift > 0.01) && (i < MAX_CYCLES) );
diff --git a/src/post-refinement.h b/src/post-refinement.h
index a10fafe4..3c8d8587 100644
--- a/src/post-refinement.h
+++ b/src/post-refinement.h
@@ -24,7 +24,7 @@
 #include "utils.h"
 
 
-extern void pr_refine(struct image *image, const double *i_full,
+extern void pr_refine(struct image *image, const RefList *full,
                       const char *sym);
 
 
diff --git a/src/reflist.c b/src/reflist.c
index a81c738e..30d9f696 100644
--- a/src/reflist.c
+++ b/src/reflist.c
@@ -133,7 +133,7 @@ void reflist_free(RefList *list)
 /********************************** Search ************************************/
 
 /* Return the first reflection in 'list' with the given indices, or NULL */
-Reflection *find_refl(RefList *list, INDICES)
+Reflection *find_refl(const RefList *list, INDICES)
 {
 	unsigned int search = SERIAL(h, k, l);
 	Reflection *refl = list->head->child[0];
@@ -185,20 +185,21 @@ Reflection *next_found_refl(Reflection *refl)
 
 /********************************** Getters ***********************************/
 
-double get_excitation_error(Reflection *refl)
+double get_excitation_error(const Reflection *refl)
 {
 	return refl->data.excitation_error;
 }
 
 
-void get_detector_pos(Reflection *refl, double *x, double *y)
+void get_detector_pos(const Reflection *refl, double *x, double *y)
 {
 	*x = refl->data.x;
 	*y = refl->data.y;
 }
 
 
-void get_indices(Reflection *refl, signed int *h, signed int *k, signed int *l)
+void get_indices(const Reflection *refl,
+                 signed int *h, signed int *k, signed int *l)
 {
 	*h = refl->data.h;
 	*k = refl->data.k;
@@ -206,19 +207,19 @@ void get_indices(Reflection *refl, signed int *h, signed int *k, signed int *l)
 }
 
 
-double get_partiality(Reflection *refl)
+double get_partiality(const Reflection *refl)
 {
 	return refl->data.p;
 }
 
 
-double get_intensity(Reflection *refl)
+double get_intensity(const Reflection *refl)
 {
 	return refl->data.intensity;
 }
 
 
-void get_partial(Reflection *refl, double *r1, double *r2, double *p,
+void get_partial(const Reflection *refl, double *r1, double *r2, double *p,
                  int *clamp_low, int *clamp_high)
 {
 	*r1 = refl->data.r1;
@@ -229,31 +230,31 @@ void get_partial(Reflection *refl, double *r1, double *r2, double *p,
 }
 
 
-int get_scalable(Reflection *refl)
+int get_scalable(const Reflection *refl)
 {
 	return refl->data.scalable;
 }
 
 
-int get_redundancy(Reflection *refl)
+int get_redundancy(const Reflection *refl)
 {
 	return refl->data.redundancy;
 }
 
 
-double get_sum_squared_dev(Reflection *refl)
+double get_sum_squared_dev(const Reflection *refl)
 {
 	return refl->data.sum_squared_dev;
 }
 
 
-double get_esd_intensity(Reflection *refl)
+double get_esd_intensity(const Reflection *refl)
 {
 	return refl->data.esd_i;
 }
 
 
-double get_phase(Reflection *refl)
+double get_phase(const Reflection *refl)
 {
 	return refl->data.phase;
 }
@@ -261,7 +262,7 @@ double get_phase(Reflection *refl)
 
 /********************************** Setters ***********************************/
 
-void copy_data(Reflection *to, Reflection *from)
+void copy_data(Reflection *to, const Reflection *from)
 {
 	memcpy(&to->data, &from->data, sizeof(struct _refldata));
 }
diff --git a/src/reflist.h b/src/reflist.h
index 5b647102..462a19af 100644
--- a/src/reflist.h
+++ b/src/reflist.h
@@ -28,26 +28,26 @@ extern RefList *reflist_new(void);
 extern void reflist_free(RefList *list);
 
 /* Search */
-extern Reflection *find_refl(RefList *list, INDICES);
+extern Reflection *find_refl(const RefList *list, INDICES);
 extern Reflection *next_found_refl(Reflection *refl);
 
 /* Get */
-extern double get_excitation_error(Reflection *refl);
-extern void get_detector_pos(Reflection *refl, double *x, double *y);
-extern void get_indices(Reflection *refl,
+extern double get_excitation_error(const Reflection *refl);
+extern void get_detector_pos(const Reflection *refl, double *x, double *y);
+extern void get_indices(const Reflection *refl,
                         signed int *h, signed int *k, signed int *l);
-extern double get_partiality(Reflection *refl);
-extern double get_intensity(Reflection *refl);
-extern void get_partial(Reflection *refl, double *r1, double *r2, double *p,
+extern double get_partiality(const Reflection *refl);
+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);
-extern int get_scalable(Reflection *refl);
-extern int get_redundancy(Reflection *refl);
-extern double get_sum_squared_dev(Reflection *refl);
-extern double get_esd_intensity(Reflection *refl);
-extern double get_phase(Reflection *refl);
+extern int get_scalable(const Reflection *refl);
+extern int get_redundancy(const Reflection *refl);
+extern double get_sum_squared_dev(const Reflection *refl);
+extern double get_esd_intensity(const Reflection *refl);
+extern double get_phase(const Reflection *refl);
 
 /* Set */
-extern void copy_data(Reflection *to, Reflection *from);
+extern void copy_data(Reflection *to, const Reflection *from);
 extern void set_detector_pos(Reflection *refl, double exerr,
                              double x, double y);
 extern void set_partial(Reflection *refl, double r1, double r2, double p,