diff options
Diffstat (limited to 'src')
-rw-r--r-- | src/merge.c | 60 | ||||
-rw-r--r-- | src/partialator.c | 13 | ||||
-rw-r--r-- | src/post-refinement.c | 4 | ||||
-rw-r--r-- | src/post-refinement.h | 2 | ||||
-rw-r--r-- | src/rejection.c | 210 | ||||
-rw-r--r-- | src/rejection.h | 6 | ||||
-rw-r--r-- | src/scaling.c | 2 |
7 files changed, 232 insertions, 65 deletions
diff --git a/src/merge.c b/src/merge.c index 14460953..9d36a359 100644 --- a/src/merge.c +++ b/src/merge.c @@ -3,11 +3,11 @@ * * Parallel weighted merging of intensities * - * Copyright © 2012-2015 Deutsches Elektronen-Synchrotron DESY, + * Copyright © 2012-2018 Deutsches Elektronen-Synchrotron DESY, * a research centre of the Helmholtz Association. * * Authors: - * 2010-2015 Thomas White <taw@physics.org> + * 2010-2018 Thomas White <taw@physics.org> * * This file is part of CrystFEL. * @@ -46,6 +46,7 @@ #include "cell.h" #include "utils.h" #include "reflist.h" +#include "reflist-utils.h" #include "cell-utils.h" @@ -88,6 +89,17 @@ static void *create_merge_job(void *vqargs) } +static int alloc_contribs(struct reflection_contributions *c) +{ + c->contribs = realloc(c->contribs, c->max_contrib*sizeof(Reflection *)); + c->contrib_crystals = realloc(c->contrib_crystals, + c->max_contrib*sizeof(Crystal *)); + if ( c->contribs == NULL ) return 1; + if ( c->contrib_crystals == NULL ) return 1; + return 0; +} + + /* Find reflection hkl in 'list', creating it if it's not there, under * protection of 'lock' and returning a locked reflection */ static Reflection *get_locked_reflection(RefList *list, pthread_rwlock_t *lock, @@ -108,12 +120,31 @@ static Reflection *get_locked_reflection(RefList *list, pthread_rwlock_t *lock, * So, we must check again */ f = find_refl(list, h, k, l); if ( f == NULL ) { + + struct reflection_contributions *c; + f = add_refl(list, h, k, l); lock_reflection(f); pthread_rwlock_unlock(lock); set_intensity(f, 0.0); set_temp1(f, 0.0); set_temp2(f, 0.0); + + c = malloc(sizeof(struct reflection_contributions)); + if ( c != NULL ) { + c->n_contrib = 0; + c->max_contrib = 32; + c->contribs = NULL; + c->contrib_crystals = NULL; + if ( alloc_contribs(c) ) { + set_contributions(f, NULL); + } else { + set_contributions(f, c); + } + } else { + set_contributions(f, NULL); + } + } else { /* Someone else created it */ lock_reflection(f); @@ -158,6 +189,7 @@ static void run_merge_job(void *vwargs, int cookie) signed int h, k, l; double mean, sumweight, M2, temp, delta, R; double corr, res, w; + struct reflection_contributions *c; if ( get_partiality(refl) < MIN_PART_MERGE ) continue; @@ -209,6 +241,18 @@ static void run_merge_job(void *vwargs, int cookie) set_temp2(f, M2 + sumweight * delta * R); set_temp1(f, temp); set_redundancy(f, get_redundancy(f)+1); + + /* Record this contribution */ + c = get_contributions(f); + if ( c != NULL ) { + c->contribs[c->n_contrib] = refl; + c->contrib_crystals[c->n_contrib++] = cr; + if ( c->n_contrib == c->max_contrib ) { + c->max_contrib += 64; + alloc_contribs(c); + } + } /* else, too bad! */ + unlock_reflection(f); wargs->n_reflections++; @@ -274,8 +318,18 @@ RefList *merge_intensities(Crystal **crystals, int n, int n_threads, Reflection *r2; get_indices(refl, &h, &k, &l); - r2 = add_refl(full2, h, k, l); + r2 = add_refl(full2, h, k, l); copy_data(r2, refl); + + } else { + + /* We do not need the contribution list any more */ + struct reflection_contributions *c; + c = get_contributions(refl); + free(c->contribs); + free(c->contrib_crystals); + free(c); + } } diff --git a/src/partialator.c b/src/partialator.c index 113535e0..78b7173b 100644 --- a/src/partialator.c +++ b/src/partialator.c @@ -199,12 +199,14 @@ static void write_split(Crystal **crystals, int n_crystals, const char *outfile, STATUS("Writing two-way split to %s ", tmp); write_reflist_2(tmp, split, sym); + free_contribs(split); reflist_free(split); snprintf(tmp, 1024, "%s2", outfile); split = merge_intensities(crystals2, n_crystals2, nthreads, min_measurements, push_res, 1); STATUS("and %s\n", tmp); write_reflist_2(tmp, split, sym); + free_contribs(split); reflist_free(split); } @@ -291,6 +293,7 @@ static void write_custom_split(struct custom_split *csplit, int dsn, split = merge_intensities(crystalsn, n_crystalsn, nthreads, min_measurements, push_res, 1); write_reflist_2(tmp, split, sym); + free_contribs(split); reflist_free(split); write_split(crystalsn, n_crystalsn, tmp, nthreads, pmodel, @@ -318,6 +321,7 @@ static void show_help(const char *s) " --no-scale Disable scale factor (G, B) refinement.\n" " --no-Bscale Disable B factor scaling.\n" " --no-pr Disable orientation/physics refinement.\n" +" --no-deltacchalf Disable rejection based on deltaCChalf.\n" " -m, --model=<model> Specify partiality model.\n" " --min-measurements=<n> Minimum number of measurements to require.\n" " --no-polarisation Disable polarisation correction.\n" @@ -912,6 +916,7 @@ int main(int argc, char *argv[]) int scaleflags = 0; double min_res = 0.0; int do_write_logs = 0; + int no_deltacchalf = 0; /* Long options */ const struct option longopts[] = { @@ -950,6 +955,7 @@ int main(int argc, char *argv[]) {"no-free", 0, &no_free, 1}, {"output-every-cycle", 0, &output_everycycle, 1}, {"no-logs", 0, &no_logs, 1}, + {"no-deltacchalf", 0, &no_deltacchalf, 1}, {0, 0, NULL, 0} }; @@ -1412,7 +1418,7 @@ int main(int argc, char *argv[]) } /* Check rejection and write figures of merit */ - check_rejection(crystals, n_crystals, full, max_B); + check_rejection(crystals, n_crystals, full, max_B, no_deltacchalf); show_all_residuals(crystals, n_crystals, full); if ( do_write_logs ) { @@ -1433,6 +1439,7 @@ int main(int argc, char *argv[]) /* Create new reference if needed */ if ( reference == NULL ) { + free_contribs(full); reflist_free(full); if ( !no_scale ) { scale_all(crystals, n_crystals, nthreads, @@ -1443,7 +1450,7 @@ int main(int argc, char *argv[]) push_res, 1); } /* else full still equals reference */ - check_rejection(crystals, n_crystals, full, max_B); + check_rejection(crystals, n_crystals, full, max_B, no_deltacchalf); show_all_residuals(crystals, n_crystals, full); if ( do_write_logs ) { @@ -1481,6 +1488,7 @@ int main(int argc, char *argv[]) /* Final merge */ STATUS("Final merge...\n"); if ( reference == NULL ) { + free_contribs(full); reflist_free(full); if ( !no_scale ) { scale_all(crystals, n_crystals, nthreads, scaleflags); @@ -1527,6 +1535,7 @@ int main(int argc, char *argv[]) reflist_free(crystal_get_reflections(crystals[i])); crystal_free(crystals[i]); } + free_contribs(full); reflist_free(full); free_symoplist(sym); free(outfile); diff --git a/src/post-refinement.c b/src/post-refinement.c index bde9c02c..0a799b9e 100644 --- a/src/post-refinement.c +++ b/src/post-refinement.c @@ -69,8 +69,8 @@ const char *str_prflag(enum prflag flag) case PRFLAG_EARLY : return "early rejection"; - case PRFLAG_CC : - return "low CC"; + case PRFLAG_DELTACCHALF : + return "negative delta CC½"; case PRFLAG_BIGB : return "B too big"; diff --git a/src/post-refinement.h b/src/post-refinement.h index 71a6d7f3..e2873b76 100644 --- a/src/post-refinement.h +++ b/src/post-refinement.h @@ -50,7 +50,7 @@ enum prflag PRFLAG_FEWREFL = 16, PRFLAG_SOLVEFAIL = 17, PRFLAG_EARLY = 18, - PRFLAG_CC = 19, + PRFLAG_DELTACCHALF = 19, PRFLAG_BIGB = 20, PRFLAG_SCALEBAD = 21, }; diff --git a/src/rejection.c b/src/rejection.c index 9d41f9b4..230093e6 100644 --- a/src/rejection.c +++ b/src/rejection.c @@ -3,11 +3,11 @@ * * Crystal rejection for scaling * - * Copyright © 2012-2015 Deutsches Elektronen-Synchrotron DESY, + * Copyright © 2012-2018 Deutsches Elektronen-Synchrotron DESY, * a research centre of the Helmholtz Association. * * Authors: - * 2010-2015 Thomas White <taw@physics.org> + * 2010-2018 Thomas White <taw@physics.org> * * This file is part of CrystFEL. * @@ -94,74 +94,174 @@ void early_rejection(Crystal **crystals, int n) } -static void check_cc(Crystal *cr, RefList *full) +static double calculate_cchalf(RefList *template, RefList *full, + Crystal *exclude, int *pnref) { - RefList *list = crystal_get_reflections(cr); - Reflection *refl; + Reflection *trefl; RefListIterator *iter; - double G = crystal_get_osf(cr); - double B = crystal_get_Bfac(cr); - UnitCell *cell = crystal_get_cell(cr); - double *vec1, *vec2; - int n, i; - double cc; - - n = num_reflections(list); - vec1 = malloc(n*sizeof(double)); - vec2 = malloc(n*sizeof(double)); - if ( (vec1 == NULL) || (vec2 == NULL) ) { - ERROR("Not enough memory to check CCs\n"); - return; - } - - i = 0; - for ( refl = first_refl(list, &iter); - refl != NULL; - refl = next_refl(refl, iter) ) + double sig2E, sig2Y; + int n = 0; + signed int oh = 0; + signed int ok = 0; + signed int ol = 0; + double wSum = 0.0; + double wSum2 = 0.0; + double mean = 0.0; + double S = 0.0; + double all_sum_var = 0.0; + long int total_contribs = 0; + + /* Iterate over all reflections */ + for ( trefl = first_refl(template, &iter); + trefl != NULL; + trefl = next_refl(trefl, iter) ) { + struct reflection_contributions *c; + int j; + double refl_sum; + double refl_sumsq; + double refl_mean, refl_var; signed int h, k, l; - double pobs, pcalc; - double res, corr, Ipart; - Reflection *match; + int nc = 0; + Reflection *refl; + double res; + + get_indices(trefl, &h, &k, &l); + + /* next_refl() will iterate over multiple copies of the same + * unique reflection, but we are only interested in seeing each + * unique reflection once. */ + if ( (h==oh) && (k==ok) && (l==ol) ) continue; + oh = h; ok = k; ol = l; + + refl = find_refl(full, h, k, l); + if ( refl == NULL ) continue; + + c = get_contributions(refl); + assert(c != NULL); - if ( !get_flag(refl) ) continue; /* Not free-flagged */ + /* Calculate the resolution just once, using the cell from the + * first crystal to contribute, otherwise it takes too long */ + res = resolution(crystal_get_cell(c->contrib_crystals[0]), + h, k, l); - get_indices(refl, &h, &k, &l); - res = resolution(cell, h, k, l); - if ( 2.0*res > crystal_get_resolution_limit(cr) ) continue; + /* Mean of contributions */ + refl_sum = 0.0; + for ( j=0; j<c->n_contrib; j++ ) { - match = find_refl(full, h, k, l); - if ( match == NULL ) continue; + double Ii, G, B; - /* Calculated partiality */ - pcalc = get_partiality(refl); + if ( c->contrib_crystals[j] == exclude ) { + continue; + } + + G = crystal_get_osf(c->contrib_crystals[j]); + B = crystal_get_Bfac(c->contrib_crystals[j]); + Ii = get_intensity(c->contribs[j]); + + /* Total (multiplicative) correction factor */ + Ii *= 1.0/G * exp(B*res*res) * get_lorentz(c->contribs[j]) + / get_partiality(c->contribs[j]); + + refl_sum += Ii; + nc++; + + } + + if ( nc < 2 ) continue; + refl_mean = refl_sum / nc; + + /* Variance of contributions */ + refl_sumsq = 0.0; + for ( j=0; j<c->n_contrib; j++ ) { + + double Ii, G, B; + + if ( c->contrib_crystals[j] == exclude ) { + continue; + } - /* Observed partiality */ - corr = G * exp(B*res*res) * get_lorentz(refl); - Ipart = get_intensity(refl) * corr; - pobs = Ipart / get_intensity(match); + G = crystal_get_osf(c->contrib_crystals[j]); + B = crystal_get_Bfac(c->contrib_crystals[j]); + Ii = get_intensity(c->contribs[j]); + + /* Total (multiplicative) correction factor */ + Ii *= 1.0/G * exp(B*res*res) * get_lorentz(c->contribs[j]) + / get_partiality(c->contribs[j]); + + refl_sumsq += (Ii-refl_mean)*(Ii-refl_mean); + /* (nc already summed above) */ + + } + + refl_var = refl_sumsq/(nc-1.0); + refl_var /= (nc/2.0); + + all_sum_var += refl_var; + n++; + total_contribs += nc; + + double w = 1.0; + wSum += w; + wSum2 += w*w; + double meanOld = mean; + mean = meanOld + (w/wSum) * (refl_mean - meanOld); + S += w * (refl_mean - meanOld) * (refl_mean - mean); - vec1[i] = pobs; - vec2[i] = pcalc; - i++; } - cc = gsl_stats_correlation(vec1, 1, vec2, 1, i); - //printf("%f\n", cc); - if ( cc < 0.5 ) crystal_set_user_flag(cr, PRFLAG_CC); + sig2E = all_sum_var / n; + sig2Y = S / (wSum - 1.0); - free(vec1); - free(vec2); + if ( pnref != NULL ) { + *pnref = n; + } + return (sig2Y - 0.5*sig2E) / (sig2Y + 0.5*sig2E); } -static UNUSED void check_ccs(Crystal **crystals, int n_crystals, RefList *full) +static void check_deltacchalf(Crystal **crystals, int n, RefList *full) { + double cchalf; int i; + int nref; + double *vals; + double mean, sd; + + cchalf = calculate_cchalf(full, full, NULL, &nref); + STATUS("Overall CChalf = %f %% (%i reflections)\n", cchalf*100.0, nref); + + vals = malloc(n*sizeof(double)); + if ( vals == NULL ) { + ERROR("Not enough memory for deltaCChalf check\n"); + return; + } - for ( i=0; i<n_crystals; i++ ) { - check_cc(crystals[i], full); + for ( i=0; i<n; i++ ) { + double cchalf, cchalfi; + RefList *template = crystal_get_reflections(crystals[i]); + cchalf = calculate_cchalf(template, full, NULL, NULL); + cchalfi = calculate_cchalf(template, full, crystals[i], &nref); + //STATUS("Frame %i:", i); + //STATUS(" With = %f ", cchalf*100.0); + //STATUS("Without = %f", cchalfi*100.0); + //STATUS(" Delta = %f ", (cchalf - cchalfi)*100.0); + //STATUS("(nref = %i)\n", nref); + vals[i] = cchalf - cchalfi; + progress_bar(i, n-1, "Calculating deltaCChalf"); } + + mean = gsl_stats_mean(vals, 1, n); + sd = gsl_stats_sd_m(vals, 1, n, mean); + STATUS("deltaCChalf = %f ± %f %%\n", mean*100.0, sd*100.0); + + for ( i=0; i<n; i++ ) { + if ( vals[i] < mean-2.0*sd ) { + crystal_set_user_flag(crystals[i], PRFLAG_DELTACCHALF); + } + } + + free(vals); } @@ -192,14 +292,16 @@ static void show_duds(Crystal **crystals, int n_crystals) } -void check_rejection(Crystal **crystals, int n, RefList *full, double max_B) +void check_rejection(Crystal **crystals, int n, RefList *full, double max_B, + int no_deltacchalf) { int i; int n_acc = 0; - /* Check according to CCs FIXME: Disabled */ - //if ( full != NULL ) check_ccs(crystals, n, full); - + /* Check according to delta CC½ */ + if ( !no_deltacchalf && (full != NULL) ) { + check_deltacchalf(crystals, n, full); + } for ( i=0; i<n; i++ ) { if ( fabs(crystal_get_Bfac(crystals[i])) > max_B ) { diff --git a/src/rejection.h b/src/rejection.h index 8979313b..5a0c3e9f 100644 --- a/src/rejection.h +++ b/src/rejection.h @@ -3,11 +3,11 @@ * * Crystal rejection for scaling * - * Copyright © 2012-2015 Deutsches Elektronen-Synchrotron DESY, + * Copyright © 2012-2018 Deutsches Elektronen-Synchrotron DESY, * a research centre of the Helmholtz Association. * * Authors: - * 2010-2015 Thomas White <taw@physics.org> + * 2010-2018 Thomas White <taw@physics.org> * * This file is part of CrystFEL. * @@ -39,6 +39,6 @@ extern void early_rejection(Crystal **crystals, int n); extern void check_rejection(Crystal **crystals, int n, RefList *full, - double max_B); + double max_B, int no_deltacchalf); #endif /* REJECTION_H */ diff --git a/src/scaling.c b/src/scaling.c index 21bd0db2..68456900 100644 --- a/src/scaling.c +++ b/src/scaling.c @@ -46,6 +46,7 @@ #include "cell.h" #include "cell-utils.h" #include "scaling.h" +#include "reflist-utils.h" struct scale_args @@ -166,6 +167,7 @@ void scale_all(Crystal **crystals, int n_crystals, int nthreads, int scaleflags) meanB /= n_crystals; STATUS("Mean B = %e\n", meanB); + free_contribs(full); reflist_free(full); niter++; |