diff options
Diffstat (limited to 'src')
| -rw-r--r-- | src/merge.c | 23 | ||||
| -rw-r--r-- | src/merge.h | 3 | ||||
| -rw-r--r-- | src/partialator.c | 14 | ||||
| -rw-r--r-- | src/scaling.c | 2 |
4 files changed, 30 insertions, 12 deletions
diff --git a/src/merge.c b/src/merge.c index 9d36a359..1dec13fd 100644 --- a/src/merge.c +++ b/src/merge.c @@ -63,6 +63,7 @@ struct merge_queue_args double push_res; int use_weak; long long int n_reflections; + int ln_merge; }; @@ -168,6 +169,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; @@ -193,7 +195,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; @@ -235,7 +237,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); @@ -272,7 +275,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; @@ -290,6 +293,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, @@ -305,10 +309,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)); 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 78b7173b..86431e4c 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"); @@ -203,7 +203,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); free_contribs(split); @@ -291,7 +291,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); free_contribs(split); reflist_free(split); @@ -1412,7 +1412,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; } @@ -1447,7 +1447,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, no_deltacchalf); @@ -1495,10 +1495,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 68456900..35481982 100644 --- a/src/scaling.c +++ b/src/scaling.c @@ -145,7 +145,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); |