Restore old scaling code

2 files changed, 301 insertions, 12 deletions

diff --git a/src/scaling.c b/src/scaling.c
index 4f92f501..0ad50634 100644
--- a/src/scaling.c
+++ b/src/scaling.c
@@ -48,6 +48,179 @@
 #include "scaling.h"
 
 
+/* 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 = -log(crystal_get_osf(cr));
+		t += shift;
+		crystal_set_osf(cr, exp(-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, 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 = -log(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)
 {
@@ -114,11 +287,35 @@ double log_residual(Crystal *cr, const RefList *full, int free,
 }
 
 
+static void do_scale_refine(Crystal *cr, const RefList *full, 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, 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 < 10 && !done );
+}
+
+
 struct scale_args
 {
 	RefList *full;
 	Crystal *crystal;
-	PartialityModel pmodel;
 	int n_reflections;
 };
 
@@ -137,16 +334,7 @@ struct queue_args
 static void scale_crystal(void *task, int id)
 {
 	struct scale_args *pargs = task;
-	int r;
-	double G;
-
-	/* Simple iterative algorithm */
-	r = linear_scale(pargs->full, crystal_get_reflections(pargs->crystal), &G, 0);
-	if ( r == 0 ) {
-		crystal_set_osf(pargs->crystal, G);
-	} else {
-		crystal_set_user_flag(pargs->crystal, PRFLAG_SCALEBAD);
-	}
+	do_scale_refine(pargs->crystal, pargs->full, &pargs->n_reflections);
 }
 
 
@@ -179,6 +367,104 @@ static void done_crystal(void *vqargs, void *task)
 }
 
 
+static double total_log_r(Crystal **crystals, int n_crystals, RefList *full,
+                          int *ninc)
+{
+	int i;
+	double total = 0.0;
+	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;
+	return total;
+}
+
+
+/* Perform iterative scaling, all the way to convergence */
+void scale_all(Crystal **crystals, int n_crystals, int nthreads,
+               PartialityModel pmodel)
+{
+	struct scale_args task_defaults;
+	struct queue_args qargs;
+	double old_res, new_res;
+	int niter = 0;
+
+	task_defaults.crystal = NULL;
+
+	qargs.task_defaults = task_defaults;
+	qargs.n_crystals = n_crystals;
+	qargs.crystals = crystals;
+
+	/* Don't have threads which are doing nothing */
+	if ( n_crystals < nthreads ) nthreads = n_crystals;
+
+	new_res = INFINITY;
+	do {
+		RefList *full;
+		int ninc;
+		double bef_res;
+
+		full = merge_intensities(crystals, n_crystals, nthreads,
+		                         pmodel, 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",
+		       bef_res, new_res, ninc);
+
+		int i;
+		double meanB = 0.0;
+		for ( i=0; i<n_crystals; i++ ) {
+			meanB += crystal_get_Bfac(crystals[i]);
+		}
+		meanB /= n_crystals;
+		STATUS("Mean B = %e\n", meanB);
+
+		reflist_free(full);
+		niter++;
+
+	} while ( (fabs(new_res-old_res) >= 0.01*old_res) && (niter < 10) );
+
+	if ( niter == 10 ) {
+		ERROR("Too many iterations - giving up!\n");
+	}
+}
+
+
+static void scale_crystal_linear(void *task, int id)
+{
+	struct scale_args *pargs = task;
+	int r;
+	double G;
+
+	/* Simple iterative algorithm */
+	r = linear_scale(pargs->full, crystal_get_reflections(pargs->crystal), &G, 0);
+	if ( r == 0 ) {
+		crystal_set_osf(pargs->crystal, G);
+	} else {
+		crystal_set_user_flag(pargs->crystal, PRFLAG_SCALEBAD);
+	}
+}
+
+
 /* Calculates G, by which list2 should be multiplied to fit list1 */
 int linear_scale(const RefList *list1, const RefList *list2, double *G,
                  int complain_loudly)
@@ -333,6 +619,6 @@ void scale_all_to_reference(Crystal **crystals, int n_crystals,
 	qargs.n_started = 0;
 	qargs.n_done = 0;
 	qargs.n_reflections = 0;
-	run_threads(nthreads, scale_crystal, get_crystal, done_crystal,
+	run_threads(nthreads, scale_crystal_linear, get_crystal, done_crystal,
 	            &qargs, n_crystals, 0, 0, 0);
 }
diff --git a/src/scaling.h b/src/scaling.h
index b345029a..620518ef 100644
--- a/src/scaling.h
+++ b/src/scaling.h
@@ -44,6 +44,9 @@ extern double log_residual(Crystal *cr, const RefList *full, int free,
 extern int linear_scale(const RefList *list1, const RefList *list2, double *G,
                         int complain_loudly);
 
+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, int nthreads);