From 3456372cb8a520ba187d3d8afe187c83c93a30b4 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Tue, 10 Mar 2015 13:54:37 +0100 Subject: Use running variance formula and add weighting --- src/hrs-scaling.c | 185 +++++++++--------------------------------------------- 1 file changed, 30 insertions(+), 155 deletions(-) diff --git a/src/hrs-scaling.c b/src/hrs-scaling.c index 52850fac..21b0c3fd 100644 --- a/src/hrs-scaling.c +++ b/src/hrs-scaling.c @@ -279,9 +279,8 @@ static void run_merge_job(void *vwargs, int cookie) { Reflection *f; signed int h, k, l; - double num, den; - int red; - double Ihl, corr, res; + double mean, sumweight, M2, temp, delta, R; + double corr, res, w, esd; if ( get_partiality(refl) < MIN_PART_MERGE ) continue; @@ -302,18 +301,15 @@ static void run_merge_job(void *vwargs, int cookie) f = add_refl(full, h, k, l); lock_reflection(f); pthread_rwlock_unlock(wargs->full_lock); - num = 0.0; - den = 0.0; - red = 0; + set_intensity(f, 0.0); + set_temp1(f, 0.0); + set_temp2(f, 0.0); } else { /* Someone else created it */ lock_reflection(f); pthread_rwlock_unlock(wargs->full_lock); - num = get_temp1(f); - den = get_temp2(f); - red = get_redundancy(f); } @@ -321,24 +317,30 @@ static void run_merge_job(void *vwargs, int cookie) lock_reflection(f); pthread_rwlock_unlock(wargs->full_lock); - num = get_temp1(f); - den = get_temp2(f); - red = get_redundancy(f); } - res = resolution(crystal_get_cell(cr), h, k, l); - corr = get_partiality(refl) * get_lorentz(refl); + mean = get_intensity(f); + sumweight = get_temp1(f); + M2 = get_temp2(f); - Ihl = get_intensity(refl) / corr; - - num += Ihl * G * exp(2.0*B*res); - den += 1.0; - red++; + res = resolution(crystal_get_cell(cr), h, k, l); - set_temp1(f, num); - set_temp2(f, den); - set_redundancy(f, red); + /* Total (multiplicative) correction factor */ + corr = G * exp(2.0*B*res) * get_lorentz(refl) + / get_partiality(refl); + + esd = get_esd_intensity(refl) * corr; + w = 1.0 / (esd*esd); + + /* Running mean and variance calculation */ + temp = w + sumweight; + delta = get_intensity(refl)*corr - mean; + R = delta * w / temp; + set_intensity(f, mean + R); + set_temp2(f, M2 + sumweight * delta * R); + set_temp1(f, temp); + set_redundancy(f, get_redundancy(f)+1); unlock_reflection(f); } } @@ -375,138 +377,15 @@ static RefList *lsq_intensities(Crystal **crystals, int n, int n_threads, refl != NULL; refl = next_refl(refl, iter) ) { - double Ih; - Ih = get_temp1(refl) / get_temp2(refl); - set_intensity(refl, Ih); - } - - return full; -} - - - -struct esd_queue_args -{ - RefList *full; - Crystal **crystals; - int n_started; - PartialityModel pmodel; -}; - - -struct esd_worker_args -{ - Crystal *crystal; - RefList *full; - PartialityModel pmodel; -}; - - -static void *create_esd_job(void *vqargs) -{ - struct esd_worker_args *wargs; - struct esd_queue_args *qargs = vqargs; - - wargs = malloc(sizeof(struct esd_worker_args)); - wargs->full = qargs->full; - wargs->pmodel = qargs->pmodel; - - wargs->crystal = qargs->crystals[qargs->n_started++]; - - return wargs; -} - - -static void run_esd_job(void *vwargs, int cookie) -{ - struct esd_worker_args *wargs = vwargs; - Crystal *cr = wargs->crystal; - RefList *full = wargs->full; - Reflection *refl; - RefListIterator *iter; - double G; - - /* If this crystal's scaling was dodgy, it doesn't contribute to the - * merged intensities */ - if ( crystal_get_user_flag(cr) != 0 ) return; - - G = crystal_get_osf(cr); - - for ( refl = first_refl(crystal_get_reflections(cr), &iter); - refl != NULL; - refl = next_refl(refl, iter) ) - { - Reflection *f; - signed int h, k, l; - double num; - double Ihl, Ih, corr; - - if ( get_partiality(refl) < MIN_PART_MERGE ) continue; - - get_indices(refl, &h, &k, &l); - f = find_refl(full, h, k, l); - assert(f != NULL); - - lock_reflection(f); - - num = get_temp1(f); - - corr = get_partiality(refl) * get_lorentz(refl); - - Ih = get_intensity(f); - Ihl = get_intensity(refl) / (G*corr); - - num += pow(Ihl - Ih, 2.0); - - set_temp1(f, num); - unlock_reflection(f); - } -} - - -static void finalise_esd_job(void *vqargs, void *vwargs) -{ - free(vwargs); -} - - -static void calculate_esds(Crystal **crystals, int n, RefList *full, - int n_threads, int min_red, PartialityModel pmodel) -{ - struct esd_queue_args qargs; - Reflection *refl; - RefListIterator *iter; - - qargs.full = full; - qargs.n_started = 0; - qargs.crystals = crystals; - qargs.pmodel = pmodel; + double var; + int red; - for ( refl = first_refl(full, &iter); - refl != NULL; - refl = next_refl(refl, iter) ) - { - set_temp1(refl, 0.0); - set_temp2(refl, 0.0); + red = get_redundancy(refl); + var = get_temp2(refl) / get_temp1(refl); + set_esd_intensity(refl, sqrt(var)/sqrt(red)); } - run_threads(n_threads, run_esd_job, create_esd_job, - finalise_esd_job, &qargs, n, 0, 0, 0); - - for ( refl = first_refl(full, &iter); - refl != NULL; - refl = next_refl(refl, iter) ) - { - double esd; - int red = get_redundancy(refl); - esd = sqrt(get_temp1(refl)); - esd /= (double)red; - set_esd_intensity(refl, esd); - - if ( red < min_red ) { - set_redundancy(refl, 0); - } - } + return full; } @@ -586,8 +465,6 @@ RefList *scale_intensities(Crystal **crystals, int n, if ( noscale ) { full = lsq_intensities(crystals, n, n_threads, pmodel); - calculate_esds(crystals, n, full, n_threads, min_redundancy, - pmodel); return full; } @@ -642,8 +519,6 @@ RefList *scale_intensities(Crystal **crystals, int n, ERROR("WARNING: Scaling did not converge.\n"); } - calculate_esds(crystals, n, full, n_threads, min_redundancy, pmodel); - free(old_osfs); return full; } -- cgit v1.2.3