Switch to simpler scaling algorithm

4 files changed, 61 insertions, 367 deletions

diff --git a/src/partialator.c b/src/partialator.c
index b5f24b1d..7c3a481c 100644
--- a/src/partialator.c
+++ b/src/partialator.c
@@ -1221,22 +1221,21 @@ int main(int argc, char *argv[])
 	STATUS("Checking patterns.\n");
 	//early_rejection(crystals, n_crystals);
 
-	/* Initial rejection, figures of merit etc */
+	/* Create reference data set if we don't already have one */
 	if ( reference == NULL ) {
-		if ( !no_scale ) {
-			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);
 
+	/* Scale everything to the reference */
 	if ( !no_scale ) {
-		scale_all_to_reference(crystals, n_crystals, full);
+		scale_all_to_reference(crystals, n_crystals, full, nthreads);
 	}
 
+	/* Check rejection and write figures of merit */
+	check_rejection(crystals, n_crystals, full, max_B);
 	show_all_residuals(crystals, n_crystals, full);
 	write_pgraph(full, crystals, n_crystals, 0);
 	write_specgraph(full, crystals[0], 0);
@@ -1248,41 +1247,19 @@ int main(int argc, char *argv[])
 
 		if ( !no_pr ) {
 			refine_all(crystals, n_crystals, full, nthreads, pmodel);
+		} else if ( !no_scale ) {
+			scale_all_to_reference(crystals, n_crystals, full, nthreads);
 		}
 
-		if ( !no_scale ) {
-
-			if ( reference == NULL ) {
-
-				/* No reference -> XSCALE-like algorithm */
-				full = merge_intensities(crystals, n_crystals,
-				                         nthreads, pmodel,
-				                         min_measurements,
-				                         push_res, 1);
-
-			} else {
-
-				/* Reference -> Ginn-style linear algorithm */
-				if ( !no_scale ) {
-					scale_all_to_reference(crystals, n_crystals,
-					                       reference);
-				}
-
-			}
-		}
-
-		check_rejection(crystals, n_crystals, full, max_B);
-
+		/* Create new reference if needed */
 		if ( reference == NULL ) {
-			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);
 		} /* else full still equals reference */
 
+		check_rejection(crystals, n_crystals, full, max_B);
 		show_all_residuals(crystals, n_crystals, full);
 		write_pgraph(full, crystals, n_crystals, i+1);
 		write_specgraph(full, crystals[0], i+1);
@@ -1319,34 +1296,24 @@ int main(int argc, char *argv[])
 		}
 	}
 
-
+	/* Final merge */
+	STATUS("Final merge...\n");
 	if ( reference == NULL ) {
-		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);
-	} /* else full still equals reference */
-	if ( !no_scale ) {
-		scale_all_to_reference(crystals, n_crystals, full);
+	} else {
+		scale_all_to_reference(crystals, n_crystals, reference, nthreads);
+		full = merge_intensities(crystals, n_crystals, nthreads,
+		                         pmodel, min_measurements, push_res, 1);
 	}
+
+	/* Write final figures of merit (no rejection any more) */
 	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 */
+	STATUS("Final profile radius: %e\n", crystal_get_profile_radius(crystals[0]));
 
 	/* Output results */
 	STATUS("Writing overall results to %s\n", outfile);
diff --git a/src/post-refinement.c b/src/post-refinement.c
index af701575..7402b91c 100644
--- a/src/post-refinement.c
+++ b/src/post-refinement.c
@@ -381,6 +381,8 @@ static void do_pr_refine(Crystal *cr, const RefList *full,
 	int n_params = 0;
 	int n_iter = 0;
 	int status;
+	int r;
+	double G;
 	double residual_start, residual_free_start;
 
 	residual_start = residual(cr, full, 0, NULL, NULL);
@@ -456,6 +458,12 @@ static void do_pr_refine(Crystal *cr, const RefList *full,
 	apply_parameters(min->x, initial, rv, cr);
 	update_predictions(cr);
 	calculate_partialities(cr, PMODEL_XSPHERE);
+	r = linear_scale(full, crystal_get_reflections(cr), &G);
+	if ( r == 0 ) {
+		crystal_set_osf(cr, G);
+	} else {
+		fprintf(stderr, "Scaling failure after refinement.\n");
+	}
 
 	if ( verbose ) {
 		STATUS("PR final: dev = %10.5e, free dev = %10.5e\n",
diff --git a/src/scaling.c b/src/scaling.c
index b2f063ec..bf8857ac 100644
--- a/src/scaling.c
+++ b/src/scaling.c
@@ -48,184 +48,6 @@
 #include "scaling.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 = -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,
-                            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 = -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)
 {
@@ -292,32 +114,6 @@ double log_residual(Crystal *cr, const RefList *full, int free,
 }
 
 
-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 );
-}
-
-
 struct scale_args
 {
 	RefList *full;
@@ -341,8 +137,14 @@ 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);
+	int r;
+	double G;
+
+	/* Simple iterative algorithm */
+	r = linear_scale(pargs->full, crystal_get_reflections(pargs->crystal), &G);
+	if ( r == 0 ) {
+		crystal_set_osf(pargs->crystal, G);
+	} /* else don't change it */
 }
 
 
@@ -375,91 +177,6 @@ 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;
-	task_defaults.pmodel = pmodel;
-
-	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");
-	}
-}
-
-
 /* Calculates G, by which list2 should be multiplied to fit list1 */
 int linear_scale(const RefList *list1, const RefList *list2, double *G)
 {
@@ -490,6 +207,7 @@ int linear_scale(const RefList *list1, const RefList *list2, double *G)
 	{
 		signed int h, k, l;
 		double Ih1, Ih2;
+		double esd1, esd2;
 		nb++;
 
 		get_indices(refl1, &h, &k, &l);
@@ -500,8 +218,13 @@ int linear_scale(const RefList *list1, const RefList *list2, double *G)
 
 		Ih1 = get_intensity(refl1);
 		Ih2 = get_intensity(refl2);
+		esd1 = get_esd_intensity(refl1);
+		esd2 = get_esd_intensity(refl2);
 
+		if ( Ih1 <= 3.0*esd1 ) continue;
+		if ( Ih2 <= 3.0*esd2 ) continue;
 		if ( (Ih1 <= 0.0) || (Ih2 <= 0.0) ) continue;
+		if ( get_redundancy(refl1) < 2 ) continue;
 		if ( isnan(Ih1) || isinf(Ih1) ) continue;
 		if ( isnan(Ih2) || isinf(Ih2) ) continue;
 
@@ -537,7 +260,7 @@ int linear_scale(const RefList *list1, const RefList *list2, double *G)
 
 	if ( isnan(*G) ) {
 		ERROR("Scaling gave NaN (%i pairs)\n", n);
-		abort();
+		*G = 1.0;
 		return 1;
 	}
 
@@ -550,25 +273,24 @@ int linear_scale(const RefList *list1, const RefList *list2, double *G)
 
 
 void scale_all_to_reference(Crystal **crystals, int n_crystals,
-                            RefList *reference)
+                            RefList *reference, int nthreads)
 {
-	int i;
-
-	for ( i=0; i<n_crystals; i++ ) {
-		double G;
-		if ( linear_scale(reference,
-		                  crystal_get_reflections(crystals[i]),
-		                  &G) == 0 )
-		{
-			if ( isnan(G) ) {
-				ERROR("NaN scaling factor for crystal %i\n", i);
-			}
-			crystal_set_osf(crystals[i], G);
-			crystal_set_Bfac(crystals[i], 0.0);
-		} else {
-			ERROR("Scaling failed for crystal %i\n", i);
-		}
-		progress_bar(i, n_crystals, "Scaling to reference");
-	}
-	progress_bar(n_crystals, n_crystals, "Scaling to reference");
+	struct scale_args task_defaults;
+	struct queue_args qargs;
+
+	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;
+
+	qargs.task_defaults.full = reference;
+	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);
 }
diff --git a/src/scaling.h b/src/scaling.h
index 4d075a3b..6ff0a337 100644
--- a/src/scaling.h
+++ b/src/scaling.h
@@ -43,10 +43,7 @@ extern double log_residual(Crystal *cr, const RefList *full, int free,
 
 extern int linear_scale(const RefList *list1, const RefList *list2, double *G);
 
-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);
+                                   RefList *reference, int nthreads);
 
 #endif	/* SCALING_H */