Merge branch 'tom/deltacchalf'

7 files changed, 232 insertions, 65 deletions

diff --git a/src/merge.c b/src/merge.c
index 14460953..9d36a359 100644
--- a/src/merge.c
+++ b/src/merge.c
@@ -3,11 +3,11 @@
  *
  * Parallel weighted merging of intensities
  *
- * Copyright © 2012-2015 Deutsches Elektronen-Synchrotron DESY,
+ * Copyright © 2012-2018 Deutsches Elektronen-Synchrotron DESY,
  *                       a research centre of the Helmholtz Association.
  *
  * Authors:
- *   2010-2015 Thomas White <taw@physics.org>
+ *   2010-2018 Thomas White <taw@physics.org>
  *
  * This file is part of CrystFEL.
  *
@@ -46,6 +46,7 @@
 #include "cell.h"
 #include "utils.h"
 #include "reflist.h"
+#include "reflist-utils.h"
 #include "cell-utils.h"
 
 
@@ -88,6 +89,17 @@ static void *create_merge_job(void *vqargs)
 }
 
 
+static int alloc_contribs(struct reflection_contributions *c)
+{
+	c->contribs = realloc(c->contribs, c->max_contrib*sizeof(Reflection *));
+	c->contrib_crystals = realloc(c->contrib_crystals,
+	                              c->max_contrib*sizeof(Crystal *));
+	if ( c->contribs == NULL ) return 1;
+	if ( c->contrib_crystals == NULL ) return 1;
+	return 0;
+}
+
+
 /* Find reflection hkl in 'list', creating it if it's not there, under
  * protection of 'lock' and returning a locked reflection */
 static Reflection *get_locked_reflection(RefList *list, pthread_rwlock_t *lock,
@@ -108,12 +120,31 @@ static Reflection *get_locked_reflection(RefList *list, pthread_rwlock_t *lock,
 		 * So, we must check again */
 		f = find_refl(list, h, k, l);
 		if ( f == NULL ) {
+
+			struct reflection_contributions *c;
+
 			f = add_refl(list, h, k, l);
 			lock_reflection(f);
 			pthread_rwlock_unlock(lock);
 			set_intensity(f, 0.0);
 			set_temp1(f, 0.0);
 			set_temp2(f, 0.0);
+
+			c = malloc(sizeof(struct reflection_contributions));
+			if ( c != NULL ) {
+				c->n_contrib = 0;
+				c->max_contrib = 32;
+				c->contribs = NULL;
+				c->contrib_crystals = NULL;
+				if ( alloc_contribs(c) ) {
+					set_contributions(f, NULL);
+				} else {
+					set_contributions(f, c);
+				}
+			} else {
+				set_contributions(f, NULL);
+			}
+
 		} else {
 			/* Someone else created it */
 			lock_reflection(f);
@@ -158,6 +189,7 @@ static void run_merge_job(void *vwargs, int cookie)
 		signed int h, k, l;
 		double mean, sumweight, M2, temp, delta, R;
 		double corr, res, w;
+		struct reflection_contributions *c;
 
 		if ( get_partiality(refl) < MIN_PART_MERGE ) continue;
 
@@ -209,6 +241,18 @@ static void run_merge_job(void *vwargs, int cookie)
 		set_temp2(f, M2 + sumweight * delta * R);
 		set_temp1(f, temp);
 		set_redundancy(f, get_redundancy(f)+1);
+
+		/* Record this contribution */
+		c = get_contributions(f);
+		if ( c != NULL ) {
+			c->contribs[c->n_contrib] = refl;
+			c->contrib_crystals[c->n_contrib++] = cr;
+			if ( c->n_contrib == c->max_contrib ) {
+				c->max_contrib += 64;
+				alloc_contribs(c);
+			}
+		} /* else, too bad! */
+
 		unlock_reflection(f);
 
 		wargs->n_reflections++;
@@ -274,8 +318,18 @@ RefList *merge_intensities(Crystal **crystals, int n, int n_threads,
 			Reflection *r2;
 
 			get_indices(refl, &h, &k, &l);
-			r2 =  add_refl(full2, h, k, l);
+			r2 = add_refl(full2, h, k, l);
 			copy_data(r2, refl);
+
+		} else {
+
+			/* We do not need the contribution list any more */
+			struct reflection_contributions *c;
+			c = get_contributions(refl);
+			free(c->contribs);
+			free(c->contrib_crystals);
+			free(c);
+
 		}
 	}
 
diff --git a/src/partialator.c b/src/partialator.c
index 113535e0..78b7173b 100644
--- a/src/partialator.c
+++ b/src/partialator.c
@@ -199,12 +199,14 @@ static void write_split(Crystal **crystals, int n_crystals, const char *outfile,
 
 	STATUS("Writing two-way split to %s ", tmp);
 	write_reflist_2(tmp, split, sym);
+	free_contribs(split);
 	reflist_free(split);
 	snprintf(tmp, 1024, "%s2", outfile);
 	split = merge_intensities(crystals2, n_crystals2, nthreads,
 		                  min_measurements, push_res, 1);
 	STATUS("and %s\n", tmp);
 	write_reflist_2(tmp, split, sym);
+	free_contribs(split);
 	reflist_free(split);
 }
 
@@ -291,6 +293,7 @@ static void write_custom_split(struct custom_split *csplit, int dsn,
 	split = merge_intensities(crystalsn, n_crystalsn, nthreads,
 		                  min_measurements, push_res, 1);
 	write_reflist_2(tmp, split, sym);
+	free_contribs(split);
 	reflist_free(split);
 
 	write_split(crystalsn, n_crystalsn, tmp, nthreads, pmodel,
@@ -318,6 +321,7 @@ static void show_help(const char *s)
 "      --no-scale             Disable scale factor (G, B) refinement.\n"
 "      --no-Bscale            Disable B factor scaling.\n"
 "      --no-pr                Disable orientation/physics refinement.\n"
+"      --no-deltacchalf       Disable rejection based on deltaCChalf.\n"
 "  -m, --model=<model>        Specify partiality model.\n"
 "      --min-measurements=<n> Minimum number of measurements to require.\n"
 "      --no-polarisation      Disable polarisation correction.\n"
@@ -912,6 +916,7 @@ int main(int argc, char *argv[])
 	int scaleflags = 0;
 	double min_res = 0.0;
 	int do_write_logs = 0;
+	int no_deltacchalf = 0;
 
 	/* Long options */
 	const struct option longopts[] = {
@@ -950,6 +955,7 @@ int main(int argc, char *argv[])
 		{"no-free",            0, &no_free,            1},
 		{"output-every-cycle", 0, &output_everycycle,  1},
 		{"no-logs",            0, &no_logs,            1},
+		{"no-deltacchalf",     0, &no_deltacchalf,     1},
 
 		{0, 0, NULL, 0}
 	};
@@ -1412,7 +1418,7 @@ int main(int argc, char *argv[])
 	}
 
 	/* Check rejection and write figures of merit */
-	check_rejection(crystals, n_crystals, full, max_B);
+	check_rejection(crystals, n_crystals, full, max_B, no_deltacchalf);
 	show_all_residuals(crystals, n_crystals, full);
 
 	if ( do_write_logs ) {
@@ -1433,6 +1439,7 @@ int main(int argc, char *argv[])
 
 		/* Create new reference if needed */
 		if ( reference == NULL ) {
+			free_contribs(full);
 			reflist_free(full);
 			if ( !no_scale ) {
 				scale_all(crystals, n_crystals, nthreads,
@@ -1443,7 +1450,7 @@ int main(int argc, char *argv[])
 			                         push_res, 1);
 		} /* else full still equals reference */
 
-		check_rejection(crystals, n_crystals, full, max_B);
+		check_rejection(crystals, n_crystals, full, max_B, no_deltacchalf);
 		show_all_residuals(crystals, n_crystals, full);
 
 		if ( do_write_logs ) {
@@ -1481,6 +1488,7 @@ int main(int argc, char *argv[])
 	/* Final merge */
 	STATUS("Final merge...\n");
 	if ( reference == NULL ) {
+		free_contribs(full);
 		reflist_free(full);
 		if ( !no_scale ) {
 			scale_all(crystals, n_crystals, nthreads, scaleflags);
@@ -1527,6 +1535,7 @@ int main(int argc, char *argv[])
 		reflist_free(crystal_get_reflections(crystals[i]));
 		crystal_free(crystals[i]);
 	}
+	free_contribs(full);
 	reflist_free(full);
 	free_symoplist(sym);
 	free(outfile);
diff --git a/src/post-refinement.c b/src/post-refinement.c
index bde9c02c..0a799b9e 100644
--- a/src/post-refinement.c
+++ b/src/post-refinement.c
@@ -69,8 +69,8 @@ const char *str_prflag(enum prflag flag)
 		case PRFLAG_EARLY :
 		return "early rejection";
 
-		case PRFLAG_CC :
-		return "low CC";
+		case PRFLAG_DELTACCHALF :
+		return "negative delta CC½";
 
 		case PRFLAG_BIGB :
 		return "B too big";
diff --git a/src/post-refinement.h b/src/post-refinement.h
index 71a6d7f3..e2873b76 100644
--- a/src/post-refinement.h
+++ b/src/post-refinement.h
@@ -50,7 +50,7 @@ enum prflag
 	PRFLAG_FEWREFL = 16,
 	PRFLAG_SOLVEFAIL = 17,
 	PRFLAG_EARLY = 18,
-	PRFLAG_CC = 19,
+	PRFLAG_DELTACCHALF = 19,
 	PRFLAG_BIGB = 20,
 	PRFLAG_SCALEBAD = 21,
 };
diff --git a/src/rejection.c b/src/rejection.c
index 9d41f9b4..230093e6 100644
--- a/src/rejection.c
+++ b/src/rejection.c
@@ -3,11 +3,11 @@
  *
  * Crystal rejection for scaling
  *
- * Copyright © 2012-2015 Deutsches Elektronen-Synchrotron DESY,
+ * Copyright © 2012-2018 Deutsches Elektronen-Synchrotron DESY,
  *                       a research centre of the Helmholtz Association.
  *
  * Authors:
- *   2010-2015 Thomas White <taw@physics.org>
+ *   2010-2018 Thomas White <taw@physics.org>
  *
  * This file is part of CrystFEL.
  *
@@ -94,74 +94,174 @@ void early_rejection(Crystal **crystals, int n)
 }
 
 
-static void check_cc(Crystal *cr, RefList *full)
+static double calculate_cchalf(RefList *template, RefList *full,
+                               Crystal *exclude, int *pnref)
 {
-	RefList *list = crystal_get_reflections(cr);
-	Reflection *refl;
+	Reflection *trefl;
 	RefListIterator *iter;
-	double G = crystal_get_osf(cr);
-	double B = crystal_get_Bfac(cr);
-	UnitCell *cell = crystal_get_cell(cr);
-	double *vec1, *vec2;
-	int n, i;
-	double cc;
-
-	n = num_reflections(list);
-	vec1 = malloc(n*sizeof(double));
-	vec2 = malloc(n*sizeof(double));
-	if ( (vec1 == NULL) || (vec2 == NULL) ) {
-		ERROR("Not enough memory to check CCs\n");
-		return;
-	}
-
-	i = 0;
-	for ( refl = first_refl(list, &iter);
-	      refl != NULL;
-	      refl = next_refl(refl, iter) )
+	double sig2E, sig2Y;
+	int n = 0;
+	signed int oh = 0;
+	signed int ok = 0;
+	signed int ol = 0;
+	double wSum = 0.0;
+	double wSum2 = 0.0;
+	double mean = 0.0;
+	double S = 0.0;
+	double all_sum_var = 0.0;
+	long int total_contribs = 0;
+
+	/* Iterate over all reflections */
+	for ( trefl = first_refl(template, &iter);
+	      trefl != NULL;
+	      trefl = next_refl(trefl, iter) )
 	{
+		struct reflection_contributions *c;
+		int j;
+		double refl_sum;
+		double refl_sumsq;
+		double refl_mean, refl_var;
 		signed int h, k, l;
-		double pobs, pcalc;
-		double res, corr, Ipart;
-		Reflection *match;
+		int nc = 0;
+		Reflection *refl;
+		double res;
+
+		get_indices(trefl, &h, &k, &l);
+
+		/* next_refl() will iterate over multiple copies of the same
+		 * unique reflection, but we are only interested in seeing each
+		 * unique reflection once. */
+		if ( (h==oh) && (k==ok) && (l==ol) ) continue;
+		oh = h;  ok = k;  ol = l;
+
+		refl = find_refl(full, h, k, l);
+		if ( refl == NULL ) continue;
+
+		c = get_contributions(refl);
+		assert(c != NULL);
 
-		if ( !get_flag(refl) ) continue;  /* Not free-flagged */
+		/* Calculate the resolution just once, using the cell from the
+		 * first crystal to contribute, otherwise it takes too long */
+		res = resolution(crystal_get_cell(c->contrib_crystals[0]),
+		                 h, k, l);
 
-		get_indices(refl, &h, &k, &l);
-		res = resolution(cell, h, k, l);
-		if ( 2.0*res > crystal_get_resolution_limit(cr) ) continue;
+		/* Mean of contributions */
+		refl_sum = 0.0;
+		for ( j=0; j<c->n_contrib; j++ ) {
 
-		match = find_refl(full, h, k, l);
-		if ( match == NULL ) continue;
+			double Ii, G, B;
 
-		/* Calculated partiality */
-		pcalc = get_partiality(refl);
+			if ( c->contrib_crystals[j] == exclude ) {
+				continue;
+			}
+
+			G = crystal_get_osf(c->contrib_crystals[j]);
+			B = crystal_get_Bfac(c->contrib_crystals[j]);
+			Ii = get_intensity(c->contribs[j]);
+
+			/* Total (multiplicative) correction factor */
+			Ii *= 1.0/G * exp(B*res*res) * get_lorentz(c->contribs[j])
+			          / get_partiality(c->contribs[j]);
+
+			refl_sum += Ii;
+			nc++;
+
+		}
+
+		if ( nc < 2 ) continue;
+		refl_mean = refl_sum / nc;
+
+		/* Variance of contributions */
+		refl_sumsq = 0.0;
+		for ( j=0; j<c->n_contrib; j++ ) {
+
+			double Ii, G, B;
+
+			if ( c->contrib_crystals[j] == exclude ) {
+				continue;
+			}
 
-		/* Observed partiality */
-		corr = G * exp(B*res*res) * get_lorentz(refl);
-		Ipart = get_intensity(refl) * corr;
-		pobs = Ipart / get_intensity(match);
+			G = crystal_get_osf(c->contrib_crystals[j]);
+			B = crystal_get_Bfac(c->contrib_crystals[j]);
+			Ii = get_intensity(c->contribs[j]);
+
+			/* Total (multiplicative) correction factor */
+			Ii *= 1.0/G * exp(B*res*res) * get_lorentz(c->contribs[j])
+			          / get_partiality(c->contribs[j]);
+
+			refl_sumsq += (Ii-refl_mean)*(Ii-refl_mean);
+			/* (nc already summed above) */
+
+		}
+
+		refl_var = refl_sumsq/(nc-1.0);
+		refl_var /= (nc/2.0);
+
+		all_sum_var += refl_var;
+		n++;
+		total_contribs += nc;
+
+		double w = 1.0;
+		wSum += w;
+		wSum2 += w*w;
+		double meanOld = mean;
+		mean = meanOld + (w/wSum) * (refl_mean - meanOld);
+		S += w * (refl_mean - meanOld) * (refl_mean - mean);
 
-		vec1[i] = pobs;
-		vec2[i] = pcalc;
-		i++;
 	}
 
-	cc = gsl_stats_correlation(vec1, 1, vec2, 1, i);
-	//printf("%f\n", cc);
-	if ( cc < 0.5 ) crystal_set_user_flag(cr, PRFLAG_CC);
+	sig2E = all_sum_var / n;
+	sig2Y = S / (wSum - 1.0);
 
-	free(vec1);
-	free(vec2);
+	if ( pnref != NULL ) {
+		*pnref = n;
+	}
+	return (sig2Y - 0.5*sig2E) / (sig2Y + 0.5*sig2E);
 }
 
 
-static UNUSED void check_ccs(Crystal **crystals, int n_crystals, RefList *full)
+static void check_deltacchalf(Crystal **crystals, int n, RefList *full)
 {
+	double cchalf;
 	int i;
+	int nref;
+	double *vals;
+	double mean, sd;
+
+	cchalf = calculate_cchalf(full, full, NULL, &nref);
+	STATUS("Overall CChalf = %f %% (%i reflections)\n", cchalf*100.0, nref);
+
+	vals = malloc(n*sizeof(double));
+	if ( vals == NULL ) {
+		ERROR("Not enough memory for deltaCChalf check\n");
+		return;
+	}
 
-	for ( i=0; i<n_crystals; i++ ) {
-		check_cc(crystals[i], full);
+	for ( i=0; i<n; i++ ) {
+		double cchalf, cchalfi;
+		RefList *template = crystal_get_reflections(crystals[i]);
+		cchalf = calculate_cchalf(template, full, NULL, NULL);
+		cchalfi = calculate_cchalf(template, full, crystals[i], &nref);
+		//STATUS("Frame %i:", i);
+		//STATUS("   With = %f  ", cchalf*100.0);
+		//STATUS("Without = %f", cchalfi*100.0);
+		//STATUS("  Delta = %f  ", (cchalf - cchalfi)*100.0);
+		//STATUS("(nref = %i)\n", nref);
+		vals[i] = cchalf - cchalfi;
+		progress_bar(i, n-1, "Calculating deltaCChalf");
 	}
+
+	mean = gsl_stats_mean(vals, 1, n);
+	sd = gsl_stats_sd_m(vals, 1, n, mean);
+	STATUS("deltaCChalf = %f ± %f %%\n", mean*100.0, sd*100.0);
+
+	for ( i=0; i<n; i++ ) {
+		if ( vals[i] < mean-2.0*sd ) {
+			crystal_set_user_flag(crystals[i], PRFLAG_DELTACCHALF);
+		}
+	}
+
+	free(vals);
 }
 
 
@@ -192,14 +292,16 @@ static void show_duds(Crystal **crystals, int n_crystals)
 }
 
 
-void check_rejection(Crystal **crystals, int n, RefList *full, double max_B)
+void check_rejection(Crystal **crystals, int n, RefList *full, double max_B,
+                     int no_deltacchalf)
 {
 	int i;
 	int n_acc = 0;
 
-	/* Check according to CCs FIXME: Disabled */
-	//if ( full != NULL ) check_ccs(crystals, n, full);
-
+	/* Check according to delta CC½ */
+	if ( !no_deltacchalf && (full != NULL) ) {
+		 check_deltacchalf(crystals, n, full);
+	}
 
 	for ( i=0; i<n; i++ ) {
 		if ( fabs(crystal_get_Bfac(crystals[i])) > max_B ) {
diff --git a/src/rejection.h b/src/rejection.h
index 8979313b..5a0c3e9f 100644
--- a/src/rejection.h
+++ b/src/rejection.h
@@ -3,11 +3,11 @@
  *
  * Crystal rejection for scaling
  *
- * Copyright © 2012-2015 Deutsches Elektronen-Synchrotron DESY,
+ * Copyright © 2012-2018 Deutsches Elektronen-Synchrotron DESY,
  *                       a research centre of the Helmholtz Association.
  *
  * Authors:
- *   2010-2015 Thomas White <taw@physics.org>
+ *   2010-2018 Thomas White <taw@physics.org>
  *
  * This file is part of CrystFEL.
  *
@@ -39,6 +39,6 @@
 
 extern void early_rejection(Crystal **crystals, int n);
 extern void check_rejection(Crystal **crystals, int n, RefList *full,
-                            double max_B);
+                            double max_B, int no_deltacchalf);
 
 #endif	/* REJECTION_H */
diff --git a/src/scaling.c b/src/scaling.c
index 21bd0db2..68456900 100644
--- a/src/scaling.c
+++ b/src/scaling.c
@@ -46,6 +46,7 @@
 #include "cell.h"
 #include "cell-utils.h"
 #include "scaling.h"
+#include "reflist-utils.h"
 
 
 struct scale_args
@@ -166,6 +167,7 @@ void scale_all(Crystal **crystals, int n_crystals, int nthreads, int scaleflags)
 		meanB /= n_crystals;
 		STATUS("Mean B = %e\n", meanB);
 
+		free_contribs(full);
 		reflist_free(full);
 		niter++;