aboutsummaryrefslogtreecommitdiff
path: root/src
diff options
context:
space:
mode:
Diffstat (limited to 'src')
-rw-r--r--src/merge.c24
-rw-r--r--src/merge.h3
-rw-r--r--src/partialator.c14
-rw-r--r--src/scaling.c2
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);
generated by cgit v1.2.3 (git 2.25.1) at 2024-11-05 18:15:50 +0000