partialator: Move scaling calculation into PR proper

5 files changed, 70 insertions, 53 deletions

diff --git a/libcrystfel/src/geometry.h b/libcrystfel/src/geometry.h
index cecdc4f1..152c0e47 100644
--- a/libcrystfel/src/geometry.h
+++ b/libcrystfel/src/geometry.h
@@ -77,7 +77,9 @@ enum gparam {
 	GPARAM_DIV,
 	GPARAM_DETX,
 	GPARAM_DETY,
-	GPARAM_CLEN
+	GPARAM_CLEN,
+	GPARAM_OSF,
+	GPARAM_BFAC
 };
 
 
diff --git a/src/partialator.c b/src/partialator.c
index 9bdb684b..19449590 100644
--- a/src/partialator.c
+++ b/src/partialator.c
@@ -85,6 +85,7 @@ struct refine_args
 	RefList *full;
 	Crystal *crystal;
 	PartialityModel pmodel;
+	int no_scale;
 	struct prdata prdata;
 };
 
@@ -105,7 +106,8 @@ static void refine_image(void *task, int id)
 	struct refine_args *pargs = task;
 	Crystal *cr = pargs->crystal;
 
-	pargs->prdata = pr_refine(cr, pargs->full, pargs->pmodel);
+	pargs->prdata = pr_refine(cr, pargs->full, pargs->pmodel,
+	                          pargs->no_scale);
 }
 
 
@@ -148,10 +150,6 @@ static void refine_all(Crystal **crystals, int n_crystals,
 	struct refine_args task_defaults;
 	struct queue_args qargs;
 
-	/* If the partiality model is "p=1", this refinement is really, really
-	 * easy... */
-	if ( pmodel == PMODEL_UNITY ) return;
-
 	task_defaults.full = full;
 	task_defaults.crystal = NULL;
 	task_defaults.pmodel = pmodel;
@@ -641,15 +639,8 @@ int main(int argc, char *argv[])
 //	early_rejection(crystals, n_crystals);
 
 	/* Make initial estimates */
-	STATUS("Performing initial scaling.\n");
-	if ( noscale ) {
-		STATUS("Skipping scaling step (--no-scale).\n");
-		full = lsq_intensities(crystals, n_crystals, nthreads, pmodel,
-		                       min_measurements);
-	} else {
-		full = scale_intensities(crystals, n_crystals, nthreads, pmodel,
-		                         min_measurements);
-	}
+	full = lsq_intensities(crystals, n_crystals, nthreads, pmodel,
+	                       min_measurements);
 
 	check_rejection(crystals, n_crystals);
 
@@ -668,11 +659,11 @@ int main(int argc, char *argv[])
 	/* Iterate */
 	for ( i=0; i<n_iter; i++ ) {
 
-		STATUS("Post refinement cycle %i of %i\n", i+1, n_iter);
+		STATUS("Refinement cycle %i of %i\n", i+1, n_iter);
 
 		srdata.n_filtered = 0;
 
-		/* Refine the geometry of all patterns to get the best fit */
+		/* Refine all crystals to get the best fit */
 		refine_all(crystals, n_crystals, full, nthreads, pmodel,
 		           &srdata);
 
@@ -681,14 +672,8 @@ int main(int argc, char *argv[])
 
 		/* Re-estimate all the full intensities */
 		reflist_free(full);
-		if ( noscale ) {
-			STATUS("Skipping scaling step (--no-scale).\n");
-			full = lsq_intensities(crystals, n_crystals, nthreads,
-			                       pmodel, min_measurements);
-		} else {
-			full = scale_intensities(crystals, n_crystals, nthreads,
-			                         pmodel, min_measurements);
-		}
+		full = lsq_intensities(crystals, n_crystals, nthreads,
+		                       pmodel, min_measurements);
 
 		check_rejection(crystals, n_crystals);
 
diff --git a/src/post-refinement.c b/src/post-refinement.c
index 8e6f6108..5a1053b8 100644
--- a/src/post-refinement.c
+++ b/src/post-refinement.c
@@ -174,14 +174,16 @@ static double volume_fraction(double rlow, double rhigh, double pr,
 }
 
 
-/* Return the gradient of partiality wrt parameter 'k' given the current status
- * of 'image'. */
-double p_gradient(Crystal *cr, int k, Reflection *refl, PartialityModel pmodel)
+/* Return the gradient of p/G wrt parameter 'k' given the current status
+ * of the crystal. */
+double gradient(Crystal *cr, int k, Reflection *refl, PartialityModel pmodel)
 {
 	double glow, ghigh;
 	double rlow, rhigh, p;
 	struct image *image = crystal_get_image(cr);
 	double R = crystal_get_profile_radius(cr);
+	double osf = crystal_get_osf(cr);
+	double gr;
 
 	get_partial(refl, &rlow, &rhigh, &p);
 
@@ -196,7 +198,8 @@ double p_gradient(Crystal *cr, int k, Reflection *refl, PartialityModel pmodel)
 		Rglow = volume_fraction_rgradient(rlow, R, pmodel);
 		Rghigh = volume_fraction_rgradient(rhigh, R, pmodel);
 
-		return 4.0*psph/(3.0*D) + (4.0*R/(3.0*D))*(Rglow - Rghigh);
+		gr = 4.0*psph/(3.0*D) + (4.0*R/(3.0*D))*(Rglow - Rghigh);
+		return gr / osf;
 
 	}
 
@@ -214,11 +217,17 @@ double p_gradient(Crystal *cr, int k, Reflection *refl, PartialityModel pmodel)
 		get_symmetric_indices(refl, &hs, &ks, &ls);
 		ds = 2.0 * resolution(crystal_get_cell(cr), hs, ks, ls);
 
-		return (ds/2.0)*(glow+ghigh) - 4.0*R*psph*ds/(3.0*D*D);
+		gr = (ds/2.0)*(glow+ghigh) - 4.0*R*psph*ds/(3.0*D*D);
+		return gr / osf;
 
 	}
 
-	return r_gradient(crystal_get_cell(cr), k, refl, image) * (glow-ghigh);
+	if ( k == GPARAM_OSF ) {
+		return -p/(osf*osf);
+	}
+
+	gr = r_gradient(crystal_get_cell(cr), k, refl, image) * (glow-ghigh);
+	return gr / osf;
 }
 
 
@@ -274,6 +283,18 @@ static void apply_shift(Crystal *cr, int k, double shift)
 		crystal_set_profile_radius(cr, t);
 		break;
 
+		case GPARAM_OSF :
+		if ( isnan(shift) ) {
+			ERROR("Refusing nan shift of OSF\n");
+		} else if (shift > 100.0 ) {
+			ERROR("Refusing large (%e) OSF shift.\n", shift);
+		} else {
+			t = crystal_get_osf(cr);
+			t += shift;
+			crystal_set_osf(cr, t);
+		}
+		break;
+
 		case GPARAM_ASX :
 		case GPARAM_ASY :
 		case GPARAM_ASZ :
@@ -296,7 +317,7 @@ static void apply_shift(Crystal *cr, int k, double shift)
 
 /* Perform one cycle of post refinement on 'image' against 'full' */
 static double pr_iterate(Crystal *cr, const RefList *full,
-                         PartialityModel pmodel, int *n_filtered)
+                         PartialityModel pmodel, int no_scale, int *n_filtered)
 {
 	gsl_matrix *M;
 	gsl_vector *v;
@@ -310,6 +331,7 @@ static double pr_iterate(Crystal *cr, const RefList *full,
 	const int verbose = 0;
 	int num_params = 0;
 	enum gparam rv[32];
+	double G;
 
 	*n_filtered = 0;
 
@@ -327,13 +349,19 @@ static double pr_iterate(Crystal *cr, const RefList *full,
 		rv[num_params++] = GPARAM_CSZ;
 	}
 
-	STATUS("Refining %i parameters.\n", num_params);
+	/* If we are scaling, refine scale factors (duh) */
+	if ( !no_scale ) {
+		rv[num_params++] = GPARAM_OSF;
+		//rv[num_params++] = GPARAM_BFAC;  FIXME: Next
+	}
 
 	reflections = crystal_get_reflections(cr);
 
 	M = gsl_matrix_calloc(num_params, num_params);
 	v = gsl_vector_calloc(num_params);
 
+	G = crystal_get_osf(cr);
+
 	/* Construct the equations, one per reflection in this image */
 	for ( refl = first_refl(reflections, &iter);
 	      refl != NULL;
@@ -344,7 +372,7 @@ static double pr_iterate(Crystal *cr, const RefList *full,
 		double w;
 		double I_partial;
 		int k;
-		double p, l;
+		double p, L;
 		Reflection *match;
 		double gradients[num_params];
 
@@ -359,19 +387,22 @@ static double pr_iterate(Crystal *cr, const RefList *full,
 
 		I_full = get_intensity(match);
 
-		/* Actual measurement of this reflection from this pattern? */
-		I_partial = get_intensity(refl) / crystal_get_osf(cr);
 		p = get_partiality(refl);
-		l = get_lorentz(refl);
+		L = get_lorentz(refl);
+
+		/* Actual measurement of this reflection from this pattern? */
+		I_partial = get_intensity(refl);
 
 		/* Calculate the weight for this reflection */
-		w =  pow(get_esd_intensity(refl), 2.0);
-		w += l * p * I_full * pow(get_esd_intensity(match), 2.0);
-		w = pow(w, -1.0);
+		//w = pow(get_esd_intensity(refl), 2.0);
+		//w += (p/L) * I_full * pow(get_esd_intensity(match), 2.0);
+		//w = pow(w, -1.0);
+		w = 1.0;
 
 		/* Calculate all gradients for this reflection */
 		for ( k=0; k<num_params; k++ ) {
-			gradients[k] = p_gradient(cr, rv[k], refl, pmodel) * l;
+			gradients[k] = gradient(cr, rv[k], refl, pmodel);
+			gradients[k] *= I_full / L;
 		}
 
 		for ( k=0; k<num_params; k++ ) {
@@ -386,8 +417,7 @@ static double pr_iterate(Crystal *cr, const RefList *full,
 				/* Matrix is symmetric */
 				if ( g > k ) continue;
 
-				M_c = gradients[g] * gradients[k];
-				M_c *= w * pow(I_full, 2.0);
+				M_c = w * gradients[g] * gradients[k];
 
 				M_curr = gsl_matrix_get(M, k, g);
 				gsl_matrix_set(M, k, g, M_curr + M_c);
@@ -395,8 +425,8 @@ static double pr_iterate(Crystal *cr, const RefList *full,
 
 			}
 
-			delta_I = I_partial - (l * p * I_full);
-			v_c = w * delta_I * I_full * gradients[k];
+			delta_I = I_partial - (p * I_full / (L*G));
+			v_c = w * delta_I * gradients[k];
 			v_curr = gsl_vector_get(v, k);
 			gsl_vector_set(v, k, v_curr + v_c);
 
@@ -486,7 +516,7 @@ static double guide_dev(Crystal *cr, const RefList *full)
 
 
 struct prdata pr_refine(Crystal *cr, const RefList *full,
-                        PartialityModel pmodel)
+                        PartialityModel pmodel, int no_scale)
 {
 	double dev;
 	int i;
@@ -518,7 +548,7 @@ struct prdata pr_refine(Crystal *cr, const RefList *full,
 		cell_get_reciprocal(crystal_get_cell(cr), &asx, &asy, &asz,
 			               &bsx, &bsy, &bsz, &csx, &csy, &csz);
 
-		pr_iterate(cr, full, pmodel, &prdata.n_filtered);
+		pr_iterate(cr, full, pmodel, no_scale, &prdata.n_filtered);
 
 		update_partialities(cr, pmodel);
 
diff --git a/src/post-refinement.h b/src/post-refinement.h
index caf5f9af..844afbd1 100644
--- a/src/post-refinement.h
+++ b/src/post-refinement.h
@@ -3,11 +3,11 @@
  *
  * Post refinement
  *
- * Copyright © 2012-2014 Deutsches Elektronen-Synchrotron DESY,
+ * Copyright © 2012-2015 Deutsches Elektronen-Synchrotron DESY,
  *                       a research centre of the Helmholtz Association.
  *
  * Authors:
- *   2010-2014 Thomas White <taw@physics.org>
+ *   2010-2015 Thomas White <taw@physics.org>
  *
  * This file is part of CrystFEL.
  *
@@ -51,10 +51,10 @@ struct prdata
 
 
 extern struct prdata pr_refine(Crystal *cr, const RefList *full,
-                               PartialityModel pmodel);
+                               PartialityModel pmodel, int no_scale);
 
 /* Exported so it can be poked by tests/pr_p_gradient_check */
-extern double p_gradient(Crystal *cr, int k, Reflection *refl,
-                         PartialityModel pmodel);
+extern double gradient(Crystal *cr, int k, Reflection *refl,
+                       PartialityModel pmodel);
 
 #endif	/* POST_REFINEMENT_H */
diff --git a/tests/pr_p_gradient_check.c b/tests/pr_p_gradient_check.c
index affddeab..1ae7b049 100644
--- a/tests/pr_p_gradient_check.c
+++ b/tests/pr_p_gradient_check.c
@@ -303,7 +303,7 @@ static double test_gradients(Crystal *cr, double incr_val, int refine,
 			grad = (grad1 + grad2) / 2.0;
 			i++;
 
-			cgrad = p_gradient(cr, refine, refl, pmodel);
+			cgrad = gradient(cr, refine, refl, pmodel);
 
 			get_partial(refl, &r1, &r2, &p);