Don't constantly re-integrate peaks

5 files changed, 119 insertions, 290 deletions

diff --git a/src/partialator.c b/src/partialator.c
index a4f3381c..37c8fbd6 100644
--- a/src/partialator.c
+++ b/src/partialator.c
@@ -37,10 +37,6 @@
 #include "reflist.h"
 
 
-/* Maximum number of iterations of NLSq to do for each image per macrocycle. */
-#define MAX_CYCLES (100)
-
-
 static void show_help(const char *s)
 {
 	printf("Syntax: %s [options]\n\n", s);
@@ -57,9 +53,6 @@ static void show_help(const char *s)
 "                              initial values for parameters, and nominal\n"
 "                              wavelengths if no per-shot value is found in \n"
 "                              an HDF5 file.\n"
-"  -x, --prefix=<p>           Prefix filenames from input file with <p>.\n"
-"      --basename             Remove the directory parts of the filenames.\n"
-"      --no-check-prefix      Don't attempt to correct the --prefix.\n"
 "  -y, --symmetry=<sym>       Merge according to symmetry <sym>.\n"
 "  -n, --iterations=<n>       Run <n> cycles of scaling and post-refinement.\n"
 "\n"
@@ -83,58 +76,8 @@ static void refine_image(int mytask, void *tasks)
 	struct refine_args *all_args = tasks;
 	struct refine_args *pargs = &all_args[mytask];
 	struct image *image = pargs->image;
-	double nominal_photon_energy = pargs->image->beam->photon_energy;
-	struct hdfile *hdfile;
-	int i;
-	double dev, last_dev;
-	RefList *reflections;
-
-	hdfile = hdfile_open(image->filename);
-	if ( hdfile == NULL ) {
-		ERROR("Couldn't open '%s'\n", image->filename);
-		return;
-	} else if ( hdfile_set_image(hdfile, "/data/data0") ) {
-		ERROR("Couldn't select path\n");
-		hdfile_close(hdfile);
-		return;
-	}
-
-	if ( hdf5_read(hdfile, pargs->image, 0, nominal_photon_energy) ) {
-		ERROR("Couldn't read '%s'\n", image->filename);
-		hdfile_close(hdfile);
-		return;
-	}
-
-	double a, b, c, al, be, ga;
-	cell_get_parameters(image->indexed_cell, &a, &b, &c, &al, &be, &ga);
-	STATUS("Initial cell: %5.2f %5.2f %5.2f nm %5.2f %5.2f %5.2f deg\n",
-	       a/1.0e-9, b/1.0e-9, c/1.0e-9,
-	       rad2deg(al), rad2deg(be), rad2deg(ga));
-
-	/* FIXME: Don't do this each time */
-	reflections = find_intersections(image, image->indexed_cell, 0);
-	dev = +INFINITY;
-	i = 0;
-	do {
-		last_dev = dev;
-		dev = pr_iterate(image, pargs->i_full, pargs->sym, reflections);
-		STATUS("Iteration %2i: mean dev = %5.2f\n", i, dev);
-		i++;
-	} while ( (fabs(last_dev - dev) > 1.0) && (i < MAX_CYCLES) );
-	mean_partial_dev(image, reflections, pargs->sym,
-	                 pargs->i_full, pargs->graph);
-	if ( pargs->pgraph ) {
-		fprintf(pargs->pgraph, "%5i %5.2f\n", mytask, dev);
-	}
 
-	free(image->data);
-	if ( image->flags != NULL ) free(image->flags);
-	hdfile_close(hdfile);
-	reflist_free(reflections);
-
-	/* Muppet proofing */
-	image->data = NULL;
-	image->flags = NULL;
+	pr_refine(image, pargs->i_full, pargs->sym);
 }
 
 
@@ -165,87 +108,6 @@ static void refine_all(struct image *images, int n_total_patterns,
 }
 
 
-static void uniquify(Reflection *refl, const char *sym)
-{
-	signed int h, k, l;
-	signed int ha, ka, la;
-
-	get_indices(refl, &h, &k, &l);
-	get_asymm(h, k, l, &ha, &ka, &la, sym);
-	set_indices(refl, h, k, l);
-}
-
-
-/* FIXME: Get rid of this */
-static void integrate_image(struct image *image, ReflItemList *obs,
-                            const char *sym)
-{
-	RefList *reflections;
-	Reflection *refl;
-	struct hdfile *hdfile;
-	double nominal_photon_energy = image->beam->photon_energy;
-
-	hdfile = hdfile_open(image->filename);
-	if ( hdfile == NULL ) {
-		ERROR("Couldn't open '%s'\n", image->filename);
-		return;
-	} else if ( hdfile_set_image(hdfile, "/data/data0") ) {
-		ERROR("Couldn't select path\n");
-		hdfile_close(hdfile);
-		return;
-	}
-
-	if ( hdf5_read(hdfile, image, 0, nominal_photon_energy) ) {
-		ERROR("Couldn't read '%s'\n", image->filename);
-		hdfile_close(hdfile);
-		return;
-	}
-
-	/* Figure out which spots should appear in this pattern */
-	reflections = find_intersections(image, image->indexed_cell, 0);
-
-	/* For each reflection, estimate the partiality */
-	for ( refl = first_refl(reflections);
-	      refl != NULL;
-	      refl = next_refl(refl) ) {
-
-		signed int h, k, l;
-		float i_partial;
-		float xc, yc;
-		double x, y;
-
-		uniquify(refl, sym);
-		get_indices(refl, &h, &k, &l);
-
-		/* Don't attempt to use spots with very small
-		 * partialities, since it won't be accurate. */
-		if ( get_partiality(refl) < 0.1 ) continue;
-
-		/* Actual measurement of this reflection from this pattern? */
-		get_detector_pos(refl, &x, &y);
-		if ( integrate_peak(image, x, y,
-		                    &xc, &yc, &i_partial, NULL, NULL, 1, 0) ) {
-			delete_refl(refl);
-			continue;
-		}
-
-		set_int(refl, i_partial);
-
-		if ( !find_item(obs, h, k, l) ) add_item(obs, h, k, l);
-
-	}
-	image->reflections = reflections;
-
-	free(image->data);
-	if ( image->flags != NULL ) free(image->flags);
-	hdfile_close(hdfile);
-
-	/* Muppet proofing */
-	image->data = NULL;
-	image->flags = NULL;
-}
-
-
 /* Decide which reflections can be scaled */
 static void select_scalable_reflections(struct image *images, int n)
 {
@@ -280,12 +142,9 @@ int main(int argc, char *argv[])
 	char *infile = NULL;
 	char *outfile = NULL;
 	char *geomfile = NULL;
-	char *prefix = NULL;
 	char *sym = NULL;
 	FILE *fh;
 	int nthreads = 1;
-	int config_basename = 0;
-	int config_checkprefix = 1;
 	struct detector *det;
 	unsigned int *cts;
 	ReflItemList *obs;
@@ -303,9 +162,6 @@ int main(int argc, char *argv[])
 		{"output",             1, NULL,               'o'},
 		{"geometry",           1, NULL,               'g'},
 		{"beam",               1, NULL,               'b'},
-		{"prefix",             1, NULL,               'x'},
-		{"basename",           0, &config_basename,    1},
-		{"no-check-prefix",    0, &config_checkprefix, 0},
 		{"symmetry",           1, NULL,               'y'},
 		{"iterations",         1, NULL,               'n'},
 		{0, 0, NULL, 0}
@@ -329,10 +185,6 @@ int main(int argc, char *argv[])
 			geomfile = strdup(optarg);
 			break;
 
-		case 'x' :
-			prefix = strdup(optarg);
-			break;
-
 		case 'j' :
 			nthreads = atoi(optarg);
 			break;
@@ -387,15 +239,6 @@ int main(int argc, char *argv[])
 		outfile = strdup("facetron.hkl");
 	}
 
-	/* Sanitise prefix */
-	if ( prefix == NULL ) {
-		prefix = strdup("");
-	} else {
-		if ( config_checkprefix ) {
-			prefix = check_prefix(prefix);
-		}
-	}
-
 	if ( sym == NULL ) sym = strdup("1");
 
 	/* Get detector geometry */
@@ -427,7 +270,11 @@ int main(int argc, char *argv[])
 
 		UnitCell *cell;
 		char *filename;
-		char *fnamereal;
+		char *rval;
+		char line[1024];
+		RefList *peaks;
+		RefList *transfer;
+		Reflection *refl;
 
 		if ( find_chunk(fh, &cell, &filename) == 1 ) {
 			ERROR("Couldn't get all of the filenames and cells"
@@ -437,38 +284,83 @@ int main(int argc, char *argv[])
 
 		images[i].indexed_cell = cell;
 
-		/* Mangle the filename now */
-		if ( config_basename ) {
-			char *tmp;
-			tmp = safe_basename(filename);
-			free(filename);
-			filename = tmp;
-		}
-		fnamereal = malloc(1024);
-		snprintf(fnamereal, 1023, "%s%s", prefix, filename);
-		free(filename);
-
-		images[i].filename = fnamereal;
+		images[i].filename = filename;
 		images[i].div = beam->divergence;
 		images[i].bw = beam->bandwidth;
 		images[i].det = det;
 		images[i].beam = beam;
 		images[i].osf = 1.0;
 		images[i].profile_radius = 0.005e9;
+		images[i].reflections = reflist_new();
 
 		/* Muppet proofing */
 		images[i].data = NULL;
 		images[i].flags = NULL;
 
-		/* Get reflections from this image.
-		 * FIXME: Use the ones from the stream */
-		integrate_image(&images[i], obs, sym);
+		/* Read integrated intensities from pattern */
+		peaks = reflist_new();
+		do {
+
+			Reflection *refl;
+			signed int h, k, l;
+			float intensity;
+			int r;
+
+			rval = fgets(line, 1023, fh);
+			chomp(line);
+
+			if ( (strlen(line) == 0) || (rval == NULL) ) break;
+
+			r = sscanf(line, "%i %i %i %f", &h, &k, &l, &intensity);
+			if ( r != 4 ) continue;
+
+			/* Not interested in the central beam */
+			if ( (h==0) && (k==0) && (l==0) ) continue;
+
+			refl = add_refl(peaks, h, k, l);
+			set_int(refl, intensity);
+
+		} while ( (strlen(line) != 0) && (rval != NULL) );
+
+		/* Calculate initial partialities and fill in intensities from
+		 * the stream */
+		transfer = find_intersections(&images[i], cell, 0);
+		images[i].reflections = reflist_new();
+
+		for ( refl = first_refl(transfer);
+		      refl != NULL;
+		      refl = next_refl(refl) ) {
+
+			Reflection *peak;
+			Reflection *new;
+			signed int h, k, l, ha, ka, la;
+			double r1, r2, p, x, y;
+			int clamp1, clamp2;
+
+			get_indices(refl, &h, &k, &l);
+			peak = find_refl(peaks, h, k, l);
+			if ( peak == NULL ) {
+				if ( (h==0) && (k==0) && (l==0) ) continue;
+				STATUS("%3i %3i %3i not found\n", h, k, l);
+				continue;
+			}
+
+			get_asymm(h, k, l, &ha, &ka, &la, sym);
+			add_item(obs, ha, ka, la);
+			new = add_refl(images[i].reflections, ha, ka, la);
+			get_partial(refl, &r1, &r2, &p, &clamp1, &clamp2);
+			get_detector_pos(refl, &x, &y);
+			set_partial(new, r1, r2, p, clamp1, clamp2);
+			set_detector_pos(new, 0.0, x, y);
+
+		}
+		reflist_free(peaks);
+		reflist_free(transfer);
 
 		progress_bar(i, n_total_patterns-1, "Loading pattern data");
 
 	}
 	fclose(fh);
-	free(prefix);
 
 	cts = new_list_count();
 
diff --git a/src/post-refinement.c b/src/post-refinement.c
index 25664e69..93b65810 100644
--- a/src/post-refinement.c
+++ b/src/post-refinement.c
@@ -28,6 +28,28 @@
 #include "cell.h"
 
 
+/* Maximum number of iterations of NLSq to do for each image per macrocycle. */
+#define MAX_CYCLES (100)
+
+
+/* Refineable parameters */
+enum {
+	REF_SCALE,
+	REF_DIV,
+	REF_R,
+	REF_ASX,
+	REF_BSX,
+	REF_CSX,
+	REF_ASY,
+	REF_BSY,
+	REF_CSY,
+	REF_ASZ,
+	REF_BSZ,
+	REF_CSZ,
+	NUM_PARAMS
+};
+
+
 /* Returns dp/dr at "r" */
 static double partiality_gradient(double r, double profile_radius)
 {
@@ -65,7 +87,7 @@ static double partiality_rgradient(double r, double profile_radius)
 
 
 /* Return the gradient of parameter 'k' given the current status of 'image'. */
-double gradient(struct image *image, int k, Reflection *refl, double r)
+static double gradient(struct image *image, int k, Reflection *refl, double r)
 {
 	double ds, tt, azi;
 	double nom, den;
@@ -187,7 +209,7 @@ static void apply_cell_shift(UnitCell *cell, int k, double shift)
 
 
 /* Apply the given shift to the 'k'th parameter of 'image'. */
-void apply_shift(struct image *image, int k, double shift)
+static void apply_shift(struct image *image, int k, double shift)
 {
 	switch ( k ) {
 
@@ -226,67 +248,19 @@ void apply_shift(struct image *image, int k, double shift)
 }
 
 
-double mean_partial_dev(struct image *image, RefList *reflections,
-                        const char *sym, double *i_full, FILE *graph)
-{
-	int n_used;
-	double delta_I = 0.0;
-	Reflection *refl;
-
-	n_used = 0;
-	for ( refl = first_refl(reflections);
-	      refl != NULL;
-	      refl = next_refl(refl) ) {
-
-		signed int hind, kind, lind;
-		signed int ha, ka, la;
-		double I_full;
-		float I_partial;
-		float xc, yc;
-		double x, y;
-		double p;
-
-		get_indices(refl, &hind, &kind, &lind);
-
-		if ( !get_scalable(refl) ) continue;
-
-		/* Actual measurement of this reflection from this pattern? */
-		get_detector_pos(refl, &x, &y);
-		/* FIXME: Get rid of this */
-		if ( integrate_peak(image, x, y,
-		                    &xc, &yc, &I_partial, NULL, NULL,
-		                    1, 0) ) {
-			continue;
-		}
-
-		get_asymm(hind, kind, lind, &ha, &ka, &la, sym);
-		I_full = lookup_intensity(i_full, ha, ka, la);
-		p = get_partiality(refl);
-		delta_I += fabs(I_partial - (p * I_full / image->osf));
-		n_used++;
-
-		if ( graph != NULL ) {
-			fprintf(graph, "%3i %3i %3i %5.2f (at %5.2f,%5.2f)"
-			               " %5.2f %5.2f\n",
-			       hind, kind, lind, I_partial/p, xc, yc,
-			       p, I_partial / I_full);
-		}
-
-	}
-
-	return delta_I / (double)n_used;
-}
-
-
 /* Perform one cycle of post refinement on 'image' against 'i_full' */
-double pr_iterate(struct image *image, double *i_full, const char *sym,
-                  RefList *reflections)
+static double pr_iterate(struct image *image, const double *i_full,
+                         const char *sym)
 {
 	gsl_matrix *M;
 	gsl_vector *v;
 	gsl_vector *shifts;
 	int param;
 	Reflection *refl;
+	RefList *reflections;
+	double max_shift;
+
+	reflections = image->reflections;
 
 	M = gsl_matrix_calloc(NUM_PARAMS, NUM_PARAMS);
 	v = gsl_vector_calloc(NUM_PARAMS);
@@ -355,14 +329,31 @@ double pr_iterate(struct image *image, double *i_full, const char *sym,
 
 	shifts = gsl_vector_alloc(NUM_PARAMS);
 	gsl_linalg_HH_solve(M, v, shifts);
+
+	max_shift = 0.0;
 	for ( param=0; param<NUM_PARAMS; param++ ) {
 		double shift = gsl_vector_get(shifts, param);
 		apply_shift(image, param, shift);
+		if ( fabs(shift) > max_shift ) max_shift = fabs(shift);
 	}
 
 	gsl_matrix_free(M);
 	gsl_vector_free(v);
 	gsl_vector_free(shifts);
 
-	return mean_partial_dev(image, reflections, sym, i_full, NULL);
+	return max_shift;
+}
+
+
+void pr_refine(struct image *image, const double *i_full, const char *sym)
+{
+	double max_shift;
+	int i;
+
+	i = 0;
+	do {
+		max_shift = pr_iterate(image, i_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 b3c6c074..a10fafe4 100644
--- a/src/post-refinement.h
+++ b/src/post-refinement.h
@@ -3,7 +3,7 @@
  *
  * Post refinement
  *
- * (c) 2006-2010 Thomas White <taw@physics.org>
+ * (c) 2006-2011 Thomas White <taw@physics.org>
  *
  * Part of CrystFEL - crystallography with a FEL
  *
@@ -21,36 +21,11 @@
 #include <stdio.h>
 
 #include "image.h"
-#include "reflist.h"
+#include "utils.h"
 
 
-/* Refineable parameters */
-enum {
-	REF_SCALE,
-	REF_DIV,
-	REF_R,
-	REF_ASX,
-	REF_BSX,
-	REF_CSX,
-	REF_ASY,
-	REF_BSY,
-	REF_CSY,
-	REF_ASZ,
-	REF_BSZ,
-	REF_CSZ,
-	NUM_PARAMS
-};
-
-/* Apply the given shift to the 'k'th parameter of 'image'. */
-void apply_shift(struct image *image, int k, double shift);
-
-
-double mean_partial_dev(struct image *image, RefList *reflections,
-                        const char *sym, double *i_full, FILE *graph);
-
-
-double pr_iterate(struct image *image, double *i_full, const char *sym,
-                  RefList *reflections);
+extern void pr_refine(struct image *image, const double *i_full,
+                      const char *sym);
 
 
 #endif	/* POST_REFINEMENT_H */
diff --git a/src/reflist.c b/src/reflist.c
index d297c4d0..2bb4e5f6 100644
--- a/src/reflist.c
+++ b/src/reflist.c
@@ -19,7 +19,6 @@
 struct _reflection {
 
 	/* Listy stuff */
-	RefList *list;
 	unsigned int serial;          /* Unique serial number, key */
 	struct _reflection *child[2]; /* Child nodes */
 	struct _reflection *parent;   /* Parent node */
@@ -224,32 +223,6 @@ void set_partial(Reflection *refl, double r1, double r2, double p,
 }
 
 
-void set_indices(Reflection *refl, signed int h, signed int k, signed int l)
-{
-	/* Tricky, because the indices determine the position in the tree */
-	Reflection copy;
-	Reflection *new;
-	Reflection transfer;
-
-	/* Copy all data */
-	copy = *refl;
-
-	/* Delete and re-add with new indices */
-	delete_refl(refl);
-	new = add_refl(copy.list, h, k, l);
-
-	/* Transfer data back */
-	transfer = *new;
-	*new = copy;
-	new->list = transfer.list;
-	new->parent = transfer.parent;
-	new->child[0] = transfer.child[0];
-	new->child[1] = transfer.child[1];
-	new->h = transfer.h;  new->k = transfer.k;  new->l = transfer.l;
-	new->serial = transfer.serial;
-}
-
-
 void set_int(Reflection *refl, double intensity)
 {
 	refl->intensity = intensity;
diff --git a/src/reflist.h b/src/reflist.h
index c7b2ad42..78aadf2a 100644
--- a/src/reflist.h
+++ b/src/reflist.h
@@ -46,8 +46,6 @@ 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,
                         double clamp_low, double clamp_high);
-extern void set_indices(Reflection *refl,
-                        signed int h, signed int k, signed int l);
 extern void set_int(Reflection *refl, double intensity);
 extern void set_scalable(Reflection *refl, int scalable);