Preparation for adjusting B factor during post-refinement

Add --no-Bscale option to partialator, and pass down as far as needed
residual() no longer does scaling: call scale_one_crystal() first if necessary
scale_one() replaces old linear_scale() function to scale a pair of RefLists
 (but so far does the same as the old function)

6 files changed, 99 insertions, 50 deletions

diff --git a/src/partialator.c b/src/partialator.c
index 0f3bd104..e775caef 100644
--- a/src/partialator.c
+++ b/src/partialator.c
@@ -317,6 +317,7 @@ static void show_help(const char *s)
 "      --stop-after=<n>       Stop after merging <n> crystals.\n"
 "  -n, --iterations=<n>       Run <n> cycles of scaling and post-refinement.\n"
 "      --no-scale             Disable scale factor (G, B) refinement.\n"
+"      --no-Bscale            Disable B factor scaling.\n"
 "      --no-pr                Disable orientation/physics refinement.\n"
 "  -m, --model=<model>        Specify partiality model.\n"
 "      --min-measurements=<n> Minimum number of measurements to require.\n"
@@ -717,6 +718,7 @@ static void all_residuals(Crystal **crystals, int n_crystals, RefList *full,
 
 		if ( crystal_get_user_flag(crystals[i]) ) continue;
 
+		/* Scaling should have been done right before calling this */
 		r = residual(crystals[i], full, 0, NULL, NULL);
 		free_r = residual(crystals[i], full, 1, NULL, NULL);
 		log_r = log_residual(crystals[i], full, 0, NULL, NULL);
@@ -754,6 +756,7 @@ struct log_qargs
 	Crystal **crystals;
 	int n_crystals;
 	RefList *full;
+	int scaleflags;
 };
 
 
@@ -761,6 +764,7 @@ struct log_args
 {
 	Crystal *cr;
 	RefList *full;
+	int scaleflags;
 	int iter;
 	int cnum;
 };
@@ -780,6 +784,7 @@ static void *get_log_task(void *vp)
 	task->full = qargs->full;
 	task->iter = qargs->iter;
 	task->cnum = qargs->next;
+	task->scaleflags = qargs->scaleflags;
 
 	qargs->next += 20;
 	return task;
@@ -790,7 +795,8 @@ static void write_logs(void *vp, int cookie)
 {
 	struct log_args *args = vp;
 	write_specgraph(args->cr, args->full, args->iter, args->cnum);
-	write_gridscan(args->cr, args->full, args->iter, args->cnum);
+	write_gridscan(args->cr, args->full, args->iter, args->cnum,
+	               args->scaleflags);
 	write_test_logs(args->cr, args->full, args->iter, args->cnum);
 }
 
@@ -803,16 +809,17 @@ static void done_log(void *qargs, void *vp)
 
 
 static void write_logs_parallel(Crystal **crystals, int n_crystals,
-                                RefList *full, int iter, int n_threads)
+                                RefList *full, int iter, int n_threads,
+                                int scaleflags)
 {
 	struct log_qargs qargs;
 
-
 	qargs.iter = iter;
 	qargs.next = 0;
 	qargs.full = full;
 	qargs.crystals = crystals;
 	qargs.n_crystals = n_crystals;
+	qargs.scaleflags = scaleflags;
 
 	STATUS("Writing logs...\n");
 	run_threads(n_threads, write_logs, get_log_task, done_log, &qargs,
@@ -839,6 +846,7 @@ int main(int argc, char *argv[])
 	int n_crystals_seen = 0;
 	char cmdline[1024];
 	int no_scale = 0;
+	int no_Bscale = 0;
 	int no_pr = 0;
 	Stream *st;
 	Crystal **crystals;
@@ -867,6 +875,7 @@ int main(int argc, char *argv[])
 	double force_bandwidth = -1.0;
 	double force_radius = -1.0;
 	char *audit_info;
+	int scaleflags = 0;
 
 	/* Long options */
 	const struct option longopts[] = {
@@ -895,6 +904,7 @@ int main(int argc, char *argv[])
 		{"force-radius",       1, NULL,               11},
 
 		{"no-scale",           0, &no_scale,           1},
+		{"no-Bscale",          0, &no_Bscale,          1},
 		{"no-pr",              0, &no_pr,              1},
 		{"no-polarisation",    0, &polarisation,       0},
 		{"no-polarization",    0, &polarisation,       0}, /* compat */
@@ -1146,6 +1156,10 @@ int main(int argc, char *argv[])
 		       "partialities (--model=unity).\n");
 	}
 
+	if ( no_Bscale ) {
+		scaleflags |= SCALE_NO_B;
+	}
+
 	if ( !no_logs ) {
 		int r = mkdir("pr-logs", S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH);
 		if ( r ) {
@@ -1327,7 +1341,7 @@ int main(int argc, char *argv[])
 	/* Create reference data set if we don't already have one */
 	if ( reference == NULL ) {
 		STATUS("Initial scaling...\n");
-		scale_all(crystals, n_crystals, nthreads);
+		scale_all(crystals, n_crystals, nthreads, scaleflags);
 		full = merge_intensities(crystals, n_crystals, nthreads,
 		                         min_measurements, push_res, 1);
 	} else {
@@ -1339,7 +1353,8 @@ int main(int argc, char *argv[])
 	show_all_residuals(crystals, n_crystals, full);
 	if ( !no_pr && !no_logs ) {
 		write_pgraph(full, crystals, n_crystals, 0, "");
-		write_logs_parallel(crystals, n_crystals, full, 0, nthreads);
+		write_logs_parallel(crystals, n_crystals, full, 0, nthreads,
+		                    scaleflags);
 	}
 
 	/* Iterate */
@@ -1349,14 +1364,15 @@ int main(int argc, char *argv[])
 
 		if ( !no_pr ) {
 			refine_all(crystals, n_crystals, full, nthreads, pmodel,
-			           0, i+1, no_logs, sym, amb);
+			           0, i+1, no_logs, sym, amb, scaleflags);
 		}
 
 		/* Create new reference if needed */
 		if ( reference == NULL ) {
 			reflist_free(full);
 			if ( !no_scale ) {
-				scale_all(crystals, n_crystals, nthreads);
+				scale_all(crystals, n_crystals, nthreads,
+				          scaleflags);
 			}
 			full = merge_intensities(crystals, n_crystals, nthreads,
 			                         min_measurements,
@@ -1403,7 +1419,7 @@ int main(int argc, char *argv[])
 	if ( reference == NULL ) {
 		reflist_free(full);
 		if ( !no_scale ) {
-			scale_all(crystals, n_crystals, nthreads);
+			scale_all(crystals, n_crystals, nthreads, scaleflags);
 		}
 		full = merge_intensities(crystals, n_crystals, nthreads,
 		                         min_measurements,
@@ -1417,7 +1433,8 @@ int main(int argc, char *argv[])
 	show_all_residuals(crystals, n_crystals, full);
 	if ( !no_pr && !no_logs ) {
 		write_pgraph(full, crystals, n_crystals, -1, "");
-		write_logs_parallel(crystals, n_crystals, full, -1, nthreads);
+		write_logs_parallel(crystals, n_crystals, full, -1, nthreads,
+		                    scaleflags);
 	}
 
 	/* Output results */
diff --git a/src/post-refinement.c b/src/post-refinement.c
index f260da7f..293fe332 100644
--- a/src/post-refinement.c
+++ b/src/post-refinement.c
@@ -95,10 +95,6 @@ double residual(Crystal *cr, const RefList *full, int free,
 	double B = crystal_get_Bfac(cr);
 	UnitCell *cell = crystal_get_cell(cr);
 
-	if ( linear_scale(full, crystal_get_reflections(cr), &G, complain) ) {
-		return GSL_NAN;
-	}
-
 	for ( refl = first_refl(crystal_get_reflections(cr), &iter);
 	      refl != NULL;
 	      refl = next_refl(refl, iter) )
@@ -232,6 +228,7 @@ struct rf_priv {
 	int verbose;
 	int serial;
 	gsl_vector *initial;
+	int scaleflags;
 
 	/* For freeing later */
 	gsl_vector *vals;
@@ -349,6 +346,11 @@ static double residual_f(const gsl_vector *v, void *pp)
 	update_predictions(cr);
 	calculate_partialities(cr, PMODEL_XSPHERE);
 
+	if ( scale_one_crystal(cr, pv->full, pv->scaleflags) ) {
+		crystal_free(cr);
+		if ( pv->verbose ) STATUS("Bad scaling\n");
+		return GSL_NAN;
+	}
 	res = residual(cr, pv->full, 0, NULL, NULL);
 
 	cell_free(crystal_get_cell(cr));
@@ -463,8 +465,8 @@ static void reindex_cell(UnitCell *cell, SymOpList *amb, int idx)
 }
 
 
-void try_reindex(Crystal *crin, const RefList *full,
-                 SymOpList *sym, SymOpList *amb)
+static void try_reindex(Crystal *crin, const RefList *full,
+                        SymOpList *sym, SymOpList *amb, int scaleflags)
 {
 	RefList *list;
 	Crystal *cr;
@@ -474,6 +476,7 @@ void try_reindex(Crystal *crin, const RefList *full,
 
 	if ( sym == NULL || amb == NULL ) return;
 
+	if ( scale_one_crystal(crin, full, scaleflags) ) return;
 	residual_original = residual(crin, full, 0, NULL, NULL);
 
 	cr = crystal_copy(crin);
@@ -494,6 +497,7 @@ void try_reindex(Crystal *crin, const RefList *full,
 		update_predictions(cr);
 		calculate_partialities(cr, PMODEL_XSPHERE);
 
+		if ( scale_one_crystal(cr, full, scaleflags) ) return;
 		residual_flipped = residual(cr, full, 0, NULL, NULL);
 
 		if ( residual_flipped < residual_original ) {
@@ -637,6 +641,7 @@ void write_specgraph(Crystal *crystal, const RefList *full,
 
 static gsl_multimin_fminimizer *setup_minimiser(Crystal *cr, const RefList *full,
                                                 int verbose, int serial,
+                                                int scaleflags,
                                                 struct rf_priv *priv)
 {
 	gsl_multimin_fminimizer *min;
@@ -654,6 +659,7 @@ static gsl_multimin_fminimizer *setup_minimiser(Crystal *cr, const RefList *full
 	priv->full = full;
 	priv->verbose = verbose;
 	priv->serial = serial;
+	priv->scaleflags = scaleflags;
 
 	priv->f.f = residual_f;
 	priv->f.n = n_params;
@@ -678,7 +684,7 @@ static gsl_multimin_fminimizer *setup_minimiser(Crystal *cr, const RefList *full
 
 
 static void write_grid(Crystal *cr, const RefList *full,
-                       signed int cycle, int serial,
+                       signed int cycle, int serial, int scaleflags,
                        const enum gparam par1, const enum gparam par2,
                        const char *name)
 {
@@ -690,7 +696,7 @@ static void write_grid(Crystal *cr, const RefList *full,
 	int idx1, idx2;
 	int i;
 
-	min = setup_minimiser(cr, full, 0, serial, &priv);
+	min = setup_minimiser(cr, full, 0, serial, scaleflags, &priv);
 
 	idx1 = 99;
 	idx2 = 99;
@@ -747,13 +753,13 @@ static void write_grid(Crystal *cr, const RefList *full,
 
 
 void write_gridscan(Crystal *cr, const RefList *full,
-                    signed int cycle, int serial)
+                    signed int cycle, int serial, int scaleflags)
 {
-	write_grid(cr, full, cycle, serial,
+	write_grid(cr, full, cycle, serial, scaleflags,
 	           GPARAM_ANG1, GPARAM_ANG2, "ang1-ang2");
-	write_grid(cr, full, cycle, serial,
+	write_grid(cr, full, cycle, serial, scaleflags,
 	           GPARAM_ANG1, GPARAM_WAVELENGTH, "ang1-wave");
-	write_grid(cr, full, cycle, serial,
+	write_grid(cr, full, cycle, serial, scaleflags,
 	           GPARAM_R, GPARAM_WAVELENGTH, "R-wave");
 }
 
@@ -761,19 +767,21 @@ void write_gridscan(Crystal *cr, const RefList *full,
 static void do_pr_refine(Crystal *cr, const RefList *full,
                          PartialityModel pmodel, int verbose, int serial,
                          int cycle, int write_logs,
-                         SymOpList *sym, SymOpList *amb)
+                         SymOpList *sym, SymOpList *amb, int scaleflags)
 {
 	gsl_multimin_fminimizer *min;
 	struct rf_priv priv;
 	int n_iter = 0;
 	int status;
-	int r;
-	double G;
 	double residual_start, residual_free_start;
 	FILE *fh = NULL;
 
-	try_reindex(cr, full, sym, amb);
+	try_reindex(cr, full, sym, amb, scaleflags);
 
+	if ( scale_one_crystal(cr, full, scaleflags) ) {
+		ERROR("Bad scaling at start of refinement.\n");
+		return;
+	}
 	residual_start = residual(cr, full, 0, NULL, NULL);
 	residual_free_start = residual(cr, full, 1, NULL, NULL);
 
@@ -782,7 +790,7 @@ static void do_pr_refine(Crystal *cr, const RefList *full,
 		       residual_start, residual_free_start);
 	}
 
-	min = setup_minimiser(cr, full, verbose, serial, &priv);
+	min = setup_minimiser(cr, full, verbose, serial, scaleflags, &priv);
 
 	if ( verbose ) {
 		double res = residual_f(min->x, &priv);
@@ -902,10 +910,7 @@ static void do_pr_refine(Crystal *cr, const RefList *full,
 		apply_parameters(min->x, priv.initial, priv.rv, cr);
 		update_predictions(cr);
 		calculate_partialities(cr, PMODEL_XSPHERE);
-		r = linear_scale(full, crystal_get_reflections(cr), &G, 0);
-		if ( r == 0 ) {
-			crystal_set_osf(cr, G);
-		}
+		scale_one_crystal(cr, full, scaleflags);
 
 		if ( verbose ) {
 
@@ -922,7 +927,7 @@ static void do_pr_refine(Crystal *cr, const RefList *full,
 	}
 
 	if ( write_logs ) {
-		write_gridscan(cr, full, cycle, serial);
+		write_gridscan(cr, full, cycle, serial, scaleflags);
 		write_specgraph(cr, full, cycle, serial);
 		write_test_logs(cr, full, cycle, serial);
 	}
@@ -950,6 +955,7 @@ struct refine_args
 	int no_logs;
 	SymOpList *sym;
 	SymOpList *amb;
+	int scaleflags;
 };
 
 
@@ -974,7 +980,7 @@ static void refine_image(void *task, int id)
 
 	do_pr_refine(cr, pargs->full, pargs->pmodel, pargs->verbose,
 	             pargs->serial, pargs->cycle, write_logs,
-	             pargs->sym, pargs->amb);
+	             pargs->sym, pargs->amb, pargs->scaleflags);
 
 	if ( crystal_get_user_flag(cr) == 0 ) {
 		pargs->prdata.refined = 1;
@@ -1013,7 +1019,7 @@ static void done_image(void *vqargs, void *task)
 void refine_all(Crystal **crystals, int n_crystals,
                 RefList *full, int nthreads, PartialityModel pmodel,
                 int verbose, int cycle, int no_logs,
-                SymOpList *sym, SymOpList *amb)
+                SymOpList *sym, SymOpList *amb, int scaleflags)
 {
 	struct refine_args task_defaults;
 	struct queue_args qargs;
@@ -1027,6 +1033,7 @@ void refine_all(Crystal **crystals, int n_crystals,
 	task_defaults.no_logs = no_logs;
 	task_defaults.sym = sym;
 	task_defaults.amb = amb;
+	task_defaults.scaleflags = scaleflags;
 
 	qargs.task_defaults = task_defaults;
 	qargs.n_started = 0;
diff --git a/src/post-refinement.h b/src/post-refinement.h
index 8c345729..b8923d2c 100644
--- a/src/post-refinement.h
+++ b/src/post-refinement.h
@@ -61,10 +61,10 @@ extern const char *str_prflag(enum prflag flag);
 extern void refine_all(Crystal **crystals, int n_crystals,
                        RefList *full, int nthreads, PartialityModel pmodel,
                        int verbose, int cycle, int no_logs,
-                       SymOpList *sym, SymOpList *amb);
+                       SymOpList *sym, SymOpList *amb, int scaleflags);
 
 extern void write_gridscan(Crystal *cr, const RefList *full,
-                           int cycle, int serial);
+                           int cycle, int serial, int scaleflags);
 
 extern void write_specgraph(Crystal *crystal, const RefList *full,
                             signed int cycle, int serial);
diff --git a/src/scaling.c b/src/scaling.c
index 39541926..5f896e65 100644
--- a/src/scaling.c
+++ b/src/scaling.c
@@ -389,7 +389,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)
+void scale_all(Crystal **crystals, int n_crystals, int nthreads, int no_Bscale)
 {
 	struct scale_args task_defaults;
 	struct queue_args qargs;
@@ -448,10 +448,11 @@ void scale_all(Crystal **crystals, int n_crystals, int nthreads)
 }
 
 
-/* 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)
+/* Calculates G and B, by which list2 should be multiplied to fit list1 */
+int scale_one(const RefList *list1, const RefList *list2, int flags,
+              double *G, double *B)
 {
+	int complain_loudly = 0;
 	const Reflection *refl1;
 	const Reflection *refl2;
 	RefListIterator *iter;
@@ -474,6 +475,8 @@ int linear_scale(const RefList *list1, const RefList *list2, double *G,
 	int n_part = 0;
 	int n_nom = 0;
 
+	*B = 0.0;  /* FIXME */
+
 	x = malloc(max_n*sizeof(double));
 	w = malloc(max_n*sizeof(double));
 	y = malloc(max_n*sizeof(double));
@@ -556,6 +559,7 @@ int linear_scale(const RefList *list1, const RefList *list2, double *G,
 
 	if ( r ) {
 		ERROR("Scaling failed.\n");
+		*G = 1.0;
 		return 1;
 	}
 
@@ -567,7 +571,7 @@ int linear_scale(const RefList *list1, const RefList *list2, double *G,
 				int i;
 				for ( i=0; i<n; i++ ) {
 					STATUS("%i %e %e %e\n", i, x[i], y[i], w[n]);
-					}
+				}
 			}
 		}
 
@@ -582,3 +586,16 @@ int linear_scale(const RefList *list1, const RefList *list2, double *G,
 	return 0;
 }
 
+
+int scale_one_crystal(Crystal *cr, const RefList *list2, int flags)
+{
+	double G, B;
+	int r;
+
+	r = scale_one(crystal_get_reflections(cr), list2, flags, &G, &B);
+	if ( r ) return r;
+
+	crystal_set_osf(cr, G);
+	crystal_set_Bfac(cr, B);
+	return 0;
+}
diff --git a/src/scaling.h b/src/scaling.h
index 2161f69c..25aca31e 100644
--- a/src/scaling.h
+++ b/src/scaling.h
@@ -3,11 +3,11 @@
  *
  * Scaling
  *
- * Copyright © 2012-2017 Deutsches Elektronen-Synchrotron DESY,
+ * Copyright © 2012-2018 Deutsches Elektronen-Synchrotron DESY,
  *                       a research centre of the Helmholtz Association.
  *
  * Authors:
- *   2010-2017 Thomas White <taw@physics.org>
+ *   2010-2018 Thomas White <taw@physics.org>
  *
  * This file is part of CrystFEL.
  *
@@ -38,12 +38,20 @@
 #include "crystal.h"
 #include "geometry.h"
 
+enum ScaleFlags
+{
+	SCALE_NO_B,  /* Don't use Debye-Waller part */
+};
+
 extern double log_residual(Crystal *cr, const RefList *full, int free,
                            int *pn_used, const char *filename);
 
-extern int linear_scale(const RefList *list1, const RefList *list2, double *G,
-                        int complain_loudly);
+extern int scale_one(const RefList *list1, const RefList *list2, int flags,
+                     double *G, double *B);
+
+extern int scale_one_crystal(Crystal *cr, const RefList *list2, int flags);
 
-extern void scale_all(Crystal **crystals, int n_crystals, int nthreads);
+extern void scale_all(Crystal **crystals, int n_crystals, int nthreads,
+                      int flags);
 
 #endif	/* SCALING_H */
diff --git a/tests/linear_scale_check.c b/tests/linear_scale_check.c
index 4cb6a1ba..5c723849 100644
--- a/tests/linear_scale_check.c
+++ b/tests/linear_scale_check.c
@@ -3,11 +3,11 @@
  *
  * Check that linear scaling works
  *
- * Copyright © 2017 Deutsches Elektronen-Synchrotron DESY,
- *                  a research centre of the Helmholtz Association.
+ * Copyright © 2017-2018 Deutsches Elektronen-Synchrotron DESY,
+ *                       a research centre of the Helmholtz Association.
  *
  * Authors:
- *   2017 Thomas White <taw@physics.org>
+ *   2017-2018 Thomas White <taw@physics.org>
  *
  * This file is part of CrystFEL.
  *
@@ -45,7 +45,7 @@ int main(int argc, char *argv[])
 	gsl_rng *rng;
 	RefList *list1;
 	RefList *list2;
-	double G;
+	double G, B;
 	int r;
 
 	list1 = reflist_new();
@@ -70,7 +70,7 @@ int main(int argc, char *argv[])
 		set_partiality(refl2, 1.0);
 	}
 
-	r = linear_scale(list1, list2, &G, 1);
+	r = scale_one(list1, list2, SCALE_NO_B, &G, &B);
 	STATUS("Scaling result: %i, G = %f\n", r, G);
 
 	return fail;