partialator: Rationalise residual calculation and display

3 files changed, 64 insertions, 118 deletions

diff --git a/src/partialator.c b/src/partialator.c
index de08df34..81c5ab37 100644
--- a/src/partialator.c
+++ b/src/partialator.c
@@ -585,6 +585,20 @@ static void normalise_scales(Crystal **crystals, int n_crystals)
 }
 
 
+static void show_all_residuals(Crystal **crystals, int n_crystals,
+                               RefList *full)
+{
+	double dev, free_dev, log_dev, free_log_dev;
+
+	all_residuals(crystals, n_crystals, full,
+	              &dev, &free_dev, &log_dev, &free_log_dev);
+	STATUS("Residuals:"
+	       "  linear          linear free     log             log free\n");
+	STATUS("         ");
+	STATUS("%15e %15e %15e %15e\n", dev, free_dev, log_dev, free_log_dev);
+}
+
+
 int main(int argc, char *argv[])
 {
 	int c;
@@ -990,29 +1004,15 @@ int main(int argc, char *argv[])
 	/* Iterate */
 	for ( i=0; i<n_iter; i++ ) {
 
-		double init_dev, init_free_dev;
-		double init_log_dev, init_free_log_dev;
-		double final_dev, final_free_dev;
-		double final_log_dev, final_free_log_dev;
+		show_all_residuals(crystals, n_crystals, full);
 
 		STATUS("Refinement cycle %i of %i\n", i+1, n_iter);
 
 		/* Refine all crystals to get the best fit */
 		refine_all(crystals, n_crystals, full, nthreads, pmodel,
-		           no_scale, no_pr, max_B,
-		           &init_dev, &init_free_dev,
-		           &init_log_dev, &init_free_log_dev,
-		           &final_dev, &final_free_dev,
-		           &final_log_dev, &final_free_log_dev);
-
-		STATUS("Overall residual: initial = %e, final = %e\n",
-		       init_dev, final_dev);
-		STATUS("Overall free residual: initial = %e, final = %e\n",
-		       init_free_dev, final_free_dev);
-		STATUS("Overall log residual: initial = %e, final = %e\n",
-		       init_log_dev, final_log_dev);
-		STATUS("Overall log free residual: initial = %e, final = %e\n",
-		       init_free_log_dev, final_free_log_dev);
+		           no_scale, no_pr, max_B);
+
+		show_all_residuals(crystals, n_crystals, full);
 
 		check_rejection(crystals, n_crystals, full);
 		normalise_scales(crystals, n_crystals);
@@ -1026,6 +1026,8 @@ int main(int argc, char *argv[])
 
 	}
 
+	show_all_residuals(crystals, n_crystals, full);
+
 	/* Output results */
 	STATUS("Writing overall results to %s\n", outfile);
 	write_reflist_2(outfile, full, sym);
diff --git a/src/post-refinement.c b/src/post-refinement.c
index 9d1aa7a2..818911df 100644
--- a/src/post-refinement.c
+++ b/src/post-refinement.c
@@ -54,6 +54,12 @@
 /* Maximum number of iterations of NLSq to do for each image per macrocycle. */
 #define MAX_CYCLES (10)
 
+struct prdata
+{
+	int refined;
+	int n_filtered;
+};
+
 const char *str_prflag(enum prflag flag)
 {
 	switch ( flag ) {
@@ -774,6 +780,38 @@ static double residual(Crystal *cr, const RefList *full, int verbose, int free,
 }
 
 
+void all_residuals(Crystal **crystals, int n_crystals, RefList *full,
+                   double *presidual, double *pfree_residual,
+                   double *plog_residual, double *pfree_log_residual)
+{
+	int i;
+	*presidual = 0.0;
+	*pfree_residual = 0.0;
+	*plog_residual = 0.0;
+	*pfree_log_residual = 0.0;
+
+	for ( i=0; i<n_crystals; i++ ) {
+
+		double r, free_r, log_r, free_log_r;
+
+		if ( crystal_get_user_flag(crystals[i]) ) continue;
+
+		r = residual(crystals[i], full, 0, 0, NULL);
+		free_r = residual(crystals[i], full, 0, 1, NULL);
+		log_r = log_residual(crystals[i], full, 0);
+		free_log_r = log_residual(crystals[i], full, 1);
+
+		if ( isnan(r) || isnan(free_r)
+		  || isnan(log_r) || isnan(free_log_r) ) continue;
+
+		*presidual += r;
+		*pfree_residual += free_r;
+		*plog_residual += log_r;
+		*pfree_log_residual += free_log_r;
+	}
+}
+
+
 static void write_residual_graph(Crystal *cr, const RefList *full)
 {
 	int i;
@@ -928,11 +966,6 @@ static struct prdata pr_refine(Crystal *cr, const RefList *full,
 	prdata.refined = 0;
 	prdata.n_filtered = 0;
 
-	prdata.initial_residual = residual(cr, full, 0, 0, NULL);
-	prdata.initial_free_residual = residual(cr, full, 0, 1, NULL);
-	prdata.initial_log_residual = log_residual(cr, full, 0);
-	prdata.initial_free_log_residual = log_residual(cr, full, 1);
-
 	if ( !no_scale ) {
 		do_scale_refine(cr, full, pmodel, verbose);
 	}
@@ -954,15 +987,6 @@ static struct prdata pr_refine(Crystal *cr, const RefList *full,
 	if ( crystal_get_user_flag(cr) == 0 ) {
 		prdata.refined = 1;
 	}
-	prdata.final_residual = residual(cr, full, 0, 0, NULL);
-	prdata.final_free_residual = residual(cr, full, 0, 1, NULL);
-	prdata.final_log_residual = log_residual(cr, full, 0);
-	prdata.final_free_log_residual = log_residual(cr, full, 1);
-
-	if ( isnan(prdata.final_residual) ) {
-		STATUS("nan residual!  Image serial %i\n",
-		       crystal_get_image(cr)->serial);
-	}
 
 	return prdata;
 }
@@ -987,14 +1011,6 @@ struct queue_args
 	Crystal **crystals;
 	int n_crystals;
 	struct refine_args task_defaults;
-	double initial_residual;
-	double initial_free_residual;
-	double initial_log_residual;
-	double initial_free_log_residual;
-	double final_residual;
-	double final_free_residual;
-	double final_log_residual;
-	double final_free_log_residual;
 };
 
 
@@ -1027,30 +1043,9 @@ static void *get_image(void *vqargs)
 static void done_image(void *vqargs, void *task)
 {
 	struct queue_args *qa = vqargs;
-	struct refine_args *pargs = task;
-	struct prdata *prd = &pargs->prdata;
 
 	qa->n_done++;
 
-	if ( !isnan(prd->initial_residual)
-	  && !isnan(prd->initial_free_residual)
-	  && !isnan(prd->initial_log_residual)
-	  && !isnan(prd->initial_free_log_residual)
-	  && !isnan(prd->final_residual)
-	  && !isnan(prd->final_free_residual)
-	  && !isnan(prd->final_log_residual)
-	  && !isnan(prd->final_free_log_residual) )
-	{
-		qa->initial_residual += prd->initial_residual;
-		qa->initial_free_residual += prd->initial_free_residual;
-		qa->initial_log_residual += prd->initial_log_residual;
-		qa->initial_free_log_residual += prd->initial_free_log_residual;
-		qa->final_residual += prd->final_residual;
-		qa->final_free_residual += prd->final_free_residual;
-		qa->final_log_residual += prd->final_log_residual;
-		qa->final_free_log_residual += prd->final_free_log_residual;
-	}
-
 	progress_bar(qa->n_done, qa->n_crystals, "Refining");
 	free(task);
 }
@@ -1058,11 +1053,7 @@ static void done_image(void *vqargs, void *task)
 
 void refine_all(Crystal **crystals, int n_crystals,
                 RefList *full, int nthreads, PartialityModel pmodel,
-                int no_scale, int no_pr, double max_B,
-                double *initial_residual, double *initial_free_residual,
-                double *initial_log_residual, double *initial_free_log_residual,
-                double *final_residual, double *final_free_residual,
-                double *final_log_residual, double *final_free_log_residual)
+                int no_scale, int no_pr, double max_B)
 {
 	struct refine_args task_defaults;
 	struct queue_args qargs;
@@ -1075,41 +1066,16 @@ void refine_all(Crystal **crystals, int n_crystals,
 	task_defaults.no_scale = no_scale;
 	task_defaults.no_pr = no_pr;
 	task_defaults.max_B = max_B;
-	task_defaults.prdata.initial_residual = 0.0;
-	task_defaults.prdata.initial_free_residual = 0.0;
-	task_defaults.prdata.initial_log_residual = 0.0;
-	task_defaults.prdata.initial_free_log_residual = 0.0;
-	task_defaults.prdata.final_residual = 0.0;
-	task_defaults.prdata.final_free_residual = 0.0;
-	task_defaults.prdata.final_log_residual = 0.0;
-	task_defaults.prdata.final_free_log_residual = 0.0;
 
 	qargs.task_defaults = task_defaults;
 	qargs.n_started = 0;
 	qargs.n_done = 0;
 	qargs.n_crystals = n_crystals;
 	qargs.crystals = crystals;
-	qargs.initial_residual = 0.0;
-	qargs.initial_free_residual = 0.0;
-	qargs.initial_log_residual = 0.0;
-	qargs.initial_free_log_residual = 0.0;
-	qargs.final_residual = 0.0;
-	qargs.final_free_residual = 0.0;
-	qargs.final_log_residual = 0.0;
-	qargs.final_free_log_residual = 0.0;
 
 	/* Don't have threads which are doing nothing */
 	if ( n_crystals < nthreads ) nthreads = n_crystals;
 
 	run_threads(nthreads, refine_image, get_image, done_image,
 	            &qargs, n_crystals, 0, 0, 0);
-
-	*initial_residual = qargs.initial_residual;
-	*initial_free_residual = qargs.initial_free_residual;
-	*initial_log_residual = qargs.initial_log_residual;
-	*initial_free_log_residual = qargs.initial_free_log_residual;
-	*final_residual = qargs.final_residual;
-	*final_free_residual = qargs.final_free_residual;
-	*final_log_residual = qargs.final_log_residual;
-	*final_free_log_residual = qargs.final_free_log_residual;
 }
diff --git a/src/post-refinement.h b/src/post-refinement.h
index 9c0f9984..b3509501 100644
--- a/src/post-refinement.h
+++ b/src/post-refinement.h
@@ -43,24 +43,6 @@
 #include "geometry.h"
 
 
-struct prdata
-{
-	int refined;
-	int n_filtered;
-
-	/* Before refinement */
-	double initial_residual;
-	double initial_free_residual;
-	double initial_log_residual;
-	double initial_free_log_residual;
-
-	/* After refinement */
-	double final_residual;
-	double final_free_residual;
-	double final_log_residual;
-	double final_free_log_residual;
-};
-
 enum prflag
 {
 	PRFLAG_OK = 0,
@@ -76,18 +58,14 @@ extern const char *str_prflag(enum prflag flag);
 
 extern void refine_all(Crystal **crystals, int n_crystals,
                        RefList *full, int nthreads, PartialityModel pmodel,
-                       int no_scale, int no_pr, double max_B,
-                       double *initial_residual,
-                       double *initial_free_residual,
-                       double *initial_log_residual,
-                       double *initial_free_log_residual,
-                       double *final_residual,
-                       double *final_free_residual,
-                       double *final_log_residual,
-                       double *final_free_log_residual);
+                       int no_scale, int no_pr, double max_B);
 
 /* Exported so it can be poked by tests/pr_p_gradient_check */
 extern double gradient(Crystal *cr, int k, Reflection *refl,
                        PartialityModel pmodel);
 
+extern void all_residuals(Crystal **crystals, int n_crystals, RefList *full,
+                          double *residual, double *free_residual,
+                          double *log_residual, double *free_log_residual);
+
 #endif	/* POST_REFINEMENT_H */