From 4f5fe4915b36d30799f446e2e66ced25fc948d00 Mon Sep 17 00:00:00 2001
From: Thomas White <taw@physics.org>
Date: Wed, 6 May 2015 11:45:32 +0200
Subject: partialator: Always honour --min-measurements

---
 src/hrs-scaling.c | 25 +++++++++++++++++++++----
 src/hrs-scaling.h |  2 +-
 src/partialator.c |  5 +++--
 3 files changed, 25 insertions(+), 7 deletions(-)

(limited to 'src')

diff --git a/src/hrs-scaling.c b/src/hrs-scaling.c
index d5f93a3c..56a8ee16 100644
--- a/src/hrs-scaling.c
+++ b/src/hrs-scaling.c
@@ -354,9 +354,10 @@ static void finalise_merge_job(void *vqargs, void *vwargs)
 
 
 RefList *lsq_intensities(Crystal **crystals, int n, int n_threads,
-                         PartialityModel pmodel)
+                         PartialityModel pmodel, int min_meas)
 {
 	RefList *full;
+	RefList *full2;
 	struct merge_queue_args qargs;
 	Reflection *refl;
 	RefListIterator *iter;
@@ -374,6 +375,10 @@ RefList *lsq_intensities(Crystal **crystals, int n, int n_threads,
 
 	pthread_rwlock_destroy(&qargs.full_lock);
 
+	/* Calculate ESDs from variances, including only reflections with
+	 * enough measurements */
+	full2 = reflist_new();
+	if ( full2 == NULL ) return NULL;
 	for ( refl = first_refl(full, &iter);
 	      refl != NULL;
 	      refl = next_refl(refl, iter) )
@@ -384,9 +389,20 @@ RefList *lsq_intensities(Crystal **crystals, int n, int n_threads,
 		red = get_redundancy(refl);
 		var = get_temp2(refl) / get_temp1(refl);
 		set_esd_intensity(refl, sqrt(var)/sqrt(red));
+
+		if ( red >= min_meas ) {
+
+			signed int h, k, l;
+			Reflection *r2;
+
+			get_indices(refl, &h, &k, &l);
+			r2 =  add_refl(full2, h, k, l);
+			copy_data(r2, refl);
+		}
 	}
 
-	return full;
+	reflist_free(full);
+	return full2;
 }
 
 
@@ -488,7 +504,7 @@ RefList *scale_intensities(Crystal **crystals, int n, int n_threads,
 	}
 
 	/* Create an initial list to refine against */
-	full = lsq_intensities(crystals, n, n_threads, pmodel);
+	full = lsq_intensities(crystals, n, n_threads, pmodel, min_redundancy);
 
 	old_osfs = malloc(n*sizeof(double));
 	old_Bs = malloc(n*sizeof(double));
@@ -539,7 +555,8 @@ RefList *scale_intensities(Crystal **crystals, int n, int n_threads,
 
 		/* Generate list for next iteration */
 		reflist_free(full);
-		full = lsq_intensities(crystals, n, n_threads, pmodel);
+		full = lsq_intensities(crystals, n, n_threads, pmodel,
+		                       min_redundancy);
 
 		i++;
 
diff --git a/src/hrs-scaling.h b/src/hrs-scaling.h
index c41cb323..b7954527 100644
--- a/src/hrs-scaling.h
+++ b/src/hrs-scaling.h
@@ -44,6 +44,6 @@ extern RefList *scale_intensities(Crystal **crystals, int n, int n_threads,
 
 
 extern RefList *lsq_intensities(Crystal **crystals, int n, int n_threads,
-                                PartialityModel pmodel);
+                                PartialityModel pmodel, int min_meas);
 
 #endif	/* HRS_SCALING_H */
diff --git a/src/partialator.c b/src/partialator.c
index ca3a9976..9bdb684b 100644
--- a/src/partialator.c
+++ b/src/partialator.c
@@ -644,7 +644,8 @@ int main(int argc, char *argv[])
 	STATUS("Performing initial scaling.\n");
 	if ( noscale ) {
 		STATUS("Skipping scaling step (--no-scale).\n");
-		full = lsq_intensities(crystals, n_crystals, nthreads, pmodel);
+		full = lsq_intensities(crystals, n_crystals, nthreads, pmodel,
+		                       min_measurements);
 	} else {
 		full = scale_intensities(crystals, n_crystals, nthreads, pmodel,
 		                         min_measurements);
@@ -683,7 +684,7 @@ int main(int argc, char *argv[])
 		if ( noscale ) {
 			STATUS("Skipping scaling step (--no-scale).\n");
 			full = lsq_intensities(crystals, n_crystals, nthreads,
-			                       pmodel);
+			                       pmodel, min_measurements);
 		} else {
 			full = scale_intensities(crystals, n_crystals, nthreads,
 			                         pmodel, min_measurements);
-- 
cgit v1.2.3