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.

3 files changed, 15 insertions, 40 deletions

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);