partialator: Improve rejection

3 files changed, 96 insertions, 11 deletions

diff --git a/src/partialator.c b/src/partialator.c
index dc0d647f..d5e1aef5 100644
--- a/src/partialator.c
+++ b/src/partialator.c
@@ -301,6 +301,7 @@ static void show_duds(Crystal **crystals, int n_crystals)
 	int n_noref = 0;
 	int n_solve = 0;
 	int n_early = 0;
+	int n_cc = 0;
 
 	for ( j=0; j<n_crystals; j++ ) {
 		int flag;
@@ -327,6 +328,10 @@ static void show_duds(Crystal **crystals, int n_crystals)
 			n_early++;
 			break;
 
+			case 6:
+			n_cc++;
+			break;
+
 			default:
 			STATUS("Unknown flag %i\n", flag);
 			break;
@@ -339,6 +344,7 @@ static void show_duds(Crystal **crystals, int n_crystals)
 		STATUS(" %i not enough reflections.\n", n_noref);
 		STATUS(" %i solve failed.\n", n_solve);
 		STATUS(" %i early rejection.\n", n_early);
+		STATUS(" %i bad CC.\n", n_cc);
 	}
 }
 
@@ -741,14 +747,14 @@ int main(int argc, char *argv[])
 	}
 
 	/* Make a first pass at cutting out crap */
-//	STATUS("Checking patterns.\n");
-//	early_rejection(crystals, n_crystals);
+	STATUS("Checking patterns.\n");
+	early_rejection(crystals, n_crystals);
 
 	/* Make initial estimates */
 	full = merge_intensities(crystals, n_crystals, nthreads, pmodel,
 	                         min_measurements, push_res);
 
-	check_rejection(crystals, n_crystals);
+	check_rejection(crystals, n_crystals, full);
 
 	show_duds(crystals, n_crystals);
 
@@ -773,14 +779,13 @@ int main(int argc, char *argv[])
 		       init_free_dev, final_free_dev);
 
 		show_duds(crystals, n_crystals);
-		check_rejection(crystals, n_crystals);
+		check_rejection(crystals, n_crystals, full);
 
 		/* Re-estimate all the full intensities */
 		reflist_free(full);
 		full = merge_intensities(crystals, n_crystals, nthreads,
 		                         pmodel, min_measurements, push_res);
 
-		check_rejection(crystals, n_crystals);
 		write_pgraph(full, crystals, n_crystals, i+1);
 
 	}
diff --git a/src/rejection.c b/src/rejection.c
index 0b59e813..28a8d32f 100644
--- a/src/rejection.c
+++ b/src/rejection.c
@@ -33,10 +33,12 @@
 
 #include <stdlib.h>
 #include <assert.h>
+#include <gsl/gsl_statistics.h>
 
 #include "crystal.h"
 #include "reflist.h"
 #include "rejection.h"
+#include "cell-utils.h"
 
 
 static double mean_intensity(RefList *list)
@@ -88,20 +90,98 @@ void early_rejection(Crystal **crystals, int n)
 
 	STATUS("Mean intensity/peak = %f ADU\n", m/n);
 	STATUS("%i crystals flagged\n", n_flag);
-	check_rejection(crystals, n);
+	check_rejection(crystals, n, NULL);
 }
 
 
-void check_rejection(Crystal **crystals, int n)
+static void check_cc(Crystal *cr, RefList *full)
+{
+	RefList *list = crystal_get_reflections(cr);
+	Reflection *refl;
+	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) )
+	{
+		signed int h, k, l;
+		double pobs, pcalc;
+		double res, corr, Ipart;
+		Reflection *match;
+
+		if ( !get_flag(refl) ) continue;  /* Not free-flagged */
+
+		get_indices(refl, &h, &k, &l);
+		res = resolution(cell, h, k, l);
+		if ( 2.0*res > crystal_get_resolution_limit(cr) ) continue;
+
+		match = find_refl(full, h, k, l);
+		if ( match == NULL ) continue;
+
+		/* Calculated partiality */
+		pcalc = get_partiality(refl);
+
+		/* Observed partiality */
+		corr = exp(-G) * exp(B*res*res) * get_lorentz(refl);
+		Ipart = get_intensity(refl) * corr;
+		pobs = Ipart / get_intensity(match);
+
+		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, 6);
+
+	free(vec1);
+	free(vec2);
+}
+
+
+static void check_ccs(Crystal **crystals, int n_crystals, RefList *full)
+{
+	int i;
+
+	for ( i=0; i<n_crystals; i++ ) {
+		check_cc(crystals[i], full);
+	}
+}
+
+
+void check_rejection(Crystal **crystals, int n, RefList *full)
 {
 	int i;
 	int n_acc = 0;
 
+	/* Check according to CCs FIXME: Disabled */
+	//if ( full != NULL ) check_ccs(crystals, n, full);
+
 	for ( i=0; i<n; i++ ) {
-		if ( crystal_get_user_flag(crystals[i]) == 0 ) {
-			n_acc++;
-			if ( n_acc >= 2 ) break;
+
+		/* Reject if B factor modulus is very large */
+		if ( fabs(crystal_get_Bfac(crystals[i])) > 1e18 ) {
+			crystal_set_user_flag(crystals[i], 1);
 		}
+
+		if ( crystal_get_user_flag(crystals[i]) == 0 ) n_acc++;
+
 	}
 
 	if ( n_acc < 2 ) {
diff --git a/src/rejection.h b/src/rejection.h
index fb73bd41..ec529941 100644
--- a/src/rejection.h
+++ b/src/rejection.h
@@ -38,6 +38,6 @@
 #include "crystal.h"
 
 extern void early_rejection(Crystal **crystals, int n);
-extern void check_rejection(Crystal **crystals, int n);
+extern void check_rejection(Crystal **crystals, int n, RefList *full);
 
 #endif	/* REJECTION_H */