From 609dcac27949236c7a1acda0bd3e4342cdb328f6 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Mon, 20 Aug 2018 16:42:25 +0200 Subject: Initial delta CC half stuff --- src/post-refinement.c | 2 +- src/rejection.c | 211 +++++++++++++++++++++++++++++++++++++++----------- 2 files changed, 166 insertions(+), 47 deletions(-) (limited to 'src') diff --git a/src/post-refinement.c b/src/post-refinement.c index bde9c02c..ebaf66ab 100644 --- a/src/post-refinement.c +++ b/src/post-refinement.c @@ -70,7 +70,7 @@ const char *str_prflag(enum prflag flag) return "early rejection"; case PRFLAG_CC : - return "low CC"; + return "negative delta CC½"; case PRFLAG_BIGB : return "B too big"; diff --git a/src/rejection.c b/src/rejection.c index 9d41f9b4..716baa1b 100644 --- a/src/rejection.c +++ b/src/rejection.c @@ -94,74 +94,194 @@ void early_rejection(Crystal **crystals, int n) } -static void check_cc(Crystal *cr, RefList *full) +struct contribs +{ + int serial; /* h,k,l */ + int n_contrib; + int max_contrib; + Reflection **contribs; + Crystal **contrib_crystals; +}; + + +struct contributionlist +{ + struct contribs *contribs; + int n; /* Number of reflections */ +}; + + +static int alloc_contribs(struct contribs *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; +} + + +static void free_contributions(struct contributionlist *clist) +{ + int i; + for ( i=0; in; i++ ) { + free(clist->contribs[i].contribs); + free(clist->contribs[i].contrib_crystals); + } + free(clist->contribs); + free(clist); +} + + +static struct contributionlist *find_all_contributions(Crystal **crystals, + int n, RefList *full) { - RefList *list = crystal_get_reflections(cr); Reflection *refl; 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; + struct contributionlist *clist; + int nc = 0; + + clist = malloc(sizeof(struct contributionlist)); + if ( clist == NULL ) return NULL; + + clist->n = num_reflections(full); + clist->contribs = malloc(clist->n*sizeof(struct contribs)); + if ( clist->contribs == NULL ) { + free(clist); + return NULL; } - i = 0; - for ( refl = first_refl(list, &iter); + for ( refl = first_refl(full, &iter); refl != NULL; refl = next_refl(refl, iter) ) { + int i; signed int h, k, l; - double pobs, pcalc; - double res, corr, Ipart; - Reflection *match; - - if ( !get_flag(refl) ) continue; /* Not free-flagged */ + struct contribs *c; get_indices(refl, &h, &k, &l); - res = resolution(cell, h, k, l); - if ( 2.0*res > crystal_get_resolution_limit(cr) ) continue; + c = &clist->contribs[nc++]; + + c->serial = SERIAL(h, k, l); + c->n_contrib = 0; + c->max_contrib = 32; + c->contribs = NULL; + c->contrib_crystals = NULL; + if ( alloc_contribs(c) ) return NULL; - match = find_refl(full, h, k, l); - if ( match == NULL ) continue; + /* Find all versions of this reflection */ + for ( i=0; icontribs[c->n_contrib] = refl2; + c->contrib_crystals[c->n_contrib] = crystals[i]; + c->n_contrib++; - vec1[i] = pobs; - vec2[i] = pcalc; - i++; + if ( c->n_contrib == c->max_contrib ) { + c->max_contrib += 64; + if ( alloc_contribs(c) ) return NULL; + } + + refl2 = next_found_refl(refl2); + + } while ( refl2 != NULL ); + } + + progress_bar(nc, clist->n, "Collecting contributions"); } - cc = gsl_stats_correlation(vec1, 1, vec2, 1, i); - //printf("%f\n", cc); - if ( cc < 0.5 ) crystal_set_user_flag(cr, PRFLAG_CC); + return clist; +} + - free(vec1); - free(vec2); +static struct contribs *lookup_contribs(struct contributionlist *clist, + signed int h, signed int k, signed int l) +{ + int i, serial; + + serial = SERIAL(h, k, l); + + for ( i=0; in; i++ ) { + if ( clist->contribs[i].serial == serial ) { + return &clist->contribs[i]; + } + } + return NULL; } -static UNUSED void check_ccs(Crystal **crystals, int n_crystals, RefList *full) +static double calculate_cchalf(Crystal **crystals, int n, + struct contributionlist *clist, Crystal *exclude) { int i; + double acc_var = 0.0; + double acc_mean = 0.0; + int acc_n = 0; + double sigE, sigy; + + for ( i=0; in_contrib; j++ ) { + + Reflection *reflc; + double Ii; + + if ( c->contrib_crystals[j] == exclude ) { + continue; + } + + reflc = c->contribs[j]; + + /* FIXME: Apply corrections */ + Ii = get_intensity(reflc); + acc_mean += Ii; + acc_var += Ii * Ii; + acc_n++; + } + + } - for ( i=0; i max_B ) { -- cgit v1.2.3 From 5447e68a8885e565d53a30fc1c1ad6e358505a8b Mon Sep 17 00:00:00 2001 From: Thomas White Date: Mon, 20 Aug 2018 17:13:51 +0200 Subject: Working overall CChalf --- src/rejection.c | 71 +++++++++++++++++++++++++++------------------------------ 1 file changed, 34 insertions(+), 37 deletions(-) (limited to 'src') diff --git a/src/rejection.c b/src/rejection.c index 716baa1b..eaf87fb6 100644 --- a/src/rejection.c +++ b/src/rejection.c @@ -220,66 +220,63 @@ static double calculate_cchalf(Crystal **crystals, int n, struct contributionlist *clist, Crystal *exclude) { int i; - double acc_var = 0.0; - double acc_mean = 0.0; - int acc_n = 0; - double sigE, sigy; + double all_sum_mean = 0.0; + double all_sumsq_mean = 0.0; + double all_sum_var = 0.0; + double sig2E, sig2Y; - for ( i=0; in; i++ ) { - if ( cr == exclude ) continue; + struct contribs *c; + int j; + double refl_sum = 0.0; + double refl_sumsq = 0.0; + double refl_mean, refl_var; - for ( refl = first_refl(list, &iter); - refl != NULL; - refl = next_refl(refl, iter) ) - { - signed int h, k, l; - struct contribs *c; - int j; + c = &clist->contribs[i]; - get_indices(refl, &h, &k, &l); - c = lookup_contribs(clist, h, k, l); - assert(c != NULL); - for ( j=0; jn_contrib; j++ ) { + for ( j=0; jn_contrib; j++ ) { - Reflection *reflc; - double Ii; + double Ii; - if ( c->contrib_crystals[j] == exclude ) { - continue; - } + if ( c->contrib_crystals[j] == exclude ) { + continue; + } - reflc = c->contribs[j]; + /* FIXME: Apply corrections */ + Ii = get_intensity(c->contribs[j]); - /* FIXME: Apply corrections */ - Ii = get_intensity(reflc); - acc_mean += Ii; - acc_var += Ii * Ii; - acc_n++; - } + refl_sum += Ii; + refl_sumsq += Ii * Ii; } + refl_mean = refl_sum / c->n_contrib; + refl_var = refl_sumsq - refl_sum*refl_sum; + + all_sum_mean += refl_mean; + all_sumsq_mean += refl_mean*refl_mean; + all_sum_var += refl_var; + } + sig2E = all_sum_var / clist->n; + sig2Y = all_sumsq_mean - all_sum_mean*all_sum_mean; - return 0.0; /* FIXME */ + return (sig2Y - 0.5*sig2E) / (sig2Y + 0.5*sig2E); } static void check_deltacchalf(Crystal **crystals, int n, RefList *full) { struct contributionlist *clist; + double cchalf; clist = find_all_contributions(crystals, n, full); - STATUS("Overall CChalf = %f\n", - calculate_cchalf(crystals, n, clist, NULL)); + cchalf = calculate_cchalf(crystals, n, clist, NULL); + STATUS("Overall CChalf = %f\n", cchalf); free_contributions(clist); } -- cgit v1.2.3 From 214b23a62b478acddaeb8d52f6d59909fa931c1e Mon Sep 17 00:00:00 2001 From: Thomas White Date: Mon, 20 Aug 2018 17:37:08 +0200 Subject: Calculate the DeltaCChalf --- src/rejection.c | 12 +++++++++--- 1 file changed, 9 insertions(+), 3 deletions(-) (limited to 'src') diff --git a/src/rejection.c b/src/rejection.c index eaf87fb6..8cd4f718 100644 --- a/src/rejection.c +++ b/src/rejection.c @@ -216,8 +216,7 @@ static struct contribs *lookup_contribs(struct contributionlist *clist, } -static double calculate_cchalf(Crystal **crystals, int n, - struct contributionlist *clist, Crystal *exclude) +static double calculate_cchalf(struct contributionlist *clist, Crystal *exclude) { int i; double all_sum_mean = 0.0; @@ -272,12 +271,19 @@ static void check_deltacchalf(Crystal **crystals, int n, RefList *full) { struct contributionlist *clist; double cchalf; + int i; clist = find_all_contributions(crystals, n, full); - cchalf = calculate_cchalf(crystals, n, clist, NULL); + cchalf = calculate_cchalf(clist, NULL); STATUS("Overall CChalf = %f\n", cchalf); + for ( i=0; i Date: Mon, 20 Aug 2018 18:33:05 +0200 Subject: Save contribution information during merge --- src/merge.c | 43 +++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 43 insertions(+) (limited to 'src') diff --git a/src/merge.c b/src/merge.c index 14460953..9e577224 100644 --- a/src/merge.c +++ b/src/merge.c @@ -88,6 +88,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 +119,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 +188,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 +240,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++; -- cgit v1.2.3 From 0fcf1d67fdb3e4238e1c5d2e17ff23c466850224 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Mon, 20 Aug 2018 18:33:24 +0200 Subject: Calculate deltaCChalf using previously stored information --- src/rejection.c | 147 +++++--------------------------------------------------- 1 file changed, 13 insertions(+), 134 deletions(-) (limited to 'src') diff --git a/src/rejection.c b/src/rejection.c index 8cd4f718..53902c2d 100644 --- a/src/rejection.c +++ b/src/rejection.c @@ -94,146 +94,29 @@ void early_rejection(Crystal **crystals, int n) } -struct contribs -{ - int serial; /* h,k,l */ - int n_contrib; - int max_contrib; - Reflection **contribs; - Crystal **contrib_crystals; -}; - - -struct contributionlist -{ - struct contribs *contribs; - int n; /* Number of reflections */ -}; - - -static int alloc_contribs(struct contribs *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; -} - - -static void free_contributions(struct contributionlist *clist) -{ - int i; - for ( i=0; in; i++ ) { - free(clist->contribs[i].contribs); - free(clist->contribs[i].contrib_crystals); - } - free(clist->contribs); - free(clist); -} - - -static struct contributionlist *find_all_contributions(Crystal **crystals, - int n, RefList *full) +static double calculate_cchalf(RefList *full, Crystal *exclude) { Reflection *refl; RefListIterator *iter; - struct contributionlist *clist; - int nc = 0; - - clist = malloc(sizeof(struct contributionlist)); - if ( clist == NULL ) return NULL; - - clist->n = num_reflections(full); - clist->contribs = malloc(clist->n*sizeof(struct contribs)); - if ( clist->contribs == NULL ) { - free(clist); - return NULL; - } - - for ( refl = first_refl(full, &iter); - refl != NULL; - refl = next_refl(refl, iter) ) - { - int i; - signed int h, k, l; - struct contribs *c; - - get_indices(refl, &h, &k, &l); - c = &clist->contribs[nc++]; - - c->serial = SERIAL(h, k, l); - c->n_contrib = 0; - c->max_contrib = 32; - c->contribs = NULL; - c->contrib_crystals = NULL; - if ( alloc_contribs(c) ) return NULL; - - /* Find all versions of this reflection */ - for ( i=0; icontribs[c->n_contrib] = refl2; - c->contrib_crystals[c->n_contrib] = crystals[i]; - c->n_contrib++; - - if ( c->n_contrib == c->max_contrib ) { - c->max_contrib += 64; - if ( alloc_contribs(c) ) return NULL; - } - - refl2 = next_found_refl(refl2); - - } while ( refl2 != NULL ); - } - - progress_bar(nc, clist->n, "Collecting contributions"); - } - - return clist; -} - - -static struct contribs *lookup_contribs(struct contributionlist *clist, - signed int h, signed int k, signed int l) -{ - int i, serial; - - serial = SERIAL(h, k, l); - - for ( i=0; in; i++ ) { - if ( clist->contribs[i].serial == serial ) { - return &clist->contribs[i]; - } - } - return NULL; -} - - -static double calculate_cchalf(struct contributionlist *clist, Crystal *exclude) -{ - int i; double all_sum_mean = 0.0; double all_sumsq_mean = 0.0; double all_sum_var = 0.0; double sig2E, sig2Y; + int n = 0; /* Iterate over all reflections */ - for ( i=0; in; i++ ) { - - struct contribs *c; + for ( refl = first_refl(full, &iter); + refl != NULL; + refl = next_refl(refl, iter) ) + { + struct reflection_contributions *c; int j; double refl_sum = 0.0; double refl_sumsq = 0.0; double refl_mean, refl_var; - c = &clist->contribs[i]; + c = get_contributions(refl); + assert(c != NULL); for ( j=0; jn_contrib; j++ ) { @@ -257,10 +140,11 @@ static double calculate_cchalf(struct contributionlist *clist, Crystal *exclude) all_sum_mean += refl_mean; all_sumsq_mean += refl_mean*refl_mean; all_sum_var += refl_var; + n++; } - sig2E = all_sum_var / clist->n; + sig2E = all_sum_var / n; sig2Y = all_sumsq_mean - all_sum_mean*all_sum_mean; return (sig2Y - 0.5*sig2E) / (sig2Y + 0.5*sig2E); @@ -269,22 +153,17 @@ static double calculate_cchalf(struct contributionlist *clist, Crystal *exclude) static void check_deltacchalf(Crystal **crystals, int n, RefList *full) { - struct contributionlist *clist; double cchalf; int i; - clist = find_all_contributions(crystals, n, full); - - cchalf = calculate_cchalf(clist, NULL); + cchalf = calculate_cchalf(full, NULL); STATUS("Overall CChalf = %f\n", cchalf); for ( i=0; i Date: Tue, 21 Aug 2018 17:56:32 +0200 Subject: Calculate DeltaCChalf using only the reflections in the crystal under consideration --- src/rejection.c | 25 +++++++++++++++++-------- 1 file changed, 17 insertions(+), 8 deletions(-) (limited to 'src') diff --git a/src/rejection.c b/src/rejection.c index 53902c2d..46e7a2f2 100644 --- a/src/rejection.c +++ b/src/rejection.c @@ -94,9 +94,10 @@ void early_rejection(Crystal **crystals, int n) } -static double calculate_cchalf(RefList *full, Crystal *exclude) +static double calculate_cchalf(RefList *template, RefList *full, + Crystal *exclude) { - Reflection *refl; + Reflection *trefl; RefListIterator *iter; double all_sum_mean = 0.0; double all_sumsq_mean = 0.0; @@ -105,15 +106,21 @@ static double calculate_cchalf(RefList *full, Crystal *exclude) int n = 0; /* Iterate over all reflections */ - for ( refl = first_refl(full, &iter); - refl != NULL; - refl = next_refl(refl, iter) ) + for ( trefl = first_refl(template, &iter); + trefl != NULL; + trefl = next_refl(trefl, iter) ) { struct reflection_contributions *c; int j; double refl_sum = 0.0; double refl_sumsq = 0.0; double refl_mean, refl_var; + signed int h, k, l; + Reflection *refl; + + get_indices(trefl, &h, &k, &l); + refl = find_refl(full, h, k, l); + if ( refl == NULL ) continue; c = get_contributions(refl); assert(c != NULL); @@ -156,12 +163,14 @@ static void check_deltacchalf(Crystal **crystals, int n, RefList *full) double cchalf; int i; - cchalf = calculate_cchalf(full, NULL); STATUS("Overall CChalf = %f\n", cchalf); for ( i=0; i Date: Thu, 30 Aug 2018 16:50:28 +0200 Subject: Add Bessel correction --- src/rejection.c | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) (limited to 'src') diff --git a/src/rejection.c b/src/rejection.c index 46e7a2f2..cb3bede5 100644 --- a/src/rejection.c +++ b/src/rejection.c @@ -116,6 +116,7 @@ static double calculate_cchalf(RefList *template, RefList *full, double refl_sumsq = 0.0; double refl_mean, refl_var; signed int h, k, l; + double bcorr; Reflection *refl; get_indices(trefl, &h, &k, &l); @@ -142,7 +143,8 @@ static double calculate_cchalf(RefList *template, RefList *full, } refl_mean = refl_sum / c->n_contrib; - refl_var = refl_sumsq - refl_sum*refl_sum; + bcorr = c->n_contrib / (c->n_contrib - 1); + refl_var = (refl_sumsq - refl_sum*refl_sum)*bcorr; all_sum_mean += refl_mean; all_sumsq_mean += refl_mean*refl_mean; @@ -152,7 +154,7 @@ static double calculate_cchalf(RefList *template, RefList *full, } sig2E = all_sum_var / n; - sig2Y = all_sumsq_mean - all_sum_mean*all_sum_mean; + sig2Y = (all_sumsq_mean - all_sum_mean*all_sum_mean)*n/(n-1); return (sig2Y - 0.5*sig2E) / (sig2Y + 0.5*sig2E); } -- cgit v1.2.3 From 90f96d061805b604cb976d9ab5aa1cf236269cc5 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Thu, 30 Aug 2018 17:09:15 +0200 Subject: =?UTF-8?q?Actually=20bother=20to=20calculate=20overall=20CC=C2=BD?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- src/rejection.c | 1 + 1 file changed, 1 insertion(+) (limited to 'src') diff --git a/src/rejection.c b/src/rejection.c index cb3bede5..6ff7265d 100644 --- a/src/rejection.c +++ b/src/rejection.c @@ -165,6 +165,7 @@ static void check_deltacchalf(Crystal **crystals, int n, RefList *full) double cchalf; int i; + cchalf = calculate_cchalf(full, full, NULL); STATUS("Overall CChalf = %f\n", cchalf); for ( i=0; i Date: Fri, 31 Aug 2018 20:31:07 +0200 Subject: Show number of reflections used for CC calculation --- src/rejection.c | 24 ++++++++++++++++-------- 1 file changed, 16 insertions(+), 8 deletions(-) (limited to 'src') diff --git a/src/rejection.c b/src/rejection.c index 6ff7265d..47cc9347 100644 --- a/src/rejection.c +++ b/src/rejection.c @@ -95,7 +95,7 @@ void early_rejection(Crystal **crystals, int n) static double calculate_cchalf(RefList *template, RefList *full, - Crystal *exclude) + Crystal *exclude, int *pnref) { Reflection *trefl; RefListIterator *iter; @@ -156,6 +156,9 @@ static double calculate_cchalf(RefList *template, RefList *full, sig2E = all_sum_var / n; sig2Y = (all_sumsq_mean - all_sum_mean*all_sum_mean)*n/(n-1); + if ( pnref != NULL ) { + *pnref = n; + } return (sig2Y - 0.5*sig2E) / (sig2Y + 0.5*sig2E); } @@ -164,17 +167,22 @@ static void check_deltacchalf(Crystal **crystals, int n, RefList *full) { double cchalf; int i; + int nref; - cchalf = calculate_cchalf(full, full, NULL); - STATUS("Overall CChalf = %f\n", cchalf); + cchalf = calculate_cchalf(full, full, NULL, &nref); + STATUS("Overall CChalf = %f (%i reflections)\n", cchalf*100.0, nref); for ( i=0; i Date: Fri, 31 Aug 2018 20:31:30 +0200 Subject: Handle each unique reflection only once --- src/rejection.c | 10 ++++++++++ 1 file changed, 10 insertions(+) (limited to 'src') diff --git a/src/rejection.c b/src/rejection.c index 47cc9347..b4e16f17 100644 --- a/src/rejection.c +++ b/src/rejection.c @@ -104,6 +104,9 @@ static double calculate_cchalf(RefList *template, RefList *full, double all_sum_var = 0.0; double sig2E, sig2Y; int n = 0; + signed int oh = 0; + signed int ok = 0; + signed int ol = 0; /* Iterate over all reflections */ for ( trefl = first_refl(template, &iter); @@ -120,6 +123,13 @@ static double calculate_cchalf(RefList *template, RefList *full, Reflection *refl; 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; -- cgit v1.2.3 From 51ee50053c9658d75ebd6df229b054a49d040c7d Mon Sep 17 00:00:00 2001 From: Thomas White Date: Fri, 31 Aug 2018 20:32:06 +0200 Subject: CC fixes --- src/rejection.c | 18 +++++++++++------- 1 file changed, 11 insertions(+), 7 deletions(-) (limited to 'src') diff --git a/src/rejection.c b/src/rejection.c index b4e16f17..b6cb098c 100644 --- a/src/rejection.c +++ b/src/rejection.c @@ -102,7 +102,7 @@ static double calculate_cchalf(RefList *template, RefList *full, double all_sum_mean = 0.0; double all_sumsq_mean = 0.0; double all_sum_var = 0.0; - double sig2E, sig2Y; + double hsig2E, sig2Y; int n = 0; signed int oh = 0; signed int ok = 0; @@ -149,12 +149,15 @@ static double calculate_cchalf(RefList *template, RefList *full, refl_sum += Ii; refl_sumsq += Ii * Ii; + nc++; } - refl_mean = refl_sum / c->n_contrib; - bcorr = c->n_contrib / (c->n_contrib - 1); - refl_var = (refl_sumsq - refl_sum*refl_sum)*bcorr; + if ( nc < 2 ) continue; + + refl_mean = refl_sum / nc; + refl_var = (refl_sumsq/nc - refl_mean*refl_mean); + refl_var *= nc/(nc-1.0); all_sum_mean += refl_mean; all_sumsq_mean += refl_mean*refl_mean; @@ -163,13 +166,14 @@ static double calculate_cchalf(RefList *template, RefList *full, } - sig2E = all_sum_var / n; - sig2Y = (all_sumsq_mean - all_sum_mean*all_sum_mean)*n/(n-1); + hsig2E = all_sum_var / n; + sig2Y = (all_sumsq_mean/n - all_sum_mean*all_sum_mean/(n*n)); + sig2Y *= n/(n-1.0); if ( pnref != NULL ) { *pnref = n; } - return (sig2Y - 0.5*sig2E) / (sig2Y + 0.5*sig2E); + return (sig2Y - hsig2E) / (sig2Y + hsig2E); } -- cgit v1.2.3 From 9a87658c6a22f910f196038e1ea910c634719d76 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Fri, 31 Aug 2018 20:32:34 +0200 Subject: Fixup variables --- src/rejection.c | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) (limited to 'src') diff --git a/src/rejection.c b/src/rejection.c index b6cb098c..1dccd2da 100644 --- a/src/rejection.c +++ b/src/rejection.c @@ -119,7 +119,7 @@ static double calculate_cchalf(RefList *template, RefList *full, double refl_sumsq = 0.0; double refl_mean, refl_var; signed int h, k, l; - double bcorr; + int nc = 0; Reflection *refl; get_indices(trefl, &h, &k, &l); -- cgit v1.2.3 From 71933aa9868e6a6502b3ef5aff976d3c072a882d Mon Sep 17 00:00:00 2001 From: Thomas White Date: Sun, 2 Sep 2018 22:37:00 +0200 Subject: Fussiness --- src/merge.c | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) (limited to 'src') diff --git a/src/merge.c b/src/merge.c index 9e577224..3b850c30 100644 --- a/src/merge.c +++ b/src/merge.c @@ -317,7 +317,7 @@ 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); } } -- cgit v1.2.3 From ea7115f694da034319bf3a643fc335336b1d3512 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Sun, 2 Sep 2018 22:38:35 +0200 Subject: New variance formula --- src/rejection.c | 50 ++++++++++++++++++++++++++++++++++++++------------ 1 file changed, 38 insertions(+), 12 deletions(-) (limited to 'src') diff --git a/src/rejection.c b/src/rejection.c index 1dccd2da..eca55f2f 100644 --- a/src/rejection.c +++ b/src/rejection.c @@ -99,14 +99,17 @@ static double calculate_cchalf(RefList *template, RefList *full, { Reflection *trefl; RefListIterator *iter; - double all_sum_mean = 0.0; - double all_sumsq_mean = 0.0; - double all_sum_var = 0.0; double hsig2E, 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); @@ -115,8 +118,8 @@ static double calculate_cchalf(RefList *template, RefList *full, { struct reflection_contributions *c; int j; - double refl_sum = 0.0; - double refl_sumsq = 0.0; + double refl_sum; + double refl_sumsq; double refl_mean, refl_var; signed int h, k, l; int nc = 0; @@ -136,6 +139,8 @@ static double calculate_cchalf(RefList *template, RefList *full, c = get_contributions(refl); assert(c != NULL); + /* Mean of contributions */ + refl_sum = 0.0; for ( j=0; jn_contrib; j++ ) { double Ii; @@ -148,7 +153,6 @@ static double calculate_cchalf(RefList *template, RefList *full, Ii = get_intensity(c->contribs[j]); refl_sum += Ii; - refl_sumsq += Ii * Ii; nc++; } @@ -156,19 +160,41 @@ static double calculate_cchalf(RefList *template, RefList *full, if ( nc < 2 ) continue; refl_mean = refl_sum / nc; - refl_var = (refl_sumsq/nc - refl_mean*refl_mean); - refl_var *= nc/(nc-1.0); - all_sum_mean += refl_mean; - all_sumsq_mean += refl_mean*refl_mean; + /* Variance of contributions */ + refl_sumsq = 0.0; + for ( j=0; jn_contrib; j++ ) { + + double Ii; + + if ( c->contrib_crystals[j] == exclude ) { + continue; + } + + /* FIXME: Apply corrections */ + Ii = get_intensity(c->contribs[j]); + + refl_sumsq += (Ii-refl_mean)*(Ii-refl_mean); + /* (nc already summed above) */ + + } + + refl_var = refl_sumsq/(nc-1.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); } hsig2E = all_sum_var / n; - sig2Y = (all_sumsq_mean/n - all_sum_mean*all_sum_mean/(n*n)); - sig2Y *= n/(n-1.0); + sig2Y = S / (wSum - 1.0); if ( pnref != NULL ) { *pnref = n; -- cgit v1.2.3 From 4f1950b8799ed98a82f45d84a7c6ef287fbfbc5c Mon Sep 17 00:00:00 2001 From: Thomas White Date: Sun, 2 Sep 2018 22:52:45 +0200 Subject: DeltaCChalf all on one line --- src/rejection.c | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) (limited to 'src') diff --git a/src/rejection.c b/src/rejection.c index eca55f2f..fbe631ab 100644 --- a/src/rejection.c +++ b/src/rejection.c @@ -218,11 +218,11 @@ static void check_deltacchalf(Crystal **crystals, int n, RefList *full) //RefList *template = full; cchalf = calculate_cchalf(template, full, NULL, NULL); cchalfi = calculate_cchalf(template, full, crystals[i], &nref); - STATUS("Frame %i:\n", i); - STATUS(" With = %f\n", cchalf*100.0); - STATUS("Without = %f\n", cchalfi*100.0); - STATUS("Delta = %f\n", cchalf - cchalfi); - STATUS("nref = %i\n", nref); + STATUS("Frame %i:", i); + STATUS(" With = %f ", cchalf*100.0); + STATUS("Without = %f", cchalfi*100.0); + STATUS(" Delta = %f ", cchalf - cchalfi); + STATUS("(nref = %i)\n", nref); } } -- cgit v1.2.3 From 2e50a23bfab56fdfad0ffc5f6c3eef6a65346644 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Tue, 4 Sep 2018 13:44:29 +0200 Subject: Debugging stuff --- src/rejection.c | 9 ++++++++- 1 file changed, 8 insertions(+), 1 deletion(-) (limited to 'src') diff --git a/src/rejection.c b/src/rejection.c index fbe631ab..83738657 100644 --- a/src/rejection.c +++ b/src/rejection.c @@ -139,6 +139,8 @@ static double calculate_cchalf(RefList *template, RefList *full, c = get_contributions(refl); assert(c != NULL); + STATUS("reflection %4i %4i %4i\n", h, k, l); + /* Mean of contributions */ refl_sum = 0.0; for ( j=0; jn_contrib; j++ ) { @@ -152,14 +154,15 @@ static double calculate_cchalf(RefList *template, RefList *full, /* FIXME: Apply corrections */ Ii = get_intensity(c->contribs[j]); + STATUS("contrib %f (crystal %p)\n", Ii, c->contrib_crystals[j]); refl_sum += Ii; nc++; } if ( nc < 2 ) continue; - refl_mean = refl_sum / nc; + STATUS(" refl_mean = %f\n", refl_mean); /* Variance of contributions */ refl_sumsq = 0.0; @@ -180,6 +183,7 @@ static double calculate_cchalf(RefList *template, RefList *full, } refl_var = refl_sumsq/(nc-1.0); + STATUS(" refl_var = %f\n", refl_var); all_sum_var += refl_var; n++; total_contribs += nc; @@ -194,7 +198,10 @@ static double calculate_cchalf(RefList *template, RefList *full, } hsig2E = all_sum_var / n; + STATUS("mean of reflection variances = %f = 0.5*sig2E\n", hsig2E); sig2Y = S / (wSum - 1.0); + STATUS("variance of reflection means = %f = sig2Y\n", sig2Y); + STATUS("%lli total contribs\n", hsig2E, sig2Y, total_contribs); if ( pnref != NULL ) { *pnref = n; -- cgit v1.2.3 From 264abc6b6b72640d78d6c213fc3f25267a7f6c1c Mon Sep 17 00:00:00 2001 From: Thomas White Date: Wed, 5 Sep 2018 10:11:34 +0200 Subject: Fix calculation and remove debugging stuff --- src/rejection.c | 9 ++------- 1 file changed, 2 insertions(+), 7 deletions(-) (limited to 'src') diff --git a/src/rejection.c b/src/rejection.c index 83738657..5ead8703 100644 --- a/src/rejection.c +++ b/src/rejection.c @@ -139,7 +139,6 @@ static double calculate_cchalf(RefList *template, RefList *full, c = get_contributions(refl); assert(c != NULL); - STATUS("reflection %4i %4i %4i\n", h, k, l); /* Mean of contributions */ refl_sum = 0.0; @@ -154,7 +153,6 @@ static double calculate_cchalf(RefList *template, RefList *full, /* FIXME: Apply corrections */ Ii = get_intensity(c->contribs[j]); - STATUS("contrib %f (crystal %p)\n", Ii, c->contrib_crystals[j]); refl_sum += Ii; nc++; @@ -162,7 +160,6 @@ static double calculate_cchalf(RefList *template, RefList *full, if ( nc < 2 ) continue; refl_mean = refl_sum / nc; - STATUS(" refl_mean = %f\n", refl_mean); /* Variance of contributions */ refl_sumsq = 0.0; @@ -183,7 +180,8 @@ static double calculate_cchalf(RefList *template, RefList *full, } refl_var = refl_sumsq/(nc-1.0); - STATUS(" refl_var = %f\n", refl_var); + refl_var /= (nc/2.0); + all_sum_var += refl_var; n++; total_contribs += nc; @@ -198,10 +196,7 @@ static double calculate_cchalf(RefList *template, RefList *full, } hsig2E = all_sum_var / n; - STATUS("mean of reflection variances = %f = 0.5*sig2E\n", hsig2E); sig2Y = S / (wSum - 1.0); - STATUS("variance of reflection means = %f = sig2Y\n", sig2Y); - STATUS("%lli total contribs\n", hsig2E, sig2Y, total_contribs); if ( pnref != NULL ) { *pnref = n; -- cgit v1.2.3 From 6bd3df7ee39b8a8ba46e9709d74d96051a27ca6f Mon Sep 17 00:00:00 2001 From: Thomas White Date: Thu, 6 Sep 2018 14:53:47 +0200 Subject: Working and validated deltaCChalf calculation --- src/rejection.c | 8 +++----- 1 file changed, 3 insertions(+), 5 deletions(-) (limited to 'src') diff --git a/src/rejection.c b/src/rejection.c index 5ead8703..5a82a940 100644 --- a/src/rejection.c +++ b/src/rejection.c @@ -99,7 +99,7 @@ static double calculate_cchalf(RefList *template, RefList *full, { Reflection *trefl; RefListIterator *iter; - double hsig2E, sig2Y; + double sig2E, sig2Y; int n = 0; signed int oh = 0; signed int ok = 0; @@ -139,7 +139,6 @@ static double calculate_cchalf(RefList *template, RefList *full, c = get_contributions(refl); assert(c != NULL); - /* Mean of contributions */ refl_sum = 0.0; for ( j=0; jn_contrib; j++ ) { @@ -195,13 +194,13 @@ static double calculate_cchalf(RefList *template, RefList *full, } - hsig2E = all_sum_var / n; + sig2E = all_sum_var / n; sig2Y = S / (wSum - 1.0); if ( pnref != NULL ) { *pnref = n; } - return (sig2Y - hsig2E) / (sig2Y + hsig2E); + return (sig2Y - 0.5*sig2E) / (sig2Y + 0.5*sig2E); } @@ -217,7 +216,6 @@ static void check_deltacchalf(Crystal **crystals, int n, RefList *full) for ( i=0; i Date: Thu, 6 Sep 2018 16:23:39 +0200 Subject: Don't forget to multiple the deltaCChalf values by 100 --- src/rejection.c | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) (limited to 'src') diff --git a/src/rejection.c b/src/rejection.c index 5a82a940..01ed2cf0 100644 --- a/src/rejection.c +++ b/src/rejection.c @@ -221,7 +221,7 @@ static void check_deltacchalf(Crystal **crystals, int n, RefList *full) STATUS("Frame %i:", i); STATUS(" With = %f ", cchalf*100.0); STATUS("Without = %f", cchalfi*100.0); - STATUS(" Delta = %f ", cchalf - cchalfi); + STATUS(" Delta = %f ", (cchalf - cchalfi)*100.0); STATUS("(nref = %i)\n", nref); } } -- cgit v1.2.3 From a61dc5721a2896267dcd783389908bf9f54ff21a Mon Sep 17 00:00:00 2001 From: Thomas White Date: Fri, 7 Sep 2018 15:30:07 +0200 Subject: s/PRFLAG_CC/PRFLAG_DELTACCHALF/ --- src/post-refinement.c | 2 +- src/post-refinement.h | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) (limited to 'src') diff --git a/src/post-refinement.c b/src/post-refinement.c index ebaf66ab..0a799b9e 100644 --- a/src/post-refinement.c +++ b/src/post-refinement.c @@ -69,7 +69,7 @@ const char *str_prflag(enum prflag flag) case PRFLAG_EARLY : return "early rejection"; - case PRFLAG_CC : + case PRFLAG_DELTACCHALF : return "negative delta CC½"; case PRFLAG_BIGB : 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, }; -- cgit v1.2.3 From 6ea010ba017579d0197f45c940f03d04800e24e7 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Fri, 7 Sep 2018 15:50:41 +0200 Subject: partialator: Add option to disable deltaCChalf --- src/partialator.c | 7 +++++-- src/rejection.c | 11 +++++++---- src/rejection.h | 6 +++--- 3 files changed, 15 insertions(+), 9 deletions(-) (limited to 'src') diff --git a/src/partialator.c b/src/partialator.c index 113535e0..fa07eab7 100644 --- a/src/partialator.c +++ b/src/partialator.c @@ -318,6 +318,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= Specify partiality model.\n" " --min-measurements= Minimum number of measurements to require.\n" " --no-polarisation Disable polarisation correction.\n" @@ -912,6 +913,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 +952,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 +1415,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 ) { @@ -1443,7 +1446,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 ) { diff --git a/src/rejection.c b/src/rejection.c index 01ed2cf0..371216a9 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 + * 2010-2018 Thomas White * * This file is part of CrystFEL. * @@ -254,13 +254,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 delta CC½ */ - if ( full != NULL ) check_deltacchalf(crystals, n, full); + if ( !no_deltacchalf && (full != NULL) ) { + check_deltacchalf(crystals, n, full); + } for ( i=0; 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 + * 2010-2018 Thomas White * * 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 */ -- cgit v1.2.3 From 6dc774e11fa507ba0b667b47e112acab4c7e4f15 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Fri, 7 Sep 2018 15:51:23 +0200 Subject: Add the actual deltaCChalf rejection Initially set at mean deltaCChalf minus 2 sigma --- src/rejection.c | 21 +++++++++++++++++++++ 1 file changed, 21 insertions(+) (limited to 'src') diff --git a/src/rejection.c b/src/rejection.c index 371216a9..7b78803e 100644 --- a/src/rejection.c +++ b/src/rejection.c @@ -209,10 +209,18 @@ 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 Date: Mon, 10 Sep 2018 14:01:06 +0200 Subject: Apply scaling corrections when calculating deltaCChalf, and add progress bar --- src/rejection.c | 35 ++++++++++++++++++++++++++--------- 1 file changed, 26 insertions(+), 9 deletions(-) (limited to 'src') diff --git a/src/rejection.c b/src/rejection.c index 7b78803e..9f298c0b 100644 --- a/src/rejection.c +++ b/src/rejection.c @@ -124,6 +124,7 @@ static double calculate_cchalf(RefList *template, RefList *full, signed int h, k, l; int nc = 0; Reflection *refl; + double res; get_indices(trefl, &h, &k, &l); @@ -139,19 +140,29 @@ static double calculate_cchalf(RefList *template, RefList *full, c = get_contributions(refl); assert(c != NULL); + /* 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); + /* Mean of contributions */ refl_sum = 0.0; for ( j=0; jn_contrib; j++ ) { - double Ii; + double Ii, G, B; if ( c->contrib_crystals[j] == exclude ) { continue; } - /* FIXME: Apply corrections */ + 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++; @@ -164,15 +175,20 @@ static double calculate_cchalf(RefList *template, RefList *full, refl_sumsq = 0.0; for ( j=0; jn_contrib; j++ ) { - double Ii; + double Ii, G, B; if ( c->contrib_crystals[j] == exclude ) { continue; } - /* FIXME: Apply corrections */ + 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) */ @@ -226,12 +242,13 @@ static void check_deltacchalf(Crystal **crystals, int n, RefList *full) 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); + //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); -- cgit v1.2.3 From e8fa57ba93687fbd6841cb3a3d09c87037b5c713 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Mon, 10 Sep 2018 16:48:31 +0200 Subject: Make sure reflection contribution list gets freed --- src/merge.c | 15 +++++++++++++-- src/partialator.c | 6 ++++++ src/scaling.c | 2 ++ 3 files changed, 21 insertions(+), 2 deletions(-) (limited to 'src') diff --git a/src/merge.c b/src/merge.c index 3b850c30..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 + * 2010-2018 Thomas White * * 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" @@ -319,6 +320,16 @@ 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); + + } 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 fa07eab7..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, @@ -1436,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, @@ -1484,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); @@ -1530,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/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++; -- cgit v1.2.3 From f22a30956317ccdc48d29eaf8d95e32e01f0be34 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Mon, 10 Sep 2018 17:17:21 +0200 Subject: Add % for overall CChalf --- src/rejection.c | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) (limited to 'src') diff --git a/src/rejection.c b/src/rejection.c index 9f298c0b..230093e6 100644 --- a/src/rejection.c +++ b/src/rejection.c @@ -229,7 +229,7 @@ static void check_deltacchalf(Crystal **crystals, int n, RefList *full) double mean, sd; cchalf = calculate_cchalf(full, full, NULL, &nref); - STATUS("Overall CChalf = %f (%i reflections)\n", cchalf*100.0, nref); + STATUS("Overall CChalf = %f %% (%i reflections)\n", cchalf*100.0, nref); vals = malloc(n*sizeof(double)); if ( vals == NULL ) { -- cgit v1.2.3