Make panel assignments invariant during prediction- and post-refinement

Reflections appearing and disappearing are problematic when trying to do a
least-squares refinement.  Therefore, assume that reflections stay on
panel and keep them under consideration even if their partialities go to
zero (i.e. they drift off Bragg).  This should stabilise both
refinements, and simplifies quite a lot of code.

Collateral "damage": the old "select_intersection()" is now gone.

12 files changed, 156 insertions, 256 deletions

diff --git a/libcrystfel/src/geometry.c b/libcrystfel/src/geometry.c
index b55b696d..35492fde 100644
--- a/libcrystfel/src/geometry.c
+++ b/libcrystfel/src/geometry.c
@@ -47,58 +47,66 @@
 #include "geometry.h"
 
 
-static signed int locate_peak(double x, double y, double z, double k,
-                              struct detector *det, double *xdap, double *ydap)
+static int locate_peak_on_panel(double x, double y, double z, double k,
+                                struct panel *p,
+                                double *pfs, double *pss)
 {
-	int i;
-	signed int found = -1;
 	const double den = k + z;
+	double fs, ss, plx, ply, xd, yd;
 
-	*xdap = -1;  *ydap = -1;
+	/* Coordinates of peak relative to central beam, in m */
+	xd = p->clen * x / den;
+	yd = p->clen * y / den;
 
-	for ( i=0; i<det->n_panels; i++ ) {
+	/* Convert to pixels */
+	xd *= p->res;
+	yd *= p->res;
 
-		double xd, yd;
-		double fs, ss, plx, ply;
-		struct panel *p;
+	/* Convert to relative to the panel corner */
+	plx = xd - p->cnx;
+	ply = yd - p->cny;
+
+	fs = p->xfs*plx + p->yfs*ply;
+	ss = p->xss*plx + p->yss*ply;
+
+	fs += p->min_fs;
+	ss += p->min_ss;
+
+	*pfs = fs;  *pss = ss;
+
+	/* Now, is this on this panel? */
+	if ( fs < p->min_fs ) return 0;
+	if ( fs > p->max_fs ) return 0;
+	if ( ss < p->min_ss ) return 0;
+	if ( ss > p->max_ss ) return 0;
 
-		p = &det->panels[i];
 
-		/* Coordinates of peak relative to central beam, in m */
-		xd = p->clen * x / den;
-		yd = p->clen * y / den;
+	return 1;
+}
 
-		/* Convert to pixels */
-		xd *= p->res;
-		yd *= p->res;
+static signed int locate_peak(double x, double y, double z, double k,
+                              struct detector *det, double *pfs, double *pss)
+{
+	int i;
 
-		/* Convert to relative to the panel corner */
-		plx = xd - p->cnx;
-		ply = yd - p->cny;
+	*pfs = -1;  *pss = -1;
 
-		fs = p->xfs*plx + p->yfs*ply;
-		ss = p->xss*plx + p->yss*ply;
+	for ( i=0; i<det->n_panels; i++ ) {
 
-		fs += p->min_fs;
-		ss += p->min_ss;
+		struct panel *p;
+
+		p = &det->panels[i];
 
-		/* Now, is this on this panel? */
-		if ( fs < p->min_fs ) continue;
-		if ( fs > p->max_fs ) continue;
-		if ( ss < p->min_ss ) continue;
-		if ( ss > p->max_ss ) continue;
+		if ( locate_peak_on_panel(x, y, z, k, p, pfs, pss) ) {
 
-		/* If peak appears on multiple panels, reject it */
-		if ( found != -1 ) return -1;
+			/* Woohoo! */
+			return i;
 
-		/* Woohoo! */
-		found = i;
-		*xdap = fs;
-		*ydap = ss;
+		}
 
 	}
 
-	return found;
+	return -1;
 }
 
 
@@ -182,7 +190,8 @@ static double partiality(PartialityModel pmodel,
 static Reflection *check_reflection(struct image *image, Crystal *cryst,
                                     PartialityModel pmodel,
                                     signed int h, signed int k, signed int l,
-                                    double xl, double yl, double zl)
+                                    double xl, double yl, double zl,
+                                    Reflection *updateme)
 {
 	const int output = 0;
 	double tl;
@@ -195,7 +204,7 @@ static Reflection *check_reflection(struct image *image, Crystal *cryst,
 	double del;
 
 	/* Don't predict 000 */
-	if ( abs(h)+abs(k)+abs(l) == 0 ) return NULL;
+	if ( (updateme == NULL) && (abs(h)+abs(k)+abs(l) == 0) ) return NULL;
 
 	pr = crystal_get_profile_radius(cryst);
 	del = image->div + crystal_get_mosaicity(cryst);
@@ -207,7 +216,7 @@ static Reflection *check_reflection(struct image *image, Crystal *cryst,
 	khigh = 1.0/(image->lambda + image->lambda*image->bw/2.0);
 
 	/* If the point is looking "backscattery", reject it straight away */
-	if ( zl < -khigh/2.0 ) return NULL;
+	if ( (updateme == NULL) && (zl < -khigh/2.0) ) return NULL;
 
 	tl = sqrt(xl*xl + yl*yl);
 
@@ -220,28 +229,47 @@ static Reflection *check_reflection(struct image *image, Crystal *cryst,
 	rlow = klow - distance(cet, cez, tl, zl);  /* Loss of precision */
 
 	/* Condition for reflection to be excited at all */
-	if ( (signbit(rlow) == signbit(rhigh))
+	if ( (updateme == NULL)
+	     && (signbit(rlow) == signbit(rhigh))
 	     && (fabs(rlow) > pr)
 	     && (fabs(rhigh) > pr) ) return NULL;
 
 	/* Calculate partiality */
 	part = partiality(pmodel, rlow, rhigh, pr);
 
-	/* Add peak to list */
-	refl = reflection_new(h, k, l);
+	if ( updateme == NULL ) {
+		refl = reflection_new(h, k, l);
+	} else {
+		refl = updateme;
+	}
+
+	/* If we are updating a previous reflection, assume it stays
+	 * on the same panel and calculate the new position even if it's
+	 * fallen off the edge of the panel. */
+	if ( updateme != NULL ) {
+
+		double fs, ss;
+		locate_peak_on_panel(xl, yl, zl, 1.0/image->lambda,
+		                     get_panel(updateme), &fs, &ss);
+		set_detector_pos(refl, fs, ss);
+
+	}
 
-	/* If we have detector information, check the spot is measured.
-	 * Otherwise, we make do with calculating the partialiaty etc. */
-	if ( image->det != NULL ) {
-		double xda, yda;        /* Position on detector */
-		signed int p;           /* Panel number */
-		p = locate_peak(xl, yl, zl, 1.0/image->lambda, image->det,
-		                &xda, &yda);
+	/* Otherwise, calculate position if we have a detector structure, and
+	 * if we don't then just make do with partiality calculation */
+	if ( (image->det != NULL) && (updateme == NULL) ) {
+
+		double fs, ss;        /* Position on detector */
+		signed int p;         /* Panel number */
+		p = locate_peak(xl, yl, zl, 1.0/image->lambda,
+		                image->det, &fs, &ss);
 		if ( p == -1 ) {
 			reflection_free(refl);
 			return NULL;
 		}
-		set_detector_pos(refl, 0.0, xda, yda);
+		set_detector_pos(refl, fs, ss);
+		set_panel(refl, &image->det->panels[p]);
+
 	}
 
 	if ( unlikely(rlow < rhigh) ) {
@@ -249,8 +277,11 @@ static Reflection *check_reflection(struct image *image, Crystal *cryst,
 		ERROR("%3i %3i %3i  rlow = %e, rhigh = %e\n",
 		      h, k, l, rlow, rhigh);
 		ERROR("div + m = %e, R = %e, bw = %e\n", del, pr, image->bw);
-		reflection_free(refl);
-		return NULL;
+		/* If we are updating, this is (kind of) OK */
+		if ( updateme == NULL ) {
+			reflection_free(refl);
+			return NULL;
+		}
 	}
 
 	set_partial(refl, rlow, rhigh, part);
@@ -339,7 +370,7 @@ RefList *find_intersections_to_res(struct image *image, Crystal *cryst,
 		zl = h*asz + k*bsz + l*csz;
 
 		refl = check_reflection(image, cryst, pmodel,
-		                        h, k, l, xl, yl, zl);
+		                        h, k, l, xl, yl, zl, NULL);
 
 		if ( refl != NULL ) {
 			add_refl_to_list(refl, reflections);
@@ -353,67 +384,6 @@ RefList *find_intersections_to_res(struct image *image, Crystal *cryst,
 }
 
 
-/* Deprecated: select reflections using Kirian-style pixel proximity */
-RefList *select_intersections(struct image *image, Crystal *cryst)
-{
-	double ax, ay, az;
-	double bx, by, bz;
-	double cx, cy, cz;
-	const double min_dist = 0.25;
-	RefList *list;
-	int i;
-
-	/* Round towards nearest */
-	fesetround(1);
-
-	/* Cell basis vectors for this image */
-	cell_get_cartesian(crystal_get_cell(cryst), &ax, &ay, &az,
-	                   &bx, &by, &bz, &cx, &cy, &cz);
-
-	list = reflist_new();
-	if ( list == NULL ) return NULL;
-
-	/* Loop over peaks, checking proximity to nearest reflection */
-	for ( i=0; i<image_feature_count(image->features); i++ ) {
-
-		struct imagefeature *f;
-		struct rvec q;
-		double h, k, l, hd, kd, ld;
-		double dsq;
-
-		f = image_get_feature(image->features, i);
-		if ( f == NULL ) continue;
-
-		/* Reciprocal space position of found peak */
-		q = get_q(image, f->fs, f->ss, NULL, 1.0/image->lambda);
-
-		/* Decimal and fractional Miller indices of nearest
-		 * reciprocal lattice point */
-		hd = q.u * ax + q.v * ay + q.w * az;
-		kd = q.u * bx + q.v * by + q.w * bz;
-		ld = q.u * cx + q.v * cy + q.w * cz;
-		h = lrint(hd);
-		k = lrint(kd);
-		l = lrint(ld);
-
-		/* Check distance */
-		dsq = pow(h-hd, 2.0) + pow(k-kd, 2.0) + pow(l-ld, 2.0);
-
-		if ( sqrt(dsq) < min_dist ) {
-
-			Reflection *refl;
-
-			refl = add_refl(list, h, k, l);
-			set_detector_pos(refl, sqrt(dsq), f->fs, f->ss);
-
-		}
-
-	}
-
-	return list;
-}
-
-
 static void set_unity_partialities(Crystal *cryst)
 {
 	Reflection *refl;
@@ -430,8 +400,7 @@ static void set_unity_partialities(Crystal *cryst)
 
 
 /* Calculate partialities and apply them to the image's reflections */
-void update_partialities_2(Crystal *cryst, PartialityModel pmodel,
-                           int *n_gained, int *n_lost, double *mean_p_change)
+void update_partialities(Crystal *cryst, PartialityModel pmodel)
 {
 	Reflection *refl;
 	RefListIterator *iter;
@@ -439,8 +408,6 @@ void update_partialities_2(Crystal *cryst, PartialityModel pmodel,
 	double bsx, bsy, bsz;
 	double csx, csy, csz;
 	struct image *image = crystal_get_image(cryst);
-	double total_p_change = 0.0;
-	int n = 0;
 
 	if ( pmodel == PMODEL_UNITY ) {
 		set_unity_partialities(cryst);
@@ -454,69 +421,20 @@ void update_partialities_2(Crystal *cryst, PartialityModel pmodel,
 	      refl != NULL;
 	      refl = next_refl(refl, iter) )
 	{
-		Reflection *vals;
-		double r1, r2, L, p, x, y;
 		double xl, yl, zl;
 		signed int h, k, l;
-		double old_p;
 
 		get_symmetric_indices(refl, &h, &k, &l);
-		old_p = get_partiality(refl);
 
 		/* Get the coordinates of the reciprocal lattice point */
 		xl = h*asx + k*bsx + l*csx;
 		yl = h*asy + k*bsy + l*csy;
 		zl = h*asz + k*bsz + l*csz;
 
-		vals = check_reflection(image, cryst, pmodel,
-		                        h, k, l, xl, yl, zl);
-
-		if ( vals == NULL ) {
-
-			if ( get_redundancy(refl) != 0 ) {
-				(*n_lost)++;
-				set_partiality(refl, 0.0);
-				set_redundancy(refl, 0);
-			}
-
-		} else {
-
-			if ( get_redundancy(refl) == 0 ) {
-				(*n_gained)++;
-				set_redundancy(refl, 1);
-			}
-
-			/* Transfer partiality stuff */
-			get_partial(vals, &r1, &r2, &p);
-			set_partial(refl, r1, r2, p);
-			L = get_lorentz(vals);
-			set_lorentz(refl, L);
-
-			/* Transfer detector location */
-			get_detector_pos(vals, &x, &y);
-			set_detector_pos(refl, 0.0, x, y);
-
-			reflection_free(vals);
-
-			total_p_change += fabs(p - old_p);
-			n++;
-
-		}
+		check_reflection(image, cryst, pmodel,
+		                 h, k, l, xl, yl, zl, refl);
 
 	}
-
-	*mean_p_change = total_p_change / n;
-}
-
-
-/* Wrapper to maintain API compatibility */
-void update_partialities(Crystal *cryst, PartialityModel pmodel)
-{
-	int n_gained = 0;
-	int n_lost = 0;
-	double mean_p_change = 0.0;
-	update_partialities_2(cryst, pmodel, &n_gained, &n_lost,
-	                      &mean_p_change);
 }
 
 
diff --git a/libcrystfel/src/geometry.h b/libcrystfel/src/geometry.h
index f39b6f8b..162c0355 100644
--- a/libcrystfel/src/geometry.h
+++ b/libcrystfel/src/geometry.h
@@ -3,12 +3,12 @@
  *
  * Geometry of diffraction
  *
- * Copyright © 2013-2014 Deutsches Elektronen-Synchrotron DESY,
+ * Copyright © 2013-2015 Deutsches Elektronen-Synchrotron DESY,
  *                       a research centre of the Helmholtz Association.
  * Copyright © 2012 Richard Kirian
  *
  * Authors:
- *   2010-2014 Thomas White <taw@physics.org>
+ *   2010-2015 Thomas White <taw@physics.org>
  *   2012      Richard Kirian
  *
  * This file is part of CrystFEL.
@@ -67,13 +67,7 @@ extern RefList *find_intersections_to_res(struct image *image, Crystal *cryst,
                                           PartialityModel pmodel,
 					  double max_res);
 
-/* Deprecated: select reflections using Kirian-style pixel proximity */
-extern RefList *select_intersections(struct image *image, Crystal *cryst);
-
 extern void update_partialities(Crystal *cryst, PartialityModel pmodel);
-extern void update_partialities_2(Crystal *cryst, PartialityModel pmodel,
-                                  int *n_gained, int *n_lost,
-                                  double *mean_p_change);
 extern void polarisation_correction(RefList *list, UnitCell *cell,
                                     struct image *image);
 
diff --git a/libcrystfel/src/integration.c b/libcrystfel/src/integration.c
index a05f64f1..0e35037d 100644
--- a/libcrystfel/src/integration.c
+++ b/libcrystfel/src/integration.c
@@ -1193,7 +1193,7 @@ static void integrate_prof2d_once(struct intcontext *ic, struct peak_box *bx,
 		get_detector_pos(bx->refl, &pfs, &pss);
 		pfs += bx->offs_fs;
 		pss += bx->offs_ss;
-		set_detector_pos(bx->refl, 0.0, pfs, pss);
+		set_detector_pos(bx->refl, pfs, pss);
 
 		if ( bx->intensity < -5.0*bx->sigma ) {
 			ic->n_implausible++;
@@ -1447,7 +1447,7 @@ static void integrate_rings_once(Reflection *refl, struct image *image,
 	/* Update position */
 	pfs += bx->offs_fs;
 	pss += bx->offs_ss;
-	set_detector_pos(refl, 0.0, pfs, pss);
+	set_detector_pos(refl, pfs, pss);
 
 	if ( get_int_diag(ic, refl) ) show_peak_box(ic, bx, results_pipe);
 
diff --git a/libcrystfel/src/reflist-utils.c b/libcrystfel/src/reflist-utils.c
index 722f800c..82943158 100644
--- a/libcrystfel/src/reflist-utils.c
+++ b/libcrystfel/src/reflist-utils.c
@@ -390,7 +390,7 @@ RefList *read_reflections_from_file(FILE *fh)
 
 			refl = add_refl(out, h, k, l);
 			set_intensity(refl, intensity);
-			set_detector_pos(refl, 0.0, fs, ss);
+			set_detector_pos(refl, fs, ss);
 
 			set_esd_intensity(refl, sigma);
 			set_redundancy(refl, cts);
diff --git a/libcrystfel/src/reflist.c b/libcrystfel/src/reflist.c
index 33c1f948..795894a1 100644
--- a/libcrystfel/src/reflist.c
+++ b/libcrystfel/src/reflist.c
@@ -75,10 +75,7 @@ struct _refldata {
 	/* Location in image */
 	double fs;
 	double ss;
-
-	/* The distance from the exact Bragg position to the coordinates
-	 * given above. */
-	double excitation_error;
+	struct panel *panel;
 
 	/* Non-zero if this reflection can be used for scaling */
 	int scalable;
@@ -316,29 +313,31 @@ Reflection *next_found_refl(Reflection *refl)
 
 /********************************** Getters ***********************************/
 
+
 /**
- * get_excitation_error:
+ * get_detector_pos:
  * @refl: A %Reflection
+ * @fs: Location at which to store the fast scan offset of the reflection
+ * @ss: Location at which to store the slow scan offset of the reflection
  *
- * Returns: The excitation error for the reflection.
  **/
-double get_excitation_error(const Reflection *refl)
+void get_detector_pos(const Reflection *refl, double *fs, double *ss)
 {
-	return refl->data.excitation_error;
+	*fs = refl->data.fs;
+	*ss = refl->data.ss;
 }
 
 
 /**
- * get_detector_pos:
+ * get_panel:
  * @refl: A %Reflection
- * @fs: Location at which to store the fast scan offset of the reflection
- * @ss: Location at which to store the slow scan offset of the reflection
+ *
+ * Returns: the panel which the reflection appears on
  *
  **/
-void get_detector_pos(const Reflection *refl, double *fs, double *ss)
+struct panel *get_panel(const Reflection *refl)
 {
-	*fs = refl->data.fs;
-	*ss = refl->data.ss;
+	return refl->data.panel;
 }
 
 
@@ -570,20 +569,32 @@ void copy_data(Reflection *to, const Reflection *from)
 /**
  * set_detector_pos:
  * @refl: A %Reflection
- * @exerr: The excitation error for this reflection
  * @fs: The fast scan offset of the reflection
  * @ss: The slow scan offset of the reflection
  *
  **/
-void set_detector_pos(Reflection *refl, double exerr, double fs, double ss)
+void set_detector_pos(Reflection *refl, double fs, double ss)
 {
-	refl->data.excitation_error = exerr;
 	refl->data.fs = fs;
 	refl->data.ss = ss;
 }
 
 
 /**
+ * set_panel:
+ * @refl: A %Reflection
+ * @panel: Pointer to the panel structure on which the reflection appears
+ *
+ * Note that the pointer will be stored, not the contents of the structure.
+ *
+ **/
+void set_panel(Reflection *refl, struct panel *p)
+{
+	refl->data.panel = p;
+}
+
+
+/**
  * set_partial:
  * @refl: A %Reflection
  * @rlow: The "low" excitation error
diff --git a/libcrystfel/src/reflist.h b/libcrystfel/src/reflist.h
index 85a87c54..8b0d8e04 100644
--- a/libcrystfel/src/reflist.h
+++ b/libcrystfel/src/reflist.h
@@ -3,11 +3,11 @@
  *
  * Fast reflection/peak list
  *
- * Copyright © 2012-2014 Deutsches Elektronen-Synchrotron DESY,
+ * Copyright © 2012-2015 Deutsches Elektronen-Synchrotron DESY,
  *                       a research centre of the Helmholtz Association.
  *
  * Authors:
- *   2011-2014 Thomas White <taw@physics.org>
+ *   2011-2015 Thomas White <taw@physics.org>
  *
  * This file is part of CrystFEL.
  *
@@ -84,8 +84,8 @@ extern Reflection *find_refl(const RefList *list, signed int h, signed int k, si
 extern Reflection *next_found_refl(Reflection *refl);
 
 /* Get */
-extern double get_excitation_error(const Reflection *refl);
 extern void get_detector_pos(const Reflection *refl, double *fs, double *ss);
+extern struct panel *get_panel(const Reflection *refl);
 extern double get_partiality(const Reflection *refl);
 extern double get_lorentz(const Reflection *refl);
 extern void get_indices(const Reflection *refl,
@@ -106,8 +106,8 @@ extern double get_mean_bg(const Reflection *refl);
 
 /* Set */
 extern void copy_data(Reflection *to, const Reflection *from);
-extern void set_detector_pos(Reflection *refl, double exerr,
-                             double fs, double ss);
+extern void set_detector_pos(Reflection *refl, double fs, double ss);
+extern void set_panel(Reflection *refl, struct panel *p);
 extern void set_partial(Reflection *refl, double rlow, double rhigh, double p);
 extern void set_partiality(Reflection *refl, double p);
 extern void set_lorentz(Reflection *refl, double L);
diff --git a/libcrystfel/src/stream.c b/libcrystfel/src/stream.c
index 8586c850..78ac9a40 100644
--- a/libcrystfel/src/stream.c
+++ b/libcrystfel/src/stream.c
@@ -319,7 +319,7 @@ static RefList *read_stream_reflections_2_3(FILE *fh, struct detector *det)
 				p = find_panel_by_name(det,pn);
 				write_fs = fs - p->orig_min_fs + p->min_fs;
 				write_ss = ss - p->orig_min_ss + p->min_ss;
-				set_detector_pos(refl, 0.0, write_fs, write_ss);
+				set_detector_pos(refl, write_fs, write_ss);
 			}
 			set_esd_intensity(refl, sigma);
 			set_peak(refl, pk);
@@ -383,11 +383,11 @@ static RefList *read_stream_reflections_2_1(FILE *fh, struct detector *det)
 				p = find_orig_panel(det, fs, ss);
 				write_fs = fs - p->orig_min_fs + p->min_fs;
 				write_ss = ss - p->orig_min_ss + p->min_ss;
-				set_detector_pos(refl, 0.0, write_fs, write_ss);
+				set_detector_pos(refl, write_fs, write_ss);
 
 			} else {
 
-				set_detector_pos(refl, 0.0, fs, ss);
+				set_detector_pos(refl, fs, ss);
 
 			}
 			set_esd_intensity(refl, sigma);
@@ -448,11 +448,11 @@ static RefList *read_stream_reflections_2_2(FILE *fh, struct detector *det)
 				p = find_orig_panel(det, fs, ss);
 				write_fs = fs - p->orig_min_fs + p->min_fs;
 				write_ss = ss - p->orig_min_ss + p->min_ss;
-				set_detector_pos(refl, 0.0, write_fs, write_ss);
+				set_detector_pos(refl, write_fs, write_ss);
 
 			} else {
 
-				set_detector_pos(refl, 0.0, fs, ss);
+				set_detector_pos(refl, fs, ss);
 
 			}
 
diff --git a/src/partialator.c b/src/partialator.c
index 4b1d5b18..ca3a9976 100644
--- a/src/partialator.c
+++ b/src/partialator.c
@@ -626,17 +626,12 @@ int main(int argc, char *argv[])
 		for ( j=0; j<images[i].n_crystals; j++ ) {
 
 			Crystal *cryst;
-			int n_gained = 0;
-			int n_lost = 0;
-			double mean_p_change = 0.0;
 
 			cryst = images[i].crystals[j];
 			crystal_set_image(cryst, &images[i]);
 
 			/* Now it's safe to do the following */
-			update_partialities_2(cryst, pmodel, &n_gained, &n_lost,
-			                      &mean_p_change);
-			assert(n_gained == 0);  /* That'd just be silly */
+			update_partialities(cryst, pmodel);
 
 		}
 	}
diff --git a/src/post-refinement.c b/src/post-refinement.c
index 96c0de12..745e5bd7 100644
--- a/src/post-refinement.c
+++ b/src/post-refinement.c
@@ -614,33 +614,17 @@ struct prdata pr_refine(Crystal *cr, const RefList *full,
 		double bsx, bsy, bsz;
 		double csx, csy, csz;
 		double dev;
-		int n_total;
-		int n_gained = 0;
-		int n_lost = 0;
 
-		n_total = num_reflections(crystal_get_reflections(cr));
 		cell_get_reciprocal(crystal_get_cell(cr), &asx, &asy, &asz,
 			               &bsx, &bsy, &bsz, &csx, &csy, &csz);
 
 		pr_iterate(cr, full, pmodel, &prdata.n_filtered);
 
-		update_partialities_2(cr, pmodel, &n_gained, &n_lost,
-		                      &mean_p_change);
+		update_partialities(cr, pmodel);
 
 		if ( verbose ) {
 			dev = guide_dev(cr, full);
-			STATUS("PR Iteration %2i: mean p change = %10.2f"
-			       " dev = %10.5e, %i gained, %i lost, %i total\n",
-			       i+1, mean_p_change, dev,
-			       n_gained, n_lost, n_total);
-		}
-
-		if ( 3*n_lost > n_total ) {
-			revert_crystal(cr, backup);
-			update_partialities_2(cr, pmodel, &n_gained, &n_lost,
-			                      &mean_p_change);
-			crystal_set_user_flag(cr, 4);
-			break;
+			STATUS("PR Iteration %2i: dev = %10.5e\n", i+1, dev);
 		}
 
 		i++;
diff --git a/src/predict-refine.c b/src/predict-refine.c
index 0e00dc5c..56a4cc2c 100644
--- a/src/predict-refine.c
+++ b/src/predict-refine.c
@@ -51,10 +51,12 @@ struct reflpeak {
 	Reflection *refl;
 	struct imagefeature *peak;
 	double Ih;   /* normalised */
+	struct panel *panel;  /* panel the reflection appears on
+                               * (we assume this never changes) */
 };
 
 static int pair_peaks(ImageFeatureList *flist, UnitCell *cell, RefList *reflist,
-                      struct reflpeak *rps)
+                      struct reflpeak *rps, struct detector *det)
 {
 	int i;
 	const double min_dist = 0.25;
@@ -102,6 +104,7 @@ static int pair_peaks(ImageFeatureList *flist, UnitCell *cell, RefList *reflist,
 
 			rps[n_acc].refl = refl;
 			rps[n_acc].peak = f;
+			rps[n_acc].panel = find_panel(det, f->fs, f->ss);
 			n_acc++;
 		}
 
@@ -112,10 +115,9 @@ static int pair_peaks(ImageFeatureList *flist, UnitCell *cell, RefList *reflist,
 
 
 static void twod_mapping(double fs, double ss, double *px, double *py,
-                         struct detector *det)
+                         struct panel *p)
 {
 	double xs, ys;
-	struct panel *p = find_panel(det, fs, ss);
 
 	xs = fs*p->fsx + ss*p->ssx;
 	ys = fs*p->fsy + ss*p->ssy;
@@ -216,13 +218,13 @@ static double pos_gradient(int param, struct reflpeak *rp, struct detector *det,
 	double fsh, ssh;
 	double tt, clen, azi, azf;
 
-	twod_mapping(rp->peak->fs, rp->peak->ss, &xpk, &ypk, det);
+	twod_mapping(rp->peak->fs, rp->peak->ss, &xpk, &ypk, rp->panel);
 	get_detector_pos(rp->refl, &fsh, &ssh);
-	twod_mapping(fsh, ssh, &xh, &yh, det);
+	twod_mapping(fsh, ssh, &xh, &yh, rp->panel);
 	get_indices(rp->refl, &h, &k, &l);
 
 	tt = asin(lambda * resolution(cell, h, k, l));
-	clen = find_panel(det, fsh, ssh)->clen;
+	clen = rp->panel->clen;
 	azi = atan2(yh, xh);
 	azf = 2.0*(cos(azi) + sin(azi));  /* FIXME: Why factor of 2? */
 
@@ -276,9 +278,9 @@ static double pos_dev(struct reflpeak *rp, struct detector *det)
 	/* Peak position term */
 	double xpk, ypk, xh, yh;
 	double fsh, ssh;
-	twod_mapping(rp->peak->fs, rp->peak->ss, &xpk, &ypk, det);
+	twod_mapping(rp->peak->fs, rp->peak->ss, &xpk, &ypk, rp->panel);
 	get_detector_pos(rp->refl, &fsh, &ssh);
-	twod_mapping(fsh, ssh, &xh, &yh, det);
+	twod_mapping(fsh, ssh, &xh, &yh, rp->panel);
 	return (xh-xpk) + (yh-ypk);
 }
 
@@ -427,7 +429,7 @@ int refine_prediction(struct image *image, Crystal *cr)
 	if ( rps == NULL ) return 1;
 
 	n = pair_peaks(image->features, crystal_get_cell(cr),
-	               crystal_get_reflections(cr), rps);
+	               crystal_get_reflections(cr), rps, image->det);
 	STATUS("%i peaks\n", n);
 	if ( n < 10 ) {
 		ERROR("Too few paired peaks to refine orientation.\n");
@@ -453,16 +455,10 @@ int refine_prediction(struct image *image, Crystal *cr)
 	/* Refine */
 	STATUS("Initial residual = %e\n", residual(rps, n, image->det));
 	for ( i=0; i<MAX_CYCLES; i++ ) {
-		int n_gain = 0;
-		int n_lost = 0;
-		double mpc;
 		iterate(rps, n, crystal_get_cell(cr), image);
-		update_partialities_2(cr, PMODEL_SCSPHERE, &n_gain, &n_lost,
-		                      &mpc);
+		update_partialities(cr, PMODEL_SCSPHERE);
 		STATUS("Residual after iteration %i = %e\n",
 		        i, residual(rps, n, image->det));
-		STATUS("%i gained, %i lost, mean p change = %e\n", n_gain,
-		       n_lost, mpc);
 	}
 
 	free(rps);
diff --git a/tests/integration_check.c b/tests/integration_check.c
index 2d0a6c7a..ccb613c3 100644
--- a/tests/integration_check.c
+++ b/tests/integration_check.c
@@ -120,7 +120,7 @@ int main(int argc, char *argv[])
 
 		list = reflist_new();
 		refl = add_refl(list, 0, 0, 0);
-		set_detector_pos(refl, 0.0, 64, 64);
+		set_detector_pos(refl, 64, 64);
 		cell = cell_new();
 		cell_set_lattice_type(cell, L_CUBIC);
 		cell_set_centering(cell, 'P');
diff --git a/tests/prediction_gradient_check.c b/tests/prediction_gradient_check.c
index 09b4a868..15b4d803 100644
--- a/tests/prediction_gradient_check.c
+++ b/tests/prediction_gradient_check.c
@@ -48,8 +48,7 @@ int checkrxy;
 
 
 static void scan(RefList *reflections, RefList *compare,
-                 int *valid, long double *vals[3], int idx,
-                 struct detector *det)
+                 int *valid, long double *vals[3], int idx)
 {
 	int i;
 	Reflection *refl;
@@ -64,6 +63,7 @@ static void scan(RefList *reflections, RefList *compare,
 		Reflection *refl2;
 		double rlow, rhigh, p;
 		double fs, ss, xh, yh;
+		struct panel *panel;
 
 		get_indices(refl, &h, &k, &l);
 		refl2 = find_refl(compare, h, k, l);
@@ -75,7 +75,8 @@ static void scan(RefList *reflections, RefList *compare,
 
 		get_partial(refl2, &rlow, &rhigh, &p);
 		get_detector_pos(refl2, &fs, &ss);
-		twod_mapping(fs, ss, &xh, &yh, det);
+		panel = get_panel(refl2);
+		twod_mapping(fs, ss, &xh, &yh, panel);
 
 		switch ( checkrxy ) {
 
@@ -169,16 +170,14 @@ static void calc_either_side(Crystal *cr, double incr_val,
 
 	cr_new = new_shifted_crystal(cr, refine, -incr_val);
 	compare = find_intersections(image, cr_new, PMODEL_SCSPHERE);
-	scan(crystal_get_reflections(cr), compare, valid, vals, 0,
-	     crystal_get_image(cr)->det);
+	scan(crystal_get_reflections(cr), compare, valid, vals, 0);
 	cell_free(crystal_get_cell(cr_new));
 	crystal_free(cr_new);
 	reflist_free(compare);
 
 	cr_new = new_shifted_crystal(cr, refine, +incr_val);
 	compare = find_intersections(image, cr_new, PMODEL_SCSPHERE);
-	scan(crystal_get_reflections(cr), compare, valid, vals, 2,
-	     crystal_get_image(cr)->det);
+	scan(crystal_get_reflections(cr), compare, valid, vals, 2);
 	cell_free(crystal_get_cell(cr_new));
 	crystal_free(cr_new);
 	reflist_free(compare);
@@ -231,8 +230,7 @@ static double test_gradients(Crystal *cr, double incr_val, int refine,
 	}
 	for ( i=0; i<nref; i++ ) valid[i] = 1;
 
-	scan(reflections, reflections, valid, vals, 1,
-	     crystal_get_image(cr)->det);
+	scan(reflections, reflections, valid, vals, 1);
 
 	calc_either_side(cr, incr_val, valid, vals, refine);
 
@@ -282,11 +280,15 @@ static double test_gradients(Crystal *cr, double incr_val, int refine,
 			} else {
 
 				struct imagefeature pk;
+				struct image *image;
 
 				pk.fs = 0.0;
 				pk.ss = 0.0;
+
+				image = crystal_get_image(cr);
 				rp.refl = refl;
 				rp.peak = &pk;
+				rp.panel = &image->det->panels[0];
 
 				cgrad = pos_gradient(refine, &rp,
 				                  crystal_get_image(cr)->det,