Scaling fixes

5 files changed, 71 insertions, 36 deletions

diff --git a/src/merge.c b/src/merge.c
index 33221bda..099bad2d 100644
--- a/src/merge.c
+++ b/src/merge.c
@@ -189,7 +189,7 @@ static void run_merge_job(void *vwargs, int cookie)
 		}
 
 		/* Total (multiplicative) correction factor */
-		corr = (1.0/G) * exp(B*res*res) * get_lorentz(refl)
+		corr = G * exp(B*res*res) * get_lorentz(refl)
 		        / get_partiality(refl);
 		if ( isnan(corr) ) {
 			ERROR("NaN corr:\n");
diff --git a/src/partialator.c b/src/partialator.c
index 46a5a25f..124ce2e4 100644
--- a/src/partialator.c
+++ b/src/partialator.c
@@ -637,7 +637,7 @@ static void write_to_pgraph(FILE *fh, RefList *list, RefList *full, Crystal *cr,
 		pcalc = get_partiality(refl);
 
 		/* Observed partiality */
-		corr = (1.0/G) * exp(B*res*res) * get_lorentz(refl);
+		corr = G * exp(B*res*res) * get_lorentz(refl);
 		Ipart = get_intensity(refl) * corr;
 		pobs = Ipart / get_intensity(match);
 
@@ -689,7 +689,7 @@ static void write_specgraph(RefList *full, Crystal *crystal, int in)
 		match = find_refl(full, h, k, l);
 		if ( match == NULL ) continue;
 
-		corr = (1.0/G) * exp(B*res*res) * get_lorentz(refl);
+		corr = G * exp(B*res*res) * get_lorentz(refl);
 		Ipart = get_intensity(refl) * corr;
 		Ifull = get_intensity(match);
 		esd = get_esd_intensity(match);
@@ -847,7 +847,6 @@ int main(int argc, char *argv[])
 	double max_B = 1e-18;
 	char *rfile = NULL;
 	RefList *reference = NULL;
-	RefList *r;
 
 	/* Long options */
 	const struct option longopts[] = {
@@ -1223,13 +1222,22 @@ int main(int argc, char *argv[])
 	//early_rejection(crystals, n_crystals);
 
 	/* Initial rejection, figures of merit etc */
-	full = merge_intensities(crystals, n_crystals, nthreads, pmodel,
-	                         min_measurements, push_res, 1);
+	if ( reference == NULL ) {
+		scale_all(crystals, n_crystals, nthreads, pmodel);
+		full = merge_intensities(crystals, n_crystals, nthreads, pmodel,
+		                         min_measurements, push_res, 1);
+	} else {
+		full = reference;
+	}
 	check_rejection(crystals, n_crystals, full, max_B);
-	r = (reference != NULL) ? reference : full;
-	show_all_residuals(crystals, n_crystals, r);
-	write_pgraph(r, crystals, n_crystals, 0);
-	write_specgraph(r, crystals[0], 0);
+
+	if ( !no_scale ) {
+		scale_all_to_reference(crystals, n_crystals, full);
+	}
+
+	show_all_residuals(crystals, n_crystals, full);
+	write_pgraph(full, crystals, n_crystals, 0);
+	write_specgraph(full, crystals[0], 0);
 
 	/* Iterate */
 	for ( i=0; i<n_iter; i++ ) {
@@ -1237,10 +1245,7 @@ int main(int argc, char *argv[])
 		STATUS("Scaling and refinement cycle %i of %i\n", i+1, n_iter);
 
 		if ( !no_pr ) {
-
-			r = (reference != NULL) ? reference : full;
-			refine_all(crystals, n_crystals, r, nthreads, pmodel);
-
+			refine_all(crystals, n_crystals, full, nthreads, pmodel);
 		}
 
 		if ( !no_scale ) {
@@ -1248,8 +1253,6 @@ int main(int argc, char *argv[])
 			if ( reference == NULL ) {
 
 				/* No reference -> XSCALE-like algorithm */
-				scale_all(crystals, n_crystals, nthreads, pmodel);
-				reflist_free(full);
 				full = merge_intensities(crystals, n_crystals,
 				                         nthreads, pmodel,
 				                         min_measurements,
@@ -1258,21 +1261,27 @@ int main(int argc, char *argv[])
 			} else {
 
 				/* Reference -> Ginn-style linear algorithm */
-				scale_all_to_reference(crystals, n_crystals,
-				                       reference);
+				if ( !no_scale ) {
+					scale_all_to_reference(crystals, n_crystals,
+					                       reference);
+				}
 
 			}
 		}
 
 		check_rejection(crystals, n_crystals, full, max_B);
-		reflist_free(full);
-		full = merge_intensities(crystals, n_crystals, nthreads,
-		                         pmodel, min_measurements,
-		                         push_res, 1);
-		r = (reference != NULL) ? reference : full;
-		show_all_residuals(crystals, n_crystals, r);
-		write_pgraph(r, crystals, n_crystals, i+1);
-		write_specgraph(r, crystals[0], i+1);
+
+		if ( reference == NULL ) {
+			scale_all(crystals, n_crystals, nthreads, pmodel);
+			reflist_free(full);
+			full = merge_intensities(crystals, n_crystals, nthreads,
+			                         pmodel, min_measurements,
+			                         push_res, 1);
+		} /* else full still equals reference */
+
+		show_all_residuals(crystals, n_crystals, full);
+		write_pgraph(full, crystals, n_crystals, i+1);
+		write_specgraph(full, crystals[0], i+1);
 
 		if ( output_everycycle ) {
 
@@ -1306,12 +1315,32 @@ int main(int argc, char *argv[])
 		}
 	}
 
-	full = merge_intensities(crystals, n_crystals, nthreads, pmodel,
-	                         min_measurements, push_res, 1);
-	r = (reference != NULL) ? reference : full;
-	show_all_residuals(crystals, n_crystals, r);
-	write_pgraph(r, crystals, n_crystals, n_iter+1);
-	write_specgraph(r, crystals[0], n_iter+1);
+
+	if ( reference == NULL ) {
+		scale_all(crystals, n_crystals, nthreads, pmodel);
+		reflist_free(full);
+		full = merge_intensities(crystals, n_crystals, nthreads,
+		                         pmodel, min_measurements,
+		                         push_res, 1);
+	} /* else full still equals reference */
+	if ( !no_scale ) {
+		scale_all_to_reference(crystals, n_crystals, full);
+	}
+	show_all_residuals(crystals, n_crystals, full);
+	write_pgraph(full, crystals, n_crystals, n_iter+1);
+	write_specgraph(full, crystals[0], n_iter+1);
+	//STATUS("Final profile radius: %e\n", crystal_get_profile_radius(crystals[0]));
+
+	/* If we've been using a reference up to now, it's time to actually
+	 * scale and merge the final version */
+	if ( reference != NULL ) {
+		STATUS("Starting final reference-free merge\n");
+		scale_all_to_reference(crystals, n_crystals, reference);
+		STATUS("Scaling done\n");
+		full = merge_intensities(crystals, n_crystals, nthreads,
+		                         pmodel, min_measurements, push_res, 1);
+		STATUS("Done\n");
+	} /* else we just did it */
 
 	/* Output results */
 	STATUS("Writing overall results to %s\n", outfile);
diff --git a/src/post-refinement.c b/src/post-refinement.c
index e1a2e05a..8e2241f6 100644
--- a/src/post-refinement.c
+++ b/src/post-refinement.c
@@ -92,7 +92,7 @@ double residual(Crystal *cr, const RefList *full, int free,
 	double B = crystal_get_Bfac(cr);
 	UnitCell *cell = crystal_get_cell(cr);
 
-	if ( linear_scale(crystal_get_reflections(cr), full, &G) ) {
+	if ( linear_scale(full, crystal_get_reflections(cr), &G) ) {
 		return INFINITY;
 	}
 
diff --git a/src/rejection.c b/src/rejection.c
index a44b697f..9d41f9b4 100644
--- a/src/rejection.c
+++ b/src/rejection.c
@@ -137,7 +137,7 @@ static void check_cc(Crystal *cr, RefList *full)
 		pcalc = get_partiality(refl);
 
 		/* Observed partiality */
-		corr = (1.0/G) * exp(B*res*res) * get_lorentz(refl);
+		corr = G * exp(B*res*res) * get_lorentz(refl);
 		Ipart = get_intensity(refl) * corr;
 		pobs = Ipart / get_intensity(match);
 
diff --git a/src/scaling.c b/src/scaling.c
index 97503b89..ab913229 100644
--- a/src/scaling.c
+++ b/src/scaling.c
@@ -515,9 +515,9 @@ int linear_scale(const RefList *list1, const RefList *list2, double *G)
 			}
 		}
 
-		x[n] = Ih2;
+		x[n] = Ih2 / get_partiality(refl2);
 		y[n] = Ih1;
-		w[n] = get_partiality(refl1);
+		w[n] = get_partiality(refl2);
 		n++;
 
 	}
@@ -534,6 +534,12 @@ int linear_scale(const RefList *list1, const RefList *list2, double *G)
 		return 1;
 	}
 
+	if ( isnan(*G) ) {
+		ERROR("Scaling gave NaN (%i pairs)\n", n);
+		abort();
+		return 1;
+	}
+
 	free(x);
 	free(y);
 	free(w);