partialator: Scale (strictly) using strong reflections only

1 files changed, 51 insertions, 28 deletions

diff --git a/src/partialator.c b/src/partialator.c
index 26c72696..6de69cd2 100644
--- a/src/partialator.c
+++ b/src/partialator.c
@@ -54,6 +54,7 @@
 #include <cell-utils.h>
 
 #include "version.h"
+#include "scaling.h"
 #include "post-refinement.h"
 #include "merge.h"
 #include "rejection.h"
@@ -564,23 +565,35 @@ static void write_pgraph(RefList *full, Crystal **crystals, int n_crystals,
 }
 
 
-static void normalise_scales(Crystal **crystals, int n_crystals)
+static void all_residuals(Crystal **crystals, int n_crystals, RefList *full,
+                          double *presidual, double *pfree_residual,
+                          double *plog_residual, double *pfree_log_residual)
 {
 	int i;
-	int n = 0;
-	double sum_B = 0.0;
-	double mean_B;
+
+	*presidual = 0.0;
+	*pfree_residual = 0.0;
+	*plog_residual = 0.0;
+	*pfree_log_residual = 0.0;
 
 	for ( i=0; i<n_crystals; i++ ) {
-		if ( crystal_get_user_flag(crystals[i]) != 0 ) continue;
-		sum_B += crystal_get_Bfac(crystals[i]);
-		n++;
-	}
-	mean_B = sum_B / n;
-	STATUS("Mean B = %e A^2\n", mean_B*1e20);
-	for ( i=0; i<n_crystals; i++ ) {
-		double B = crystal_get_Bfac(crystals[i]);
-		crystal_set_Bfac(crystals[i], B-mean_B);
+
+		double r, free_r, log_r, free_log_r;
+
+		if ( crystal_get_user_flag(crystals[i]) ) continue;
+
+		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);
+		free_log_r = log_residual(crystals[i], full, 1, NULL, NULL);
+
+		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;
 	}
 }
 
@@ -993,40 +1006,50 @@ int main(int argc, char *argv[])
 	STATUS("Checking patterns.\n");
 	//early_rejection(crystals, n_crystals);
 
-	/* Make initial estimates */
+	/* Initial rejection, figures of merit etc */
 	full = merge_intensities(crystals, n_crystals, nthreads, pmodel,
 	                         min_measurements, push_res, 1);
-
 	check_rejection(crystals, n_crystals, full);
-
+	show_all_residuals(crystals, n_crystals, full);
 	write_pgraph(full, crystals, n_crystals, 0);
 
 	/* Iterate */
 	for ( i=0; i<n_iter; i++ ) {
 
-		show_all_residuals(crystals, n_crystals, full);
-
-		STATUS("Refinement cycle %i of %i\n", i+1, n_iter);
+		STATUS("Scaling and 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);
+		if ( !no_scale ) {
+			scale_all(crystals, n_crystals, nthreads, pmodel,
+			          max_B);
+			reflist_free(full);
+			full = merge_intensities(crystals, n_crystals, nthreads,
+			                         pmodel, min_measurements,
+			                         push_res, 1);
+		}
 
-		show_all_residuals(crystals, n_crystals, full);
+		if ( !no_pr ) {
+			refine_all(crystals, n_crystals, full, nthreads,
+			           pmodel);
+			reflist_free(full);
+			full = merge_intensities(crystals, n_crystals, nthreads,
+			                         pmodel, min_measurements,
+			                         push_res, 1);
+		}
 
 		check_rejection(crystals, n_crystals, full);
-		normalise_scales(crystals, n_crystals);
-
-		/* Re-estimate all the full intensities */
 		reflist_free(full);
 		full = merge_intensities(crystals, n_crystals, nthreads,
-		                         pmodel, min_measurements, push_res, 1);
-
+		                         pmodel, min_measurements,
+		                         push_res, 1);
+		show_all_residuals(crystals, n_crystals, full);
 		write_pgraph(full, crystals, n_crystals, i+1);
 
 	}
 
+	full = merge_intensities(crystals, n_crystals, nthreads, pmodel,
+	                         min_measurements, push_res, 1);
 	show_all_residuals(crystals, n_crystals, full);
+	write_pgraph(full, crystals, n_crystals, n_iter+1);
 
 	/* Output results */
 	STATUS("Writing overall results to %s\n", outfile);