Separate "refinable" and "scalable" concepts

7 files changed, 133 insertions, 56 deletions

diff --git a/doc/reference/CrystFEL-sections.txt b/doc/reference/CrystFEL-sections.txt
index 65b5eac0..af517540 100644
--- a/doc/reference/CrystFEL-sections.txt
+++ b/doc/reference/CrystFEL-sections.txt
@@ -23,6 +23,7 @@ get_symmetric_indices
 get_intensity
 get_partial
 get_scalable
+get_refinable
 get_redundancy
 get_sum_squared_dev
 get_esd_intensity
@@ -32,6 +33,7 @@ set_detector_pos
 set_partial
 set_int
 set_scalable
+set_refinable
 set_redundancy
 set_sum_squared_dev
 set_esd_intensity
diff --git a/src/geometry.c b/src/geometry.c
index e2454adf..7a11b110 100644
--- a/src/geometry.c
+++ b/src/geometry.c
@@ -275,11 +275,10 @@ RefList *find_intersections(struct image *image, UnitCell *cell)
 }
 
 
-/* Calculate partialities and apply them to the image's raw_reflections,
- * while adding to a ReflItemList of the currentl scalable (asymmetric)
- * reflections. */
-void update_partialities(struct image *image, const char *sym,
-                         int *n_expected, int *n_found, int *n_notfound)
+/* Predict reflections in "image" */
+void predict_corresponding_reflections(struct image *image, const char *sym,
+                                       int *n_expected, int *n_found,
+                                       int *n_notfound)
 {
 	Reflection *refl;
 	RefListIterator *iter;
@@ -346,3 +345,42 @@ void update_partialities(struct image *image, const char *sym,
 
 	reflist_free(predicted);
 }
+
+
+/* Calculate partialities and apply them to the image's raw_reflections */
+void update_partialities(struct image *image, const char *sym)
+{
+	Reflection *refl;
+	RefListIterator *iter;
+	RefList *predicted;
+
+	predicted = find_intersections(image, image->indexed_cell);
+
+	for ( refl = first_refl(image->reflections, &iter);
+	      refl != NULL;
+	      refl = next_refl(refl, iter) )
+	{
+
+		Reflection *p_peak;
+		double r1, r2, p;
+		signed int h, k, l;
+		int clamp1, clamp2;
+
+		/* Get predicted indices and location */
+		get_symmetric_indices(refl, &h, &k, &l);
+
+		/* Look for this reflection in the pattern */
+		p_peak = find_refl(predicted, h, k, l);
+		if ( p_peak == NULL ) {
+			set_partial(refl, 0.0, 0.0, 0.0, -1, +1);
+			continue;
+		} else {
+			/* Transfer partiality stuff */
+			get_partial(p_peak, &r1, &r2, &p, &clamp1, &clamp2);
+			set_partial(refl, r1, r2, p, clamp1, clamp2);
+		}
+
+	}
+
+	reflist_free(predicted);
+}
diff --git a/src/geometry.h b/src/geometry.h
index 892ff747..efda6c87 100644
--- a/src/geometry.h
+++ b/src/geometry.h
@@ -21,7 +21,11 @@
 
 extern RefList *find_intersections(struct image *image, UnitCell *cell);
 
-extern void update_partialities(struct image *image, const char *sym,
-                                int *n_expected, int *n_found, int *n_notfound);
+extern void predict_corresponding_reflections(struct image *image,
+                                              const char *sym, int *n_expected,
+                                              int *n_found, int *n_notfound);
+
+extern void update_partialities(struct image *image, const char *sym);
+
 
 #endif	/* GEOMETRY_H */
diff --git a/src/partialator.c b/src/partialator.c
index 167f736d..43c4bda0 100644
--- a/src/partialator.c
+++ b/src/partialator.c
@@ -191,6 +191,35 @@ static int select_scalable_reflections(RefList *list, ReflItemList *sc_l)
 }
 
 
+static void select_reflections_for_refinement(struct image *images, int n,
+                                              ReflItemList *scalable)
+{
+	int i;
+
+	for ( i=0; i<n; i++ ) {
+
+		Reflection *refl;
+		RefListIterator *iter;
+
+		for ( refl = first_refl(images[i].reflections, &iter);
+		      refl != NULL;
+		      refl = next_refl(refl, iter) )
+		{
+			signed int h, k, l;
+			int sc;
+
+			get_indices(refl, &h, &k, &l);
+			sc = get_scalable(refl);
+
+			if ( sc && find_item(scalable, h, k, l) ) {
+				set_refinable(refl, 1);
+			}
+		}
+
+	}
+}
+
+
 int main(int argc, char *argv[])
 {
 	int c;
@@ -397,8 +426,8 @@ int main(int argc, char *argv[])
 		reflist_free(cur->reflections);
 		cur->reflections = as;
 
-		update_partialities(cur, sym,
-		                    &n_expected, &n_found, &n_notfound);
+		predict_corresponding_reflections(cur, sym, &n_expected,
+		                                 &n_found, &n_notfound);
 
 		nobs += select_scalable_reflections(cur->reflections, scalable);
 
@@ -419,6 +448,8 @@ int main(int argc, char *argv[])
 	full = scale_intensities(images, n_usable_patterns, sym,
 	                         scalable, cref, reference);
 
+	select_reflections_for_refinement(images, n_usable_patterns, scalable);
+
 	for ( i=0; i<num_items(scalable); i++ ) {
 		Reflection *f;
 		struct refl_item *it = get_item(scalable, i);
@@ -459,6 +490,7 @@ int main(int argc, char *argv[])
 		FILE *fhg;
 		FILE *fhp;
 		char filename[1024];
+		int j;
 
 		STATUS("Post refinement cycle %i of %i\n", i+1, n_iter);
 
@@ -484,9 +516,14 @@ int main(int argc, char *argv[])
 
 		nobs = 0;
 		clear_items(scalable);
-		for ( i=0; i<n_usable_patterns; i++ ) {
+		for ( j=0; j<n_usable_patterns; j++ ) {
+
+			struct image *cur = &images[j];
+
+			predict_corresponding_reflections(cur, sym, &n_expected,
+		                                          &n_found,
+		                                          &n_notfound);
 
-			struct image *cur = &images[i];
 			nobs += select_scalable_reflections(cur->reflections,
 			                                    scalable);
 
diff --git a/src/post-refinement.c b/src/post-refinement.c
index 7eef5dfa..aeafcddc 100644
--- a/src/post-refinement.c
+++ b/src/post-refinement.c
@@ -31,7 +31,7 @@
 
 
 /* Maximum number of iterations of NLSq to do for each image per macrocycle. */
-#define MAX_CYCLES (5)
+#define MAX_CYCLES (10)
 
 
 /* Returns dp/dr at "r" */
@@ -347,7 +347,7 @@ static double pr_iterate(struct image *image, const RefList *full,
 		Reflection *match;
 		double gradients[NUM_PARAMS];
 
-		if ( !get_scalable(refl) ) continue;
+		if ( !get_refinable(refl) ) continue;
 
 		/* Find the full version */
 		get_indices(refl, &ha, &ka, &la);
@@ -484,66 +484,33 @@ void pr_refine(struct image *image, const RefList *full, const char *sym)
 	double max_shift, dev;
 	int i;
 	const int verbose = 1;
-	int nexp, nfound, nnotfound;
-
-	update_partialities(image, sym, &nexp, &nfound, &nnotfound);
 
 	if ( verbose ) {
 		dev = mean_partial_dev(image, full, sym);
-		STATUS("PR starting dev = %5.2f (%i out of %i found)\n",
-		       dev, nfound, nexp);
-	}
-
-	if ( (double)nfound/(double)nexp < 0.5 ) {
-		ERROR("Refusing to refine this image: %i out of %i found\n",
-		       nfound, nexp);
-		return;
+		STATUS("PR starting dev = %5.2f\n", dev);
 	}
 
 	i = 0;
+	image->pr_dud = 0;
 	do {
 
 		double asx, asy, asz;
 		double bsx, bsy, bsz;
 		double csx, csy, csz;
 		double dev;
-		int old_nexp, old_nfound;
 
 		cell_get_reciprocal(image->indexed_cell, &asx, &asy, &asz,
 			               &bsx, &bsy, &bsz, &csx, &csy, &csz);
-		old_nexp = nexp;
-		old_nfound = nfound;
 
 		max_shift = pr_iterate(image, full, sym);
 
-		update_partialities(image, sym, &nexp, &nfound, &nnotfound);
+		update_partialities(image, sym);
 
 		if ( verbose ) {
 			dev = mean_partial_dev(image, full, sym);
-			STATUS("PR Iteration %2i: max shift = %5.2f"
-			       " dev = %5.2f (%i out of %i found)\n",
-			       i+1, max_shift, dev, nfound, nexp);
-		}
-
-		if ( (double)nfound / (double)nexp < 0.5 ) {
-
-			if ( verbose ) {
-				ERROR("Bad refinement step - backtracking.\n");
-				ERROR("I'll come back to this image later.\n");
-			}
-
-			cell_set_reciprocal(image->indexed_cell, asx, asy, asz,
-				                  bsx, bsy, bsz, csx, csy, csz);
-
-			update_partialities(image, sym,
-		                            &nexp, &nfound, &nnotfound);
-
-			image->pr_dud = 1;
-
-			return;
-
-		} else {
-			image->pr_dud = 0;
+			STATUS("PR Iteration %2i: max shift = %10.2f"
+			       " dev = %10.5e\n",
+			       i+1, max_shift, dev);
 		}
 
 		i++;
diff --git a/src/reflist.c b/src/reflist.c
index a8bfb4f6..19d5b1a4 100644
--- a/src/reflist.c
+++ b/src/reflist.c
@@ -67,6 +67,10 @@ struct _refldata {
 	/* Non-zero if this reflection can be used for scaling */
 	int scalable;
 
+	/* Non-zero if this reflection should be used as a "guide star" for
+	 * post refinement */
+	int refinable;
+
 	/* Intensity */
 	double intensity;
 	double esd_i;
@@ -371,8 +375,7 @@ void get_partial(const Reflection *refl, double *r1, double *r2, double *p,
  * get_scalable:
  * @refl: A %Reflection
  *
- * Returns: non-zero if this reflection was marked as useful for scaling and
- * post refinement.
+ * Returns: non-zero if this reflection can be scaled.
  *
  **/
 int get_scalable(const Reflection *refl)
@@ -382,6 +385,19 @@ int get_scalable(const Reflection *refl)
 
 
 /**
+ * get_refinable:
+ * @refl: A %Reflection
+ *
+ * Returns: non-zero if this reflection can be used for post refinement.
+ *
+ **/
+int get_refinable(const Reflection *refl)
+{
+	return refl->data.refinable;
+}
+
+
+/**
  * get_redundancy:
  * @refl: A %Reflection
  *
@@ -523,8 +539,7 @@ void set_int(Reflection *refl, double intensity)
 /**
  * set_scalable:
  * @refl: A %Reflection
- * @scalable: Non-zero if this reflection was marked as useful for scaling and
- * post refinement.
+ * @scalable: Non-zero if this reflection should be scaled.
  *
  **/
 void set_scalable(Reflection *refl, int scalable)
@@ -534,6 +549,18 @@ void set_scalable(Reflection *refl, int scalable)
 
 
 /**
+ * set_refinable:
+ * @refl: A %Reflection
+ * @scalable: Non-zero if this reflection can be used for post refinement.
+ *
+ **/
+void set_refinable(Reflection *refl, int refinable)
+{
+	refl->data.refinable = refinable;
+}
+
+
+/**
  * set_redundancy:
  * @refl: A %Reflection
  * @red: New redundancy for the reflection
diff --git a/src/reflist.h b/src/reflist.h
index ffdef55a..0d836edf 100644
--- a/src/reflist.h
+++ b/src/reflist.h
@@ -68,6 +68,7 @@ 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(const Reflection *refl);
+extern int get_refinable(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);
@@ -81,6 +82,7 @@ extern void set_partial(Reflection *refl, double r1, double r2, double p,
                         double clamp_low, double clamp_high);
 extern void set_int(Reflection *refl, double intensity);
 extern void set_scalable(Reflection *refl, int scalable);
+extern void set_refinable(Reflection *refl, int refinable);
 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);