Merge branch 'tom/partials'

5 files changed, 102 insertions, 51 deletions

diff --git a/src/merge.c b/src/merge.c
index 56ca2020..8017182b 100644
--- a/src/merge.c
+++ b/src/merge.c
@@ -187,7 +187,7 @@ static void run_merge_job(void *vwargs, int cookie)
 		Reflection *f;
 		signed int h, k, l;
 		double mean, sumweight, M2, temp, delta, R;
-		double corr, res, w;
+		double res, w;
 		struct reflection_contributions *c;
 
 		if ( get_partiality(refl) < MIN_PART_MERGE ) continue;
@@ -217,25 +217,16 @@ static void run_merge_job(void *vwargs, int cookie)
 			continue;
 		}
 
-		/* Total (multiplicative) correction factor */
-		corr = 1.0/G * exp(B*res*res) * get_lorentz(refl)
-		        / get_partiality(refl);
-		if ( isnan(corr) ) {
-			ERROR("NaN corr:\n");
-			ERROR("G = %f, B = %e\n", G, B);
-			ERROR("res = %e\n", res);
-			ERROR("p = %f\n", get_partiality(refl));
-			unlock_reflection(f);
-			continue;
-		}
-
 		/* Reflections count less the more they have to be scaled up */
-		w = 1.0;
+		w = get_partiality(refl);
 
 		/* Running mean and variance calculation */
 		temp = w + sumweight;
-		if (ln_merge) delta = log(get_intensity(refl)*corr) - mean;
-		else delta = get_intensity(refl)*corr - mean;
+		if ( ln_merge ) {
+			delta = log(correct_reflection(refl, G, B, res)) - mean;
+		} else {
+			delta = correct_reflection(refl, G,  B, res) - mean;
+		}
 		R = delta * w / temp;
 		set_intensity(f, mean + R);
 		set_temp2(f, M2 + sumweight * delta * R);
@@ -352,6 +343,25 @@ RefList *merge_intensities(Crystal **crystals, int n, int n_threads,
 }
 
 
+/* Correct intensity in pattern for scaling and Lorentz factors,
+ *  but not partiality nor polarisation */
+double correct_reflection_nopart(Reflection *refl, double osf, double Bfac,
+                                 double res)
+{
+	double corr = osf * exp(-Bfac*res*res);
+	return (get_intensity(refl) / corr) / get_lorentz(refl);
+}
+
+
+/* Correct intensity in pattern for scaling, partiality and Lorentz factors,
+ *  but not polarisation */
+double correct_reflection(Reflection *refl, double osf, double Bfac, double res)
+{
+	double Ipart = correct_reflection_nopart(refl, osf, Bfac, res);
+	return Ipart / get_partiality(refl);
+}
+
+
 double residual(Crystal *cr, const RefList *full, int free,
                 int *pn_used, const char *filename)
 {
@@ -368,7 +378,7 @@ double residual(Crystal *cr, const RefList *full, int free,
 	      refl != NULL;
 	      refl = next_refl(refl, iter) )
 	{
-		double p, w, corr, res;
+		double w, res;
 		signed int h, k, l;
 		Reflection *match;
 		double I_full;
@@ -384,13 +394,9 @@ double residual(Crystal *cr, const RefList *full, int free,
 
 		if ( get_redundancy(match) < 2 ) continue;
 
-		p = get_partiality(refl);
-		//if ( p < 0.2 ) continue;
-
-		corr = get_lorentz(refl) / (G * exp(-B*res*res));
-		int1 = get_intensity(refl) * corr;
-		int2 = p*I_full;
-		w = p;
+		int1 = correct_reflection_nopart(refl, G, B, res);
+		int2 = get_partiality(refl)*I_full;
+		w = 1.0;
 
 		num += fabs(int1 - int2) * w;
 		den += fabs(int1 + int2) * w/2.0;
@@ -454,6 +460,7 @@ double log_residual(Crystal *cr, const RefList *full, int free,
 		fx = log(G) - B*s*s + log(p) + log(I_full) - log(I_partial) - log(L);
 		w = 1.0;
 		dev += w*fx*fx;
+		n_used++;
 
 		if ( fh != NULL ) {
 			fprintf(fh, "%4i %4i %4i %e %e\n",
diff --git a/src/merge.h b/src/merge.h
index 476f8755..2c2cfeda 100644
--- a/src/merge.h
+++ b/src/merge.h
@@ -47,6 +47,12 @@ extern RefList *merge_intensities(Crystal **crystals, int n, int n_threads,
                                   int min_meas, double push_res, int use_weak,
                                   int ln_merge);
 
+extern double correct_reflection_nopart(Reflection *refl, double osf,
+                                        double Bfac, double res);
+
+extern double correct_reflection(Reflection *refl, double osf, double Bfac,
+                                 double res);
+
 extern double residual(Crystal *cr, const RefList *full, int free,
                        int *pn_used, const char *filename);
 
diff --git a/src/partialator.c b/src/partialator.c
index 86431e4c..8cf1dad9 100644
--- a/src/partialator.c
+++ b/src/partialator.c
@@ -638,7 +638,7 @@ static void write_to_pgraph(FILE *fh, RefList *list, RefList *full, Crystal *cr,
 	{
 		signed int h, k, l;
 		double pobs, pcalc;
-		double res, corr, Ipart;
+		double res, Ipart;
 		Reflection *match;
 
 		if ( !get_flag(refl) ) continue;  /* Not free-flagged */
@@ -660,8 +660,7 @@ static void write_to_pgraph(FILE *fh, RefList *list, RefList *full, Crystal *cr,
 		pcalc = get_partiality(refl);
 
 		/* Observed partiality */
-		corr = G * exp(B*res*res) * get_lorentz(refl);
-		Ipart = get_intensity(refl) * corr;
+		Ipart = correct_reflection_nopart(refl, G, B, res);
 		pobs = Ipart / get_intensity(match);
 
 		fprintf(fh, "%5i %4i %4i %4i %e %e %8.3f %8.3f %s\n",
@@ -716,6 +715,10 @@ static void all_residuals(Crystal **crystals, int n_crystals, RefList *full,
 	int n_nan_linear_free = 0;
 	int n_nan_log = 0;
 	int n_nan_log_free = 0;
+	int n_non_linear = 0;
+	int n_non_linear_free = 0;
+	int n_non_log = 0;
+	int n_non_log_free = 0;
 
 	*presidual = 0.0;
 	*pfree_residual = 0.0;
@@ -725,19 +728,35 @@ static void all_residuals(Crystal **crystals, int n_crystals, RefList *full,
 	for ( i=0; i<n_crystals; i++ ) {
 
 		double r, free_r, log_r, free_log_r;
+		int n;
 
 		if ( crystal_get_user_flag(crystals[i]) ) continue;
 
 		/* Scaling should have been done right before calling this */
-		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) ) n_nan_linear++;
-		if ( isnan(free_r) ) n_nan_linear_free++;
-		if ( isnan(log_r) ) n_nan_log++;
-		if ( isnan(free_log_r) ) n_nan_log_free++;
+		r = residual(crystals[i], full, 0, &n, NULL);
+		if ( n == 0 ) {
+			n_non_linear++;
+		} else if ( isnan(r) ) {
+			n_nan_linear++;
+		}
+		free_r = residual(crystals[i], full, 1, &n, NULL);
+		if ( n == 0 ) {
+			n_non_linear_free++;
+		} else if ( isnan(free_r) ) {
+			n_nan_linear_free++;
+		}
+		log_r = log_residual(crystals[i], full, 0, &n, NULL);
+		if ( n == 0 ) {
+			n_non_log++;
+		} else if ( isnan(log_r) ) {
+			n_nan_log++;
+		}
+		free_log_r = log_residual(crystals[i], full, 1, &n, NULL);
+		if ( n == 0 ) {
+			n_non_log_free++;
+		} else if ( isnan(free_log_r) ) {
+			n_nan_log_free++;
+		}
 
 		if ( isnan(r) || isnan(free_r)
 		  || isnan(log_r) || isnan(free_log_r) ) continue;
@@ -750,6 +769,23 @@ static void all_residuals(Crystal **crystals, int n_crystals, RefList *full,
 		n_used++;
 	}
 
+	if ( n_non_linear ) {
+		ERROR("WARNING: %i crystals had no reflections in linear "
+		      "residual calculation\n", n_non_linear);
+	}
+	if ( n_non_linear_free ) {
+		ERROR("WARNING: %i crystals had no reflections in linear free "
+		      "residual calculation\n", n_non_linear_free);
+	}
+	if ( n_non_log ) {
+		ERROR("WARNING: %i crystals had no reflections in log "
+		      "residual calculation\n", n_non_log);
+	}
+	if ( n_non_log_free ) {
+		ERROR("WARNING: %i crystals had no reflections in log free "
+		      "residual calculation\n", n_non_log_free);
+	}
+
 	if ( n_nan_linear ) {
 		ERROR("WARNING: %i crystals had NaN linear residuals\n",
 		      n_nan_linear);
diff --git a/src/post-refinement.c b/src/post-refinement.c
index 0a799b9e..b1af92db 100644
--- a/src/post-refinement.c
+++ b/src/post-refinement.c
@@ -559,7 +559,7 @@ void write_specgraph(Crystal *crystal, const RefList *full,
 	      refl != NULL;
 	      refl = next_refl(refl, iter) )
 	{
-		double corr, Ipart, Ifull, pobs, pcalc;
+		double Ipart, Ifull, pobs, pcalc;
 		double res;
 		signed int h, k, l;
 		Reflection *match;
@@ -573,8 +573,7 @@ void write_specgraph(Crystal *crystal, const RefList *full,
 		/* Don't calculate pobs if reference reflection is weak */
 		if ( fabs(get_intensity(match)) / get_esd_intensity(match) < 3.0 ) continue;
 
-		corr = G * exp(B*res*res) * get_lorentz(refl);
-		Ipart = get_intensity(refl) * corr;
+		Ipart = correct_reflection_nopart(refl, G, B, res);
 		Ifull = get_intensity(match);
 		pobs = Ipart / Ifull;
 		pcalc = get_partiality(refl);
diff --git a/src/rejection.c b/src/rejection.c
index b3a3cb67..5d1041fa 100644
--- a/src/rejection.c
+++ b/src/rejection.c
@@ -143,11 +143,7 @@ static int calculate_refl_mean_var(RefList *full)
 
 			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]);
+			Ii = correct_reflection(c->contribs[j], G, B, res);
 
 			Ex += Ii - K;
 			Ex2 += (Ii - K) * (Ii - K);
@@ -232,16 +228,13 @@ static double calculate_cchalf(RefList *template, RefList *full,
 		while ( exrefl != NULL ) {
 
 			double G, B;
-			double Ii = get_intensity(exrefl);
 
 			G = crystal_get_osf(exclude);
 			B = crystal_get_Bfac(exclude);
 
 			if ( get_partiality(exrefl) > MIN_PART_MERGE ) {
 
-				/* Total (multiplicative) correction factor */
-				Ii *= 1.0/G * exp(B*res*res) * get_lorentz(exrefl)
-				          / get_partiality(exrefl);
+				double Ii = correct_reflection(exrefl, G, B, res);
 
 				/* Remove contribution of this reflection */
 				Ex -= Ii - K;
@@ -289,6 +282,7 @@ static void check_deltacchalf(Crystal **crystals, int n, RefList *full)
 	double mean, sd;
 	int nref = 0;
 	int nnan = 0;
+	int nnon = 0;
 
 	if ( calculate_refl_mean_var(full) ) {
 		STATUS("No reflection contributions for deltaCChalf "
@@ -315,13 +309,22 @@ static void check_deltacchalf(Crystal **crystals, int n, RefList *full)
 		//STATUS("Without = %f", cchalfi*100.0);
 		//STATUS("  Delta = %f  ", (cchalf - cchalfi)*100.0);
 		//STATUS("(nref = %i)\n", nref);
-		vals[i] = cchalf - cchalfi;
-		if ( isnan(vals[i]) || isinf(vals[i]) ) {
+		if ( nref == 0 ) {
 			vals[i] = 0.0;
-			nnan++;
+			nnon++;
+		} else {
+			vals[i] = cchalf - cchalfi;
+			if ( isnan(vals[i]) || isinf(vals[i]) ) {
+				vals[i] = 0.0;
+				nnan++;
+			}
 		}
 		progress_bar(i, n-1, "Calculating deltaCChalf");
 	}
+	if ( nnon > 0 ) {
+		STATUS("WARNING: %i patterns had no reflections in deltaCChalf "
+		       "calculation (I set deltaCChalf=zero for them)\n", nnon);
+	}
 	if ( nnan > 0 ) {
 		STATUS("WARNING: %i NaN or inf deltaCChalf values were "
 		       "replaced with zero\n", nnan);