From 516fe4e0275323b3439ceda0e8a92204b8479c6e Mon Sep 17 00:00:00 2001 From: Thomas White Date: Fri, 16 Nov 2018 16:07:00 +0100 Subject: Restore reflection weighting --- src/merge.c | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) (limited to 'src/merge.c') diff --git a/src/merge.c b/src/merge.c index 56ca2020..a655b2cc 100644 --- a/src/merge.c +++ b/src/merge.c @@ -230,7 +230,7 @@ static void run_merge_job(void *vwargs, int cookie) } /* 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; -- cgit v1.2.3 From f30cc81c78493b3e0299f15a2ebf38aad2ba8e28 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Thu, 22 Nov 2018 16:18:43 +0100 Subject: log_residual: Actually set n_used --- src/merge.c | 1 + 1 file changed, 1 insertion(+) (limited to 'src/merge.c') diff --git a/src/merge.c b/src/merge.c index a655b2cc..1aaa272f 100644 --- a/src/merge.c +++ b/src/merge.c @@ -454,6 +454,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", -- cgit v1.2.3 From 6e9655241ff531b300d38d3e2f5a1e5c7d69f7b3 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Thu, 22 Nov 2018 16:49:49 +0100 Subject: Don't weight reflections by partiality in residual Because we want under-prediction to be penalised just as much as over-prediction --- src/merge.c | 9 +++------ 1 file changed, 3 insertions(+), 6 deletions(-) (limited to 'src/merge.c') diff --git a/src/merge.c b/src/merge.c index 1aaa272f..787afca3 100644 --- a/src/merge.c +++ b/src/merge.c @@ -368,7 +368,7 @@ double residual(Crystal *cr, const RefList *full, int free, refl != NULL; refl = next_refl(refl, iter) ) { - double p, w, corr, res; + double w, corr, res; signed int h, k, l; Reflection *match; double I_full; @@ -384,13 +384,10 @@ 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; + int2 = get_partiality(refl)*I_full; + w = 1.0; num += fabs(int1 - int2) * w; den += fabs(int1 + int2) * w/2.0; -- cgit v1.2.3 From a977dfd051fb0483cf27693fc8e950c99e7c88ed Mon Sep 17 00:00:00 2001 From: Thomas White Date: Fri, 30 Nov 2018 17:04:13 +0100 Subject: Factorise correction of intensity for G, B, p and L There were no fewer than 8 different places in the code where these factors needed to be applied. They all need to agree on conventions such as whether the intensities in the pattern should be multiplied or divided by G to "correct" them. They were not. This commit reduces the number of places to three: one function (actually two functions, so that the value before partiality can also be calculated consistently), plus log_residual and scale_one_crystal, where slightly different logarithmic versions are used instead. --- src/merge.c | 45 +++++++++++++++++++++++++++------------------ 1 file changed, 27 insertions(+), 18 deletions(-) (limited to 'src/merge.c') diff --git a/src/merge.c b/src/merge.c index 787afca3..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 = 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 w, corr, res; + double w, res; signed int h, k, l; Reflection *match; double I_full; @@ -384,8 +394,7 @@ double residual(Crystal *cr, const RefList *full, int free, if ( get_redundancy(match) < 2 ) continue; - corr = get_lorentz(refl) / (G * exp(-B*res*res)); - int1 = get_intensity(refl) * corr; + int1 = correct_reflection_nopart(refl, G, B, res); int2 = get_partiality(refl)*I_full; w = 1.0; -- cgit v1.2.3