From 7c8bfecf277d065c52c3c66ec061ceaf8df638e6 Mon Sep 17 00:00:00 2001 From: Ken Beyerlein Date: Fri, 11 Nov 2016 19:04:36 +0100 Subject: Merge ln of intensity rather than intensities themselves --- src/merge.c | 24 +++++++++++++++++++++--- src/merge.h | 3 ++- src/partialator.c | 14 +++++++------- src/scaling.c | 2 +- 4 files changed, 31 insertions(+), 12 deletions(-) diff --git a/src/merge.c b/src/merge.c index 14460953..da333f6a 100644 --- a/src/merge.c +++ b/src/merge.c @@ -62,6 +62,7 @@ struct merge_queue_args double push_res; int use_weak; long long int n_reflections; + int ln_merge; }; @@ -137,6 +138,7 @@ static void run_merge_job(void *vwargs, int cookie) Crystal *cr = wargs->crystal; RefList *full = wargs->qargs->full; double push_res = wargs->qargs->push_res; + int ln_merge = wargs->qargs->ln_merge; Reflection *refl; RefListIterator *iter; double G, B; @@ -161,7 +163,7 @@ static void run_merge_job(void *vwargs, int cookie) if ( get_partiality(refl) < MIN_PART_MERGE ) continue; - if ( !wargs->qargs->use_weak ) { + if ( !wargs->qargs->use_weak || ln_merge ) { if (get_intensity(refl) < 3.0*get_esd_intensity(refl)) { continue; @@ -203,7 +205,8 @@ static void run_merge_job(void *vwargs, int cookie) /* Running mean and variance calculation */ temp = w + sumweight; - delta = get_intensity(refl)*corr - mean; + if (ln_merge) delta = log(get_intensity(refl)*corr) - mean; + else delta = get_intensity(refl)*corr - mean; R = delta * w / temp; set_intensity(f, mean + R); set_temp2(f, M2 + sumweight * delta * R); @@ -228,7 +231,7 @@ static void finalise_merge_job(void *vqargs, void *vwargs) RefList *merge_intensities(Crystal **crystals, int n, int n_threads, int min_meas, - double push_res, int use_weak) + double push_res, int use_weak, int ln_merge) { RefList *full; RefList *full2; @@ -246,6 +249,7 @@ RefList *merge_intensities(Crystal **crystals, int n, int n_threads, qargs.push_res = push_res; qargs.use_weak = use_weak; qargs.n_reflections = 0; + qargs.ln_merge = ln_merge; pthread_rwlock_init(&qargs.full_lock, NULL); run_threads(n_threads, run_merge_job, create_merge_job, @@ -261,10 +265,23 @@ RefList *merge_intensities(Crystal **crystals, int n, int n_threads, refl != NULL; refl = next_refl(refl, iter) ) { + + /* Correct for averaging log of intensities*/ + if (ln_merge){ + double ln_I, ln_temp2; + + ln_temp2 = get_temp2(refl); + set_temp2(refl, exp(ln_temp2)); + + ln_I = get_intensity(refl); + set_intensity(refl, exp(ln_I)); + } + double var; int red; red = get_redundancy(refl); + //TODO is this still correct for log averaging? var = get_temp2(refl) / get_temp1(refl); set_esd_intensity(refl, sqrt(var)/sqrt(red)); @@ -276,6 +293,7 @@ RefList *merge_intensities(Crystal **crystals, int n, int n_threads, get_indices(refl, &h, &k, &l); r2 = add_refl(full2, h, k, l); copy_data(r2, refl); + } } diff --git a/src/merge.h b/src/merge.h index b33afdea..9bc54126 100644 --- a/src/merge.h +++ b/src/merge.h @@ -40,7 +40,8 @@ #include "geometry.h" extern RefList *merge_intensities(Crystal **crystals, int n, int n_threads, - int min_meas, double push_res, int use_weak); + int min_meas, double push_res, int use_weak, + int ln_merge); 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 e152db7c..6d10bc37 100644 --- a/src/partialator.c +++ b/src/partialator.c @@ -190,7 +190,7 @@ static void write_split(Crystal **crystals, int n_crystals, const char *outfile, } snprintf(tmp, 1024, "%s1", outfile); split = merge_intensities(crystals1, n_crystals1, nthreads, - min_measurements, push_res, 1); + min_measurements, push_res, 1, 0); if ( split == NULL ) { ERROR("Not enough crystals for two way split!\n"); @@ -202,7 +202,7 @@ static void write_split(Crystal **crystals, int n_crystals, const char *outfile, reflist_free(split); snprintf(tmp, 1024, "%s2", outfile); split = merge_intensities(crystals2, n_crystals2, nthreads, - min_measurements, push_res, 1); + min_measurements, push_res, 1, 0); STATUS("and %s\n", tmp); write_reflist_2(tmp, split, sym); reflist_free(split); @@ -289,7 +289,7 @@ static void write_custom_split(struct custom_split *csplit, int dsn, STATUS("Writing dataset '%s' to %s (%i crystals)\n", csplit->dataset_names[dsn], tmp, n_crystalsn); split = merge_intensities(crystalsn, n_crystalsn, nthreads, - min_measurements, push_res, 1); + min_measurements, push_res, 1, 0); write_reflist_2(tmp, split, sym); reflist_free(split); @@ -1343,7 +1343,7 @@ int main(int argc, char *argv[]) scale_all(crystals, n_crystals, nthreads, scaleflags); } full = merge_intensities(crystals, n_crystals, nthreads, - min_measurements, push_res, 1); + min_measurements, push_res, 1, 0); } else { full = reference; } @@ -1376,7 +1376,7 @@ int main(int argc, char *argv[]) } full = merge_intensities(crystals, n_crystals, nthreads, min_measurements, - push_res, 1); + push_res, 1, 0); } /* else full still equals reference */ check_rejection(crystals, n_crystals, full, max_B); @@ -1423,10 +1423,10 @@ int main(int argc, char *argv[]) } full = merge_intensities(crystals, n_crystals, nthreads, min_measurements, - push_res, 1); + push_res, 1, 0); } else { full = merge_intensities(crystals, n_crystals, nthreads, - min_measurements, push_res, 1); + min_measurements, push_res, 1, 0); } /* Write final figures of merit (no rejection any more) */ diff --git a/src/scaling.c b/src/scaling.c index 21bd0db2..dcfa7a59 100644 --- a/src/scaling.c +++ b/src/scaling.c @@ -144,7 +144,7 @@ void scale_all(Crystal **crystals, int n_crystals, int nthreads, int scaleflags) double bef_res; full = merge_intensities(crystals, n_crystals, nthreads, - 2, INFINITY, 0); + 2, INFINITY, 0, 1); old_res = new_res; bef_res = total_log_r(crystals, n_crystals, full, NULL); -- cgit v1.2.3 From 3563ec2b2dbc7d2090ed795a47e35a1b30b8b653 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Tue, 26 Jun 2018 14:07:25 +0200 Subject: Fix tests for log merging --- tests/partialator_merge_check_2 | 4 ++-- tests/partialator_merge_check_3 | 6 +++--- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/tests/partialator_merge_check_2 b/tests/partialator_merge_check_2 index a8ff3ddd..9938f64f 100755 --- a/tests/partialator_merge_check_2 +++ b/tests/partialator_merge_check_2 @@ -44,8 +44,8 @@ cat > partialator_merge_check_2_ans.hkl << EOF CrystFEL reflection list version 2.0 Symmetry: 1 h k l I phase sigma(I) nmeas - 1 0 0 150.00 - 0.00 2 - 19 0 0 75.00 - 0.00 2 + 1 0 0 141.42 - 0.00 2 + 19 0 0 70.71 - 0.00 2 EOF $PARTIALATOR -i partialator_merge_check_2.stream \ diff --git a/tests/partialator_merge_check_3 b/tests/partialator_merge_check_3 index 933422b0..2d237b02 100755 --- a/tests/partialator_merge_check_3 +++ b/tests/partialator_merge_check_3 @@ -45,9 +45,9 @@ cat > partialator_merge_check_3_ans.hkl << EOF CrystFEL reflection list version 2.0 Symmetry: 4 h k l I phase sigma(I) nmeas - 1 0 0 75.00 - 0.00 2 - 2 0 0 150.00 - 0.00 2 - 20 0 0 75.00 - 0.00 2 + 1 0 0 70.71 - 0.00 2 + 2 0 0 141.42 - 0.00 2 + 20 0 0 70.71 - 0.00 2 EOF $PARTIALATOR -i partialator_merge_check_3.stream \ -- cgit v1.2.3