aboutsummaryrefslogtreecommitdiff
path: root/src/hrs-scaling.c
diff options
context:
space:
mode:
authorThomas White <taw@physics.org>2015-03-10 13:54:37 +0100
committerThomas White <taw@physics.org>2015-03-12 16:37:25 +0100
commit3456372cb8a520ba187d3d8afe187c83c93a30b4 (patch)
treea0c56698906e60a1de24ac5f3364bf5e9159a346 /src/hrs-scaling.c
parent5e34b26a7965eafbec680148367697a371de35e8 (diff)
Use running variance formula and add weighting
Diffstat (limited to 'src/hrs-scaling.c')
-rw-r--r--src/hrs-scaling.c185
1 files 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;
}