1 files changed, 171 insertions, 290 deletions

diff --git a/src/scaling.c b/src/scaling.c
index 77a64384..cabc5952 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,290 +38,22 @@
 #include <gsl/gsl_linalg.h>
 #include <gsl/gsl_eigen.h>
 #include <gsl/gsl_blas.h>
+#include <gsl/gsl_fit.h>
+#include <gsl/gsl_statistics_double.h>
 
 #include "merge.h"
 #include "post-refinement.h"
 #include "symmetry.h"
 #include "cell.h"
 #include "cell-utils.h"
-
-
-/* Maximum number of iterations of NLSq to do for each image per macrocycle. */
-#define MAX_CYCLES (10)
-
-
-/* Apply the given shift to the 'k'th parameter of 'image'. */
-static void apply_shift(Crystal *cr, int k, double shift)
-{
-	double t;
-
-	switch ( k ) {
-
-		case GPARAM_BFAC :
-		t = crystal_get_Bfac(cr);
-		t += shift;
-		crystal_set_Bfac(cr, t);
-		break;
-
-		case GPARAM_OSF :
-		t = crystal_get_osf(cr);
-		t += shift;
-		crystal_set_osf(cr, t);
-		break;
-
-		default :
-		ERROR("No shift defined for parameter %i\n", k);
-		abort();
-
-	}
-}
-
-
-/* Perform one cycle of scaling of 'cr' against 'full' */
-static double scale_iterate(Crystal *cr, const RefList *full,
-                            PartialityModel pmodel, int *nr)
-{
-	gsl_matrix *M;
-	gsl_vector *v;
-	gsl_vector *shifts;
-	int param;
-	Reflection *refl;
-	RefListIterator *iter;
-	RefList *reflections;
-	double max_shift;
-	int nref = 0;
-	int num_params = 0;
-	enum gparam rv[32];
-	double G, B;
-
-	*nr = 0;
-
-	rv[num_params++] = GPARAM_OSF;
-	rv[num_params++] = GPARAM_BFAC;
-
-	M = gsl_matrix_calloc(num_params, num_params);
-	v = gsl_vector_calloc(num_params);
-
-	reflections = crystal_get_reflections(cr);
-	G = crystal_get_osf(cr);
-	B = crystal_get_Bfac(cr);
-
-	/* Scaling terms */
-	for ( refl = first_refl(reflections, &iter);
-	      refl != NULL;
-	      refl = next_refl(refl, iter) )
-	{
-		signed int ha, ka, la;
-		double I_full, delta_I, esd;
-		double w;
-		double I_partial;
-		int k;
-		double p, L, s;
-		double fx;
-		Reflection *match;
-		double gradients[num_params];
-
-		/* If reflection is free-flagged, don't use it here */
-		if ( get_flag(refl) ) continue;
-
-		/* Find the full version */
-		get_indices(refl, &ha, &ka, &la);
-		match = find_refl(full, ha, ka, la);
-		if ( match == NULL ) continue;
-
-		/* Merged intensitty */
-		I_full = get_intensity(match);
-
-		/* Actual measurement of this reflection from this pattern */
-		I_partial = get_intensity(refl);
-		esd = get_esd_intensity(refl);
-		p = get_partiality(refl);
-
-		/* Scale only using strong reflections */
-		if ( I_partial <= 3.0*esd ) continue; /* Also because of log */
-		if ( get_redundancy(match) < 2 ) continue;
-		if ( I_full <= 0 ) continue;  /* Because log */
-		if ( p <= 0.0 ) continue; /* Because of log */
-
-		L = get_lorentz(refl);
-		s = resolution(crystal_get_cell(cr), ha, ka, la);
-
-		/* Calculate the weight for this reflection */
-		w = 1.0;
-
-		/* Calculate all gradients for this reflection */
-		for ( k=0; k<num_params; k++ ) {
-
-			if ( rv[k] == GPARAM_OSF ) {
-				gradients[k] = 1.0;
-			} else if ( rv[k] == GPARAM_BFAC ) {
-				gradients[k] = -s*s;
-			} else {
-				ERROR("Unrecognised scaling gradient.\n");
-				abort();
-			}
-		}
-
-		for ( k=0; k<num_params; k++ ) {
-
-			int g;
-			double v_c, v_curr;
-
-			for ( g=0; g<num_params; g++ ) {
-
-				double M_c, M_curr;
-
-				/* Matrix is symmetric */
-				if ( g > k ) continue;
-
-				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);
-				gsl_matrix_set(M, g, k, M_curr + M_c);
-
-			}
-
-			fx = G + log(p) - log(L) - B*s*s + log(I_full);
-			delta_I = log(I_partial) - fx;
-			v_c = w * delta_I * gradients[k];
-			v_curr = gsl_vector_get(v, k);
-			gsl_vector_set(v, k, v_curr + v_c);
-
-		}
-
-		nref++;
-	}
-
-	*nr = nref;
-
-	if ( nref < num_params ) {
-		crystal_set_user_flag(cr, PRFLAG_FEWREFL);
-		gsl_matrix_free(M);
-		gsl_vector_free(v);
-		return 0.0;
-	}
-
-	max_shift = 0.0;
-	shifts = solve_svd(v, M, NULL, 0);
-	if ( shifts != NULL ) {
-
-		for ( param=0; param<num_params; param++ ) {
-			double shift = gsl_vector_get(shifts, param);
-			apply_shift(cr, rv[param], shift);
-			if ( fabs(shift) > max_shift ) max_shift = fabs(shift);
-		}
-
-	} else {
-		crystal_set_user_flag(cr, PRFLAG_SOLVEFAIL);
-	}
-
-	gsl_matrix_free(M);
-	gsl_vector_free(v);
-	gsl_vector_free(shifts);
-
-	return max_shift;
-}
-
-
-double log_residual(Crystal *cr, const RefList *full, int free,
-                    int *pn_used, const char *filename)
-{
-	double dev = 0.0;
-	double G, B;
-	Reflection *refl;
-	RefListIterator *iter;
-	int n_used = 0;
-	FILE *fh = NULL;
-
-	G = crystal_get_osf(cr);
-	B = crystal_get_Bfac(cr);
-	if ( filename != NULL ) {
-		fh = fopen(filename, "a");
-		if ( fh == NULL ) {
-			ERROR("Failed to open '%s'\n", filename);
-		}
-	}
-
-	for ( refl = first_refl(crystal_get_reflections(cr), &iter);
-	      refl != NULL;
-	      refl = next_refl(refl, iter) )
-	{
-		double p, L, s, w;
-		signed int h, k, l;
-		Reflection *match;
-		double esd, I_full, I_partial;
-		double fx, dc;
-
-		if ( free && !get_flag(refl) ) continue;
-
-		get_indices(refl, &h, &k, &l);
-		match = find_refl(full, h, k, l);
-		if ( match == NULL ) continue;
-
-		p = get_partiality(refl);
-		L = get_lorentz(refl);
-		I_partial = get_intensity(refl);
-		I_full = get_intensity(match);
-		esd = get_esd_intensity(refl);
-		s = resolution(crystal_get_cell(cr), h, k, l);
-
-		if ( I_partial <= 3.0*esd ) continue; /* Also because of log */
-		if ( get_redundancy(match) < 2 ) continue;
-		if ( I_full <= 0 ) continue;  /* Because log */
-		if ( p <= 0.0 ) continue; /* Because of log */
-
-		fx = G + log(p) - log(L) - B*s*s + log(I_full);
-		dc = log(I_partial) - fx;
-		w = 1.0;
-		dev += w*dc*dc;
-
-		if ( fh != NULL ) {
-			fprintf(fh, "%4i %4i %4i %e %e\n",
-			        h, k, l, s, dev);
-		}
-
-	}
-
-	if ( fh != NULL ) fclose(fh);
-
-	if ( pn_used != NULL ) *pn_used = n_used;
-	return dev;
-}
-
-
-static void do_scale_refine(Crystal *cr, const RefList *full,
-                            PartialityModel pmodel, int *nr)
-{
-	int i, done;
-	double old_dev;
-
-	old_dev = log_residual(cr, full, 0, NULL, NULL);
-
-	i = 0;
-	done = 0;
-	do {
-
-		double dev;
-
-		scale_iterate(cr, full, pmodel, nr);
-
-		dev = log_residual(cr, full, 0, 0, NULL);
-		if ( fabs(dev - old_dev) < dev*0.01 ) done = 1;
-
-		i++;
-		old_dev = dev;
-
-	} while ( i < MAX_CYCLES && !done );
-}
+#include "scaling.h"
 
 
 struct scale_args
 {
 	RefList *full;
 	Crystal *crystal;
-	PartialityModel pmodel;
-	int n_reflections;
+	int flags;
 };
 
 
@@ -331,7 +63,6 @@ struct queue_args
 	int n_done;
 	Crystal **crystals;
 	int n_crystals;
-	long long int n_reflections;
 	struct scale_args task_defaults;
 };
 
@@ -339,8 +70,7 @@ struct queue_args
 static void scale_crystal(void *task, int id)
 {
 	struct scale_args *pargs = task;
-	do_scale_refine(pargs->crystal, pargs->full, pargs->pmodel,
-	                &pargs->n_reflections);
+	scale_one_crystal(pargs->crystal, pargs->full, pargs->flags);
 }
 
 
@@ -363,11 +93,7 @@ static void *get_crystal(void *vqargs)
 static void done_crystal(void *vqargs, void *task)
 {
 	struct queue_args *qa = vqargs;
-	struct scale_args *wargs = task;
-
 	qa->n_done++;
-	qa->n_reflections += wargs->n_reflections;
-
 	progress_bar(qa->n_done, qa->n_crystals, "Scaling");
 	free(task);
 }
@@ -381,14 +107,12 @@ static double total_log_r(Crystal **crystals, int n_crystals, RefList *full,
 	int n = 0;
 
 	for ( i=0; i<n_crystals; i++ ) {
-
 		double r;
 		if ( crystal_get_user_flag(crystals[i]) ) continue;
 		r = log_residual(crystals[i], full, 0, NULL, NULL);
 		if ( isnan(r) ) continue;
 		total += r;
 		n++;
-
 	}
 
 	if ( ninc != NULL ) *ninc = n;
@@ -397,8 +121,7 @@ static double total_log_r(Crystal **crystals, int n_crystals, RefList *full,
 
 
 /* Perform iterative scaling, all the way to convergence */
-void scale_all(Crystal **crystals, int n_crystals, int nthreads,
-               PartialityModel pmodel)
+void scale_all(Crystal **crystals, int n_crystals, int nthreads, int scaleflags)
 {
 	struct scale_args task_defaults;
 	struct queue_args qargs;
@@ -406,7 +129,7 @@ void scale_all(Crystal **crystals, int n_crystals, int nthreads,
 	int niter = 0;
 
 	task_defaults.crystal = NULL;
-	task_defaults.pmodel = pmodel;
+	task_defaults.flags = scaleflags;
 
 	qargs.task_defaults = task_defaults;
 	qargs.n_crystals = n_crystals;
@@ -422,18 +145,15 @@ void scale_all(Crystal **crystals, int n_crystals, int nthreads,
 		double bef_res;
 
 		full = merge_intensities(crystals, n_crystals, nthreads,
-		                         pmodel, 2, INFINITY, 0);
+		                         2, INFINITY, 0);
 		old_res = new_res;
 		bef_res = total_log_r(crystals, n_crystals, full, NULL);
 
 		qargs.task_defaults.full = full;
 		qargs.n_started = 0;
 		qargs.n_done = 0;
-		qargs.n_reflections = 0;
 		run_threads(nthreads, scale_crystal, get_crystal, done_crystal,
 		            &qargs, n_crystals, 0, 0, 0);
-		STATUS("%lli reflections went into the scaling.\n",
-		       qargs.n_reflections);
 
 		new_res = total_log_r(crystals, n_crystals, full, &ninc);
 		STATUS("Log residual went from %e to %e, %i crystals\n",
@@ -456,3 +176,164 @@ void scale_all(Crystal **crystals, int n_crystals, int nthreads,
 		ERROR("Too many iterations - giving up!\n");
 	}
 }
+
+
+/* Calculates G and B, by which cr's reflections should be multiplied to fit reference */
+int scale_one_crystal(Crystal *cr, const RefList *listR, int flags)
+{
+	const Reflection *reflS;
+	RefListIterator *iter;
+	int max_n = 256;
+	int n = 0;
+	double *x;
+	double *y;
+	double *w;
+	int r;
+	double cov00, cov01, cov11, chisq;
+	int n_esdS = 0;
+	int n_esdR = 0;
+	int n_ihS = 0;
+	int n_ihR = 0;
+	int n_nanS = 0;
+	int n_nanR = 0;
+	int n_infS = 0;
+	int n_infR = 0;
+	int n_part = 0;
+	int n_nom = 0;
+	int n_red = 0;
+	RefList *listS = crystal_get_reflections(cr);
+	UnitCell *cell = crystal_get_cell(cr);
+	double G, B;
+
+	assert(cell != NULL);
+	assert(listR != NULL);
+	assert(listS != NULL);
+
+	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 ( reflS = first_refl_const(listS, &iter);
+	      reflS != NULL;
+	      reflS = next_refl_const(reflS, iter) )
+	{
+		signed int h, k, l;
+		const Reflection *reflR;
+		double IhR, IhS, esdS, pS, LS;
+		double s;
+
+		nb++;
+
+		get_indices(reflS, &h, &k, &l);
+		reflR = find_refl(listR, h, k, l);
+		if ( reflR == NULL ) {
+			n_nom++;
+			continue;
+		}
+
+		s = resolution(cell, h, k, l);
+
+		IhR = get_intensity(reflR);
+		IhS = get_intensity(reflS);
+		esdS = get_esd_intensity(reflS);
+		pS = get_partiality(reflS);
+		LS = get_lorentz(reflS);
+
+		/* Problem cases in approximate descending order of severity */
+		if ( isnan(IhR) ) { n_nanR++; continue; }
+		if ( isinf(IhR) ) { n_infR++; continue; }
+		if ( isnan(IhS) ) { n_nanS++; continue; }
+		if ( isinf(IhS) ) { n_infS++; continue; }
+		if ( pS < 0.3 ) { n_part++; continue; }
+		if ( IhS <= 0.0 ) { n_ihS++; continue; }
+		if ( IhS <= 3.0*esdS ) { n_esdS++; continue; }
+		if ( IhR <= 0.0 ) { n_ihR++; continue; }
+		if ( get_redundancy(reflR) < 2 ) { n_red++; 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] = s*s;
+		y[n] = log(LS) + log(IhS) -log(pS) - log(IhR);
+		w[n] = pS;
+		n++;
+
+	}
+
+	if ( n < 2 ) {
+		if ( flags & SCALE_VERBOSE_ERRORS ) {
+			ERROR("Not enough reflections for scaling (had %i, but %i remain)\n", nb, n);
+			if ( n_esdR ) ERROR("%i reference reflection esd\n", n_esdR);
+			if ( n_esdS ) ERROR("%i subject reflection esd\n", n_esdS);
+			if ( n_ihR ) ERROR("%i reference reflection intensity\n", n_ihR);
+			if ( n_red ) ERROR("%i reference reflection redundancy\n", n_red);
+			if ( n_ihS ) ERROR("%i subject reflection intensity\n", n_ihS);
+			if ( n_nanR ) ERROR("%i reference reflection nan\n", n_nanR);
+			if ( n_nanS ) ERROR("%i subject reflection nan\n", n_nanS);
+			if ( n_infR ) ERROR("%i reference reflection inf\n", n_infR);
+			if ( n_infS ) ERROR("%i subject reflection inf\n", n_infS);
+			if ( n_part ) ERROR("%i subject reflection partiality\n", n_part);
+			if ( n_nom ) ERROR("%i no match in reference list\n", n_nom);
+		}
+		free(x);
+		free(y);
+		free(w);
+		return 1;
+	}
+
+	if ( flags & SCALE_NO_B ) {
+		G = gsl_stats_wmean(w, 1, y, 1, n);
+		B = 0.0;
+		r = 0;
+	} else {
+		r = gsl_fit_wlinear(x, 1, w, 1, y, 1, n, &G, &B, &cov00, &cov01, &cov11, &chisq);
+	}
+
+	if ( r ) {
+		ERROR("Scaling failed.\n");
+		free(x);
+		free(y);
+		free(w);
+		return 1;
+	}
+
+	if ( isnan(G) ) {
+
+		if ( flags & SCALE_VERBOSE_ERRORS ) {
+			ERROR("Scaling gave NaN (%i pairs)\n", n);
+			if ( n < 10 ) {
+				int i;
+				for ( i=0; i<n; i++ ) {
+					STATUS("%3i %e %e %e\n", i, x[i], y[i], w[i]);
+				}
+			}
+		}
+
+		free(x);
+		free(y);
+		free(w);
+		return 1;
+	}
+
+	crystal_set_osf(cr, exp(G));
+	crystal_set_Bfac(cr, -B);
+
+	free(x);
+	free(y);
+	free(w);
+
+	return 0;
+}