Use linear scale when scaling against reference

4 files changed, 147 insertions, 106 deletions

diff --git a/src/partialator.c b/src/partialator.c
index 49605bf9..af04d935 100644
--- a/src/partialator.c
+++ b/src/partialator.c
@@ -654,7 +654,8 @@ static void write_specgraph(RefList *full, Crystal *crystal, int in)
 	char tmp[256];
 	Reflection *refl;
 	RefListIterator *iter;
-	double G, B;
+	double G = crystal_get_osf(crystal);
+	double B = crystal_get_Bfac(crystal);
 	UnitCell *cell;
 	struct image *image = crystal_get_image(crystal);
 
@@ -671,8 +672,6 @@ static void write_specgraph(RefList *full, Crystal *crystal, int in)
 
 	fprintf(fh, "khalf/m   pcalc    pobs\n");
 
-	G = crystal_get_osf(crystal);
-	B = crystal_get_Bfac(crystal);
 	cell = crystal_get_cell(crystal);
 
 	for ( refl = first_refl(crystal_get_reflections(crystal), &iter);
@@ -848,6 +847,7 @@ 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[] = {
@@ -1226,8 +1226,10 @@ int main(int argc, char *argv[])
 	full = merge_intensities(crystals, n_crystals, nthreads, pmodel,
 	                         min_measurements, push_res, 1);
 	check_rejection(crystals, n_crystals, full, max_B);
-	show_all_residuals(crystals, n_crystals, full);
-	write_pgraph(full, crystals, n_crystals, 0);
+	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);
 
 	/* Iterate */
 	for ( i=0; i<n_iter; i++ ) {
@@ -1235,22 +1237,36 @@ int main(int argc, char *argv[])
 		STATUS("Scaling and refinement cycle %i of %i\n", i+1, n_iter);
 
 		if ( !no_scale ) {
-			scale_all(crystals, n_crystals, nthreads, pmodel);
-			reflist_free(full);
-			full = merge_intensities(crystals, n_crystals, nthreads,
-			                         pmodel, min_measurements,
-			                         push_res, 1);
+
+			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,
+				                         push_res, 1);
+
+			} else {
+
+				/* Reference -> Ginn-style linear algorithm */
+				scale_all_to_reference(crystals, n_crystals,
+				                       reference);
+
+			}
 		}
 
 		if ( !no_pr ) {
 
-			RefList *r = (reference != NULL) ? reference : full;
+			r = (reference != NULL) ? reference : full;
 			refine_all(crystals, n_crystals, r, nthreads, pmodel);
 
 			reflist_free(full);
 			full = merge_intensities(crystals, n_crystals, nthreads,
 			                         pmodel, min_measurements,
 			                         push_res, 1);
+
 		}
 
 		check_rejection(crystals, n_crystals, full, max_B);
@@ -1258,9 +1274,10 @@ int main(int argc, char *argv[])
 		full = merge_intensities(crystals, n_crystals, nthreads,
 		                         pmodel, min_measurements,
 		                         push_res, 1);
-		show_all_residuals(crystals, n_crystals, full);
-		write_pgraph(full, crystals, n_crystals, i+1);
-		write_specgraph(full, crystals[0], i+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 ( output_everycycle ) {
 
@@ -1296,8 +1313,10 @@ int main(int argc, char *argv[])
 
 	full = merge_intensities(crystals, n_crystals, nthreads, pmodel,
 	                         min_measurements, push_res, 1);
-	show_all_residuals(crystals, n_crystals, full);
-	write_pgraph(full, crystals, n_crystals, n_iter+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);
 
 	/* Output results */
 	STATUS("Writing overall results to %s\n", outfile);
diff --git a/src/post-refinement.c b/src/post-refinement.c
index 165080b4..5624b15d 100644
--- a/src/post-refinement.c
+++ b/src/post-refinement.c
@@ -44,6 +44,7 @@
 #include "cell.h"
 #include "cell-utils.h"
 #include "reflist-utils.h"
+#include "scaling.h"
 
 
 struct prdata
@@ -79,104 +80,23 @@ const char *str_prflag(enum prflag flag)
 }
 
 
-static int linear_scale(RefList *list1, const RefList *list2, double *G)
-{
-	Reflection *refl1;
-	Reflection *refl2;
-	RefListIterator *iter;
-	int max_n = 256;
-	int n = 0;
-	double *x;
-	double *y;
-	double *w;
-	int r;
-	double cov11;
-	double sumsq;
-
-	x = malloc(max_n*sizeof(double));
-	w = malloc(max_n*sizeof(double));
-	y = malloc(max_n*sizeof(double));
-	if ( (x==NULL) || (y==NULL) || (w==NULL) ) {
-		ERROR("Failed to allocate memory for scaling.\n");
-		return 1;
-	}
-
-	for ( refl1 = first_refl(list1, &iter);
-	      refl1 != NULL;
-	      refl1 = next_refl(refl1, iter) )
-	{
-		signed int h, k, l;
-		double Ih1, Ih2;
-
-		get_indices(refl1, &h, &k, &l);
-		refl2 = find_refl(list2, h, k, l);
-		if ( refl2 == NULL ) continue;
-
-		Ih1 = get_intensity(refl1);
-		Ih2 = get_intensity(refl2);
-
-		if ( (Ih1 <= 0.0) || (Ih2 <= 0.0) ) continue;
-		if ( isnan(Ih1) || isinf(Ih1) ) continue;
-		if ( isnan(Ih2) || isinf(Ih2) ) continue;
-
-		if ( n == max_n ) {
-			max_n *= 2;
-			x = realloc(x, max_n*sizeof(double));
-			y = realloc(y, max_n*sizeof(double));
-			w = realloc(w, max_n*sizeof(double));
-			if ( (x==NULL) || (y==NULL) || (w==NULL) ) {
-				ERROR("Failed to allocate memory for scaling.\n");
-				return 1;
-			}
-		}
-
-		x[n] = Ih2;
-		y[n] = Ih1;
-		w[n] = get_partiality(refl1);
-		n++;
-
-	}
-
-	if ( n < 2 ) {
-		ERROR("Not enough reflections for scaling\n");
-		return 1;
-	}
-
-	r = gsl_fit_wmul(x, 1, w, 1, y, 1, n, G, &cov11, &sumsq);
-
-	if ( r ) {
-		ERROR("Scaling failed.\n");
-		return 1;
-	}
-
-	free(x);
-	free(y);
-	free(w);
-
-	return 0;
-}
-
-
-
 double residual(Crystal *cr, const RefList *full, int free,
                 int *pn_used, const char *filename)
 {
-	double G;
 	Reflection *refl;
 	RefListIterator *iter;
 	int n_used = 0;
 	double num = 0.0;
 	double den = 0.0;
-
-	if ( linear_scale(crystal_get_reflections(cr), full, &G) ) {
-		return INFINITY;
-	}
+	double G = crystal_get_osf(cr);
+	double B = crystal_get_Bfac(cr);
+	UnitCell *cell = crystal_get_cell(cr);
 
 	for ( refl = first_refl(crystal_get_reflections(cr), &iter);
 	      refl != NULL;
 	      refl = next_refl(refl, iter) )
 	{
-		double p, w;
+		double p, w, corr, res;
 		signed int h, k, l;
 		Reflection *match;
 		double I_full;
@@ -185,6 +105,7 @@ double residual(Crystal *cr, const RefList *full, int free,
 		if ( free && !get_flag(refl) ) continue;
 
 		get_indices(refl, &h, &k, &l);
+		res = resolution(cell, h, k, l);
 		match = find_refl(full, h, k, l);
 		if ( match == NULL ) continue;
 		I_full = get_intensity(match);
@@ -194,7 +115,8 @@ double residual(Crystal *cr, const RefList *full, int free,
 		p = get_partiality(refl);
 		//if ( p < 0.2 ) continue;
 
-		int1 = get_intensity(refl) / G;
+		corr = exp(-G) * exp(B*res*res) * get_lorentz(refl);
+		int1 = get_intensity(refl) * corr;
 		int2 = p*I_full;
 		w = 1.0;
 
diff --git a/src/scaling.c b/src/scaling.c
index 77a64384..0908cfb3 100644
--- a/src/scaling.c
+++ b/src/scaling.c
@@ -3,11 +3,11 @@
  *
  * Scaling
  *
- * Copyright © 2012-2015 Deutsches Elektronen-Synchrotron DESY,
+ * Copyright © 2012-2017 Deutsches Elektronen-Synchrotron DESY,
  *                       a research centre of the Helmholtz Association.
  *
  * Authors:
- *   2010-2015 Thomas White <taw@physics.org>
+ *   2010-2017 Thomas White <taw@physics.org>
  *
  * This file is part of CrystFEL.
  *
@@ -38,6 +38,7 @@
 #include <gsl/gsl_linalg.h>
 #include <gsl/gsl_eigen.h>
 #include <gsl/gsl_blas.h>
+#include <gsl/gsl_fit.h>
 
 #include "merge.h"
 #include "post-refinement.h"
@@ -456,3 +457,99 @@ void scale_all(Crystal **crystals, int n_crystals, int nthreads,
 		ERROR("Too many iterations - giving up!\n");
 	}
 }
+
+
+int linear_scale(RefList *list1, const RefList *list2, double *G)
+{
+	Reflection *refl1;
+	Reflection *refl2;
+	RefListIterator *iter;
+	int max_n = 256;
+	int n = 0;
+	double *x;
+	double *y;
+	double *w;
+	int r;
+	double cov11;
+	double sumsq;
+
+	x = malloc(max_n*sizeof(double));
+	w = malloc(max_n*sizeof(double));
+	y = malloc(max_n*sizeof(double));
+	if ( (x==NULL) || (y==NULL) || (w==NULL) ) {
+		ERROR("Failed to allocate memory for scaling.\n");
+		return 1;
+	}
+
+	int nb = 0;
+	for ( refl1 = first_refl(list1, &iter);
+	      refl1 != NULL;
+	      refl1 = next_refl(refl1, iter) )
+	{
+		signed int h, k, l;
+		double Ih1, Ih2;
+		nb++;
+
+		get_indices(refl1, &h, &k, &l);
+		refl2 = find_refl(list2, h, k, l);
+		if ( refl2 == NULL ) {
+			continue;
+		}
+
+		Ih1 = get_intensity(refl1);
+		Ih2 = get_intensity(refl2);
+
+		if ( (Ih1 <= 0.0) || (Ih2 <= 0.0) ) continue;
+		if ( isnan(Ih1) || isinf(Ih1) ) continue;
+		if ( isnan(Ih2) || isinf(Ih2) ) continue;
+
+		if ( n == max_n ) {
+			max_n *= 2;
+			x = realloc(x, max_n*sizeof(double));
+			y = realloc(y, max_n*sizeof(double));
+			w = realloc(w, max_n*sizeof(double));
+			if ( (x==NULL) || (y==NULL) || (w==NULL) ) {
+				ERROR("Failed to allocate memory for scaling.\n");
+				return 1;
+			}
+		}
+
+		x[n] = Ih2;
+		y[n] = Ih1;
+		w[n] = get_partiality(refl1);
+		n++;
+
+	}
+
+	if ( n < 2 ) {
+		ERROR("Not enough reflections for scaling (had %i, but %i remain)\n", nb, n);
+		return 1;
+	}
+
+	r = gsl_fit_wmul(x, 1, w, 1, y, 1, n, G, &cov11, &sumsq);
+
+	if ( r ) {
+		ERROR("Scaling failed.\n");
+		return 1;
+	}
+
+	free(x);
+	free(y);
+	free(w);
+
+	return 0;
+}
+
+
+void scale_all_to_reference(Crystal **crystals, int n_crystals,
+                            RefList *reference)
+{
+	int i;
+
+	for ( i=0; i<n_crystals; i++ ) {
+		double G;
+		linear_scale(crystal_get_reflections(crystals[i]), reference, &G);
+		crystal_set_osf(crystals[i], -log(G));
+		crystal_set_Bfac(crystals[i], 0.0);
+	}
+}
diff --git a/src/scaling.h b/src/scaling.h
index 534043d9..7b055003 100644
--- a/src/scaling.h
+++ b/src/scaling.h
@@ -3,11 +3,11 @@
  *
  * Scaling
  *
- * Copyright © 2012-2015 Deutsches Elektronen-Synchrotron DESY,
+ * Copyright © 2012-2017 Deutsches Elektronen-Synchrotron DESY,
  *                       a research centre of the Helmholtz Association.
  *
  * Authors:
- *   2010-2015 Thomas White <taw@physics.org>
+ *   2010-2017 Thomas White <taw@physics.org>
  *
  * This file is part of CrystFEL.
  *
@@ -44,4 +44,7 @@ extern double log_residual(Crystal *cr, const RefList *full, int free,
 extern void scale_all(Crystal **crystals, int n_crystals, int nthreads,
                       PartialityModel pmodel);
 
+extern void scale_all_to_reference(Crystal **crystals, int n_crystals,
+                                   RefList *reference);
+
 #endif	/* SCALING_H */