diff options
-rw-r--r-- | libcrystfel/src/fom.c | 612 | ||||
-rw-r--r-- | libcrystfel/src/fom.h | 89 | ||||
-rw-r--r-- | src/check_hkl.c | 384 | ||||
-rw-r--r-- | src/compare_hkl.c | 118 |
4 files changed, 715 insertions, 488 deletions
diff --git a/libcrystfel/src/fom.c b/libcrystfel/src/fom.c index 19f35a86..c2de704a 100644 --- a/libcrystfel/src/fom.c +++ b/libcrystfel/src/fom.c @@ -43,6 +43,33 @@ #include "reflist.h" #include "reflist-utils.h" +struct fom_context +{ + enum fom_type fom; + int nshells; + int *cts; + + /* For R-factors */ + double *num; + double *den; + + /* For "double" R-factors */ + double *num2; + double *den2; + + /* For CCs */ + double **vec1; + double **vec2; + int *n; + int nmax; + + /* For "counting" things e.g. d1sig or d2sig */ + int *n_within; + + long int *n_meas; + long int *possible; +}; + enum fom_type fom_type_from_string(const char *s) { if ( strcasecmp(s, "r1i") == 0 ) return FOM_R1I; @@ -78,6 +105,12 @@ static struct fom_context *init_fom(enum fom_type fom, int nmax, int nshells) fctx->cts[i] = 0; } + fctx->num = NULL; + fctx->den = NULL; + fctx->num2 = NULL; + fctx->den2 = NULL; + fctx->possible = NULL; + switch ( fctx->fom ) { case FOM_RANORSPLIT : @@ -95,6 +128,10 @@ static struct fom_context *init_fom(enum fom_type fom, int nmax, int nshells) case FOM_R2 : case FOM_RSPLIT : case FOM_RANO : + case FOM_MEAN_INTENSITY : + case FOM_STDDEV_INTENSITY : + case FOM_SNR : + case FOM_REDUNDANCY : fctx->num = malloc(nshells*sizeof(double)); fctx->den = malloc(nshells*sizeof(double)); if ( (fctx->num == NULL) || (fctx->den == NULL) ) return NULL; @@ -104,6 +141,15 @@ static struct fom_context *init_fom(enum fom_type fom, int nmax, int nshells) } break; + case FOM_COMPLETENESS : + /* Uses 'cts' and 'possible' only */ + break; + + case FOM_NUM_MEASUREMENTS : + fctx->n_meas = calloc(nshells, sizeof(long int)); + if ( fctx->n_meas == NULL ) return NULL; + break; + case FOM_CC : case FOM_CCSTAR : case FOM_CCANO : @@ -140,37 +186,45 @@ static struct fom_context *init_fom(enum fom_type fom, int nmax, int nshells) } -static void add_to_fom(struct fom_context *fctx, double i1, double i2, - double i1bij, double i2bij, double sig1, double sig2, - int bin) +static int add_to_fom(struct fom_context *fctx, + Reflection *refl1, + Reflection *refl2, + Reflection *refl1bij, + Reflection *refl2bij, + int bin) { - double f1, f2; + double i1, i2, i1bij, i2bij, sig1, sig2; double im, imbij; + int bad = 0; fctx->cts[bin]++; - /* Negative intensities have already been weeded out. */ - f1 = sqrt(i1); - f2 = sqrt(i2); - switch ( fctx->fom ) { case FOM_R1I : + i1 = get_intensity(refl1); + i2 = get_intensity(refl2); fctx->num[bin] += fabs(i1 - i2); fctx->den[bin] += i1; break; case FOM_R1F : - fctx->num[bin] += fabs(f1 - f2); - fctx->den[bin] += f1; + i1 = get_intensity(refl1); + i2 = get_intensity(refl2); + fctx->num[bin] += fabs(sqrt(i1) - sqrt(i2)); + fctx->den[bin] += sqrt(i1); break; case FOM_R2 : + i1 = get_intensity(refl1); + i2 = get_intensity(refl2); fctx->num[bin] += pow(i1 - i2, 2.0); fctx->den[bin] += pow(i1, 2.0); break; case FOM_RSPLIT : + i1 = get_intensity(refl1); + i2 = get_intensity(refl2); fctx->num[bin] += fabs(i1 - i2); fctx->den[bin] += i1 + i2; break; @@ -178,6 +232,8 @@ static void add_to_fom(struct fom_context *fctx, double i1, double i2, case FOM_CC : case FOM_CCSTAR : assert(fctx->n[bin] < fctx->nmax); + i1 = get_intensity(refl1); + i2 = get_intensity(refl2); fctx->vec1[bin][fctx->n[bin]] = i1; fctx->vec2[bin][fctx->n[bin]] = i2; fctx->n[bin]++; @@ -186,17 +242,27 @@ static void add_to_fom(struct fom_context *fctx, double i1, double i2, case FOM_CCANO : case FOM_CRDANO : assert(fctx->n[bin] < fctx->nmax); + i1 = get_intensity(refl1); + i2 = get_intensity(refl2); + i1bij = get_intensity(refl1bij); + i2bij = get_intensity(refl2bij); fctx->vec1[bin][fctx->n[bin]] = i1 - i1bij; fctx->vec2[bin][fctx->n[bin]] = i2 - i2bij; fctx->n[bin]++; break; case FOM_RANORSPLIT : + i1 = get_intensity(refl1); + i2 = get_intensity(refl2); fctx->num2[bin] += fabs(i1 - i2); fctx->den2[bin] += i1 + i2; /* Intentional fall-through (no break) */ case FOM_RANO : + i1 = get_intensity(refl1); + i2 = get_intensity(refl2); + i1bij = get_intensity(refl1bij); + i2bij = get_intensity(refl2bij); im = (i1 + i2)/2.0; imbij = (i1bij + i2bij)/2.0; fctx->num[bin] += fabs(im - imbij); @@ -204,22 +270,69 @@ static void add_to_fom(struct fom_context *fctx, double i1, double i2, break; case FOM_D1SIG : + i1 = get_intensity(refl1); + i2 = get_intensity(refl2); + sig1 = get_esd_intensity(refl1); + sig2 = get_esd_intensity(refl2); if ( fabs(i1-i2) < sqrt(sig1*sig1 + sig2*sig2) ) { fctx->n_within[bin]++; } break; case FOM_D2SIG : + i1 = get_intensity(refl1); + i2 = get_intensity(refl2); + sig1 = get_esd_intensity(refl1); + sig2 = get_esd_intensity(refl2); if ( fabs(i1-i2) < 2.0*sqrt(sig1*sig1 + sig2*sig2) ) { fctx->n_within[bin]++; } break; + case FOM_NUM_MEASUREMENTS : + fctx->n_meas[bin] += get_redundancy(refl1); + break; + + case FOM_REDUNDANCY : + fctx->num[bin] += get_redundancy(refl1); + fctx->den[bin] += 1.0; + break; + + case FOM_SNR : + i1 = get_intensity(refl1); + sig1 = get_esd_intensity(refl1); + if ( isfinite(i1/sig1) ) { + fctx->num[bin] += i1/sig1; + fctx->den[bin] += 1.0; + } else { + bad = 1; + } + break; + + case FOM_MEAN_INTENSITY : + i1 = get_intensity(refl1); + fctx->num[bin] += i1; + fctx->den[bin] += 1.0; + break; + + /* FIXME: Delete this FoM? */ + case FOM_STDDEV_INTENSITY : + fctx->num[bin] += 1.0; + fctx->den[bin] += 1.0; + break; + + case FOM_COMPLETENESS : + /* fctx->cts already incremented, as needed. + * Will calculate possible reflections later */ + break; + } + + return bad; } -double fom_overall(struct fom_context *fctx) +double fom_overall_value(struct fom_context *fctx) { double overall_num = INFINITY; double overall_den = 0.0; @@ -234,6 +347,9 @@ double fom_overall(struct fom_context *fctx) double variance_signal; double variance_error; double cc = INFINITY; + long int total_meas = 0; + long int overall_cts = 0; + long int overall_possible = 0; switch ( fctx->fom ) { @@ -242,6 +358,10 @@ double fom_overall(struct fom_context *fctx) case FOM_R2 : case FOM_RSPLIT : case FOM_RANO : + case FOM_REDUNDANCY : + case FOM_SNR : + case FOM_MEAN_INTENSITY : + case FOM_STDDEV_INTENSITY : overall_num = 0.0; overall_den = 0.0; for ( i=0; i<fctx->nshells; i++ ) { @@ -321,14 +441,38 @@ double fom_overall(struct fom_context *fctx) } break; + case FOM_NUM_MEASUREMENTS : + total_meas = 0; + for ( i=0; i<fctx->nshells; i++ ) { + total_meas += fctx->n_meas[i]; + } + break; + + case FOM_COMPLETENESS : + for ( i=0; i<fctx->nshells; i++ ) { + overall_cts += fctx->cts[i]; + overall_possible += fctx->possible[i]; + } + break; + } switch ( fctx->fom ) { case FOM_R1I : case FOM_R1F : + case FOM_REDUNDANCY : + case FOM_SNR : + case FOM_MEAN_INTENSITY : + case FOM_STDDEV_INTENSITY : return overall_num/overall_den; + case FOM_COMPLETENESS : + return (double)overall_cts / overall_possible; + + case FOM_NUM_MEASUREMENTS : + return total_meas; + case FOM_R2 : return sqrt(overall_num/overall_den); @@ -361,7 +505,7 @@ double fom_overall(struct fom_context *fctx) } -double fom_shell(struct fom_context *fctx, int i) +double fom_shell_value(struct fom_context *fctx, int i) { double cc; int j; @@ -374,6 +518,10 @@ double fom_shell(struct fom_context *fctx, int i) case FOM_R1I : case FOM_R1F : + case FOM_REDUNDANCY : + case FOM_SNR : + case FOM_MEAN_INTENSITY : + case FOM_STDDEV_INTENSITY : return fctx->num[i]/fctx->den[i]; case FOM_R2 : @@ -420,6 +568,12 @@ double fom_shell(struct fom_context *fctx, int i) case FOM_D2SIG : return (double)fctx->n_within[i] / fctx->cts[i]; + case FOM_NUM_MEASUREMENTS : + return fctx->n_meas[i]; + + case FOM_COMPLETENESS : + return (double)fctx->cts[i] / fctx->possible[i]; + } ERROR("This point is never reached.\n"); @@ -469,7 +623,7 @@ struct fom_shells *fom_make_resolution_shells(double rmin, double rmax, } -double fom_shell_label(struct fom_shells *s, int i) +double fom_shell_centre(struct fom_shells *s, int i) { return s->rmins[i] + (s->rmaxs[i] - s->rmins[i])/2.0; } @@ -605,15 +759,143 @@ static int wilson_scale(RefList *list1, RefList *list2, UnitCell *cell) } +static int calculate_possible(struct fom_context *fctx, + struct fom_shells *shells, + UnitCell *cell, + const SymOpList *sym) +{ + RefList *counted; + int hmax, kmax, lmax; + double ax, ay, az; + double bx, by, bz; + double cx, cy, cz; + signed int h, k, l; + + fctx->possible = calloc(fctx->nshells, sizeof(long int)); + if ( fctx->possible == NULL ) return 1; + + counted = reflist_new(); + if ( counted == NULL ) { + free(fctx->possible); + return 1; + } + + cell_get_cartesian(cell, &ax, &ay, &az, + &bx, &by, &bz, + &cx, &cy, &cz); + hmax = shells->rmaxs[fctx->nshells-1] * modulus(ax, ay, az); + kmax = shells->rmaxs[fctx->nshells-1] * modulus(bx, by, bz); + lmax = shells->rmaxs[fctx->nshells-1] * modulus(cx, cy, cz); + for ( h=-hmax; h<=hmax; h++ ) { + for ( k=-kmax; k<=kmax; k++ ) { + for ( l=-lmax; l<=lmax; l++ ) { + + double d; + signed int hs, ks, ls; + int bin; + int i; + + get_asymm(sym, h, k, l, &hs, &ks, &ls); + d = 2.0 * resolution(cell, hs, ks, ls); + + if ( forbidden_reflection(cell, h, k, l) ) continue; + + bin = -1; + for ( i=0; i<fctx->nshells; i++ ) { + if ( (d>shells->rmins[i]) && (d<=shells->rmaxs[i]) ) { + bin = i; + break; + } + } + if ( bin == -1 ) continue; + + if ( find_refl(counted, hs, ks, ls) != NULL ) continue; + add_refl(counted, hs, ks, ls); + + fctx->possible[bin]++; + + } + } + } + reflist_free(counted); + + return 0; +} + + +static int is_anomalous(enum fom_type fom) +{ + switch ( fom ) { + + case FOM_CCANO: + case FOM_RANO: + case FOM_CRDANO: + case FOM_RANORSPLIT: + return 1; + + case FOM_R1I: + case FOM_R1F: + case FOM_R2: + case FOM_RSPLIT: + case FOM_CC: + case FOM_CCSTAR: + case FOM_D1SIG: + case FOM_D2SIG: + case FOM_NUM_MEASUREMENTS: + case FOM_REDUNDANCY: + case FOM_SNR: + case FOM_MEAN_INTENSITY: + case FOM_STDDEV_INTENSITY: + case FOM_COMPLETENESS: + return 0; + } + + ERROR("This point never reached\n"); + abort(); +} + + +static int is_single_list(enum fom_type fom) +{ + switch ( fom ) { + + case FOM_CCANO: + case FOM_RANO: + case FOM_CRDANO: + case FOM_RANORSPLIT: + case FOM_R1I: + case FOM_R1F: + case FOM_R2: + case FOM_RSPLIT: + case FOM_CC: + case FOM_CCSTAR: + case FOM_D1SIG: + case FOM_D2SIG: + return 0; + + case FOM_NUM_MEASUREMENTS: + case FOM_REDUNDANCY: + case FOM_SNR: + case FOM_MEAN_INTENSITY: + case FOM_STDDEV_INTENSITY: + case FOM_COMPLETENESS: + return 1; + } + + ERROR("This point never reached\n"); + abort(); +} + struct fom_context *fom_calculate(RefList *list1, RefList *list2, UnitCell *cell, struct fom_shells *shells, enum fom_type fom, - int noscale, SymOpList *sym) + int noscale, const SymOpList *sym) { Reflection *refl1; RefListIterator *iter; struct fom_context *fctx; - int n_out; + long int n_out = 0; + long int n_rej = 0; fctx = init_fom(fom, num_reflections(list1), shells->nshells); @@ -622,40 +904,44 @@ struct fom_context *fom_calculate(RefList *list1, RefList *list2, UnitCell *cell return NULL; } - if ( !noscale && wilson_scale(list1, list2, cell) ) { - ERROR("Error with scaling.\n"); - return NULL; - } + if ( !is_single_list(fom) ) { + if ( !noscale && wilson_scale(list1, list2, cell) ) { + ERROR("Error with scaling.\n"); + return NULL; + } - for ( refl1 = first_refl(list1, &iter); - refl1 != NULL; - refl1 = next_refl(refl1, iter) ) - { - Reflection *refl2; - signed int h, k, l; - set_flag(refl1, 0); - get_indices(refl1, &h, &k, &l); - refl2 = find_refl(list2, h, k, l); - assert(refl2 != NULL); - set_flag(refl2, 0); + for ( refl1 = first_refl(list1, &iter); + refl1 != NULL; + refl1 = next_refl(refl1, iter) ) + { + Reflection *refl2; + signed int h, k, l; + set_flag(refl1, 0); + get_indices(refl1, &h, &k, &l); + refl2 = find_refl(list2, h, k, l); + assert(refl2 != NULL); + set_flag(refl2, 0); + } } - n_out = 0; for ( refl1 = first_refl(list1, &iter); refl1 != NULL; refl1 = next_refl(refl1, iter) ) { signed int h, k, l; int bin; - double i1, i2; - double i1bij, i2bij; - double sig1, sig2; Reflection *refl2; + Reflection *refl1_bij = NULL; + Reflection *refl2_bij = NULL; get_indices(refl1, &h, &k, &l); - refl2 = find_refl(list2, h, k, l); - if ( refl2 == NULL ) continue; + if ( is_single_list(fom) ) { + refl2 = NULL; + } else { + refl2 = find_refl(list2, h, k, l); + if ( refl2 == NULL ) continue; + } bin = get_bin(shells, refl1, cell); if ( bin == -1 ) { @@ -663,17 +949,8 @@ struct fom_context *fom_calculate(RefList *list1, RefList *list2, UnitCell *cell continue; } - i1 = get_intensity(refl1); - i2 = get_intensity(refl2); - sig1 = get_esd_intensity(refl1); - sig2 = get_esd_intensity(refl2); + if ( is_anomalous(fom) ) { - if ( (fom == FOM_CCANO) || (fom == FOM_CRDANO) - || (fom == FOM_RANO) || (fom == FOM_RANORSPLIT) ) - { - - Reflection *refl1_bij = NULL; - Reflection *refl2_bij = NULL; signed int hb, kb, lb; if ( find_equiv_in_list(list1, -h, -k, -l, sym, @@ -700,46 +977,59 @@ struct fom_context *fom_calculate(RefList *list1, RefList *list2, UnitCell *cell assert(refl1_bij != NULL); assert(refl2_bij != NULL); - i1bij = get_intensity(refl1_bij); - i2bij = get_intensity(refl2_bij); - - } else { - - /* Make it obvious if these get used by mistake */ - i1bij = +INFINITY; - i2bij = +INFINITY; - } - add_to_fom(fctx, i1, i2, i1bij, i2bij, sig1, sig2, bin); + n_rej += add_to_fom(fctx, refl1, refl2, refl1_bij, refl2_bij, bin); + } - if ( n_out) { + if ( n_out ) { ERROR("WARNING: %i reflection pairs outside range.\n", n_out); } + if ( n_rej ) { + if ( fom == FOM_SNR ) { + ERROR("WARNING: %li reflections had infinite or " + "invalid values of I/sigma(I).\n", n_rej); + } else { + ERROR("WARNING: %li reflections rejected by add_to_fom\n", + n_rej); + } + } + + if ( fom == FOM_COMPLETENESS ) { + calculate_possible(fctx, shells, cell, sym); + } return fctx; } -int fom_select_reflections(RefList *list1, RefList *list2, - RefList *list1_acc, RefList *list2_acc, - UnitCell *cell, SymOpList *sym, - int anom, double rmin_fix, double rmax_fix, - double sigma_cutoff, int ignore_negs, - int zero_negs, int mul_cutoff) +struct fom_rejections fom_select_reflection_pairs(RefList *list1, RefList *list2, + RefList **plist1_acc, + RefList **plist2_acc, + UnitCell *cell, SymOpList *sym, + int anom, double rmin_fix, double rmax_fix, + double sigma_cutoff, int ignore_negs, + int zero_negs, int mul_cutoff) { Reflection *refl1; RefListIterator *iter; - int ncom, nrej, nmul, nneg, nres, nbij, ncen; - - /* Select reflections to be used */ - ncom = 0; - nrej = 0; - nmul = 0; - nneg = 0; - nres = 0; - nbij = 0; - ncen = 0; + struct fom_rejections rej; + RefList *list1_acc; + RefList *list2_acc; + + rej.common = 0; + rej.low_snr = 0; + rej.negative_deleted = 0; + rej.negative_zeroed = 0; + rej.few_measurements = 0; + rej.outside_resolution_range = 0; + rej.no_bijvoet = 0; + rej.centric = 0; + rej.nan_inf_value = 0; + + list1_acc = reflist_new(); + list2_acc = reflist_new(); + for ( refl1 = first_refl(list1, &iter); refl1 != NULL; refl1 = next_refl(refl1, iter) ) @@ -766,20 +1056,27 @@ int fom_select_reflections(RefList *list1, RefList *list2, mul1 = get_redundancy(refl1); mul2 = get_redundancy(refl2); + if ( !isfinite(val1) || !isfinite(val2) + || !isfinite(esd1) || !isfinite(esd2) ) + { + rej.nan_inf_value++; + continue; + } + if ( (val1 < sigma_cutoff * esd1) || (val2 < sigma_cutoff * esd2) ) { - nrej++; + rej.low_snr++; continue; } if ( ignore_negs && ((val1 < 0.0) || (val2 < 0.0)) ) { - nneg++; + rej.negative_deleted++; continue; } if ( (mul1 < mul_cutoff) || (mul2 < mul_cutoff) ) { - nmul++; + rej.few_measurements++; continue; } @@ -793,13 +1090,14 @@ int fom_select_reflections(RefList *list1, RefList *list2, val2 = 0.0; d = 1; } - if ( d ) nneg++; + if ( d ) rej.negative_zeroed++; + continue; } if ( rmin_fix > 0.0 ) { double res = 2.0*resolution(cell, h, k, l); if ( res < rmin_fix ) { - nres++; + rej.outside_resolution_range++; continue; } } @@ -807,7 +1105,7 @@ int fom_select_reflections(RefList *list1, RefList *list2, if ( rmax_fix > 0.0 ) { double res = 2.0*resolution(cell, h, k, l); if ( res > rmax_fix ) { - nres++; + rej.outside_resolution_range++; continue; } } @@ -820,7 +1118,7 @@ int fom_select_reflections(RefList *list1, RefList *list2, copy_data(refl2_acc, refl2); set_intensity(refl2_acc, val2); - ncom++; + rej.common++; } @@ -833,7 +1131,7 @@ int fom_select_reflections(RefList *list1, RefList *list2, list1_acc = reflist_new(); list2_acc = reflist_new(); - ncom = 0; + rej.common = 0; for ( refl1 = first_refl(list1, &iter); refl1 != NULL; @@ -857,7 +1155,7 @@ int fom_select_reflections(RefList *list1, RefList *list2, val2 = get_intensity(refl2); if ( is_centric(h, k, l, sym) ) { - ncen++; + rej.centric++; continue; } @@ -874,7 +1172,7 @@ int fom_select_reflections(RefList *list1, RefList *list2, } if ( (refl1_bij == NULL) || (refl2_bij == NULL) ) { - nbij++; + rej.no_bijvoet++; continue; } @@ -886,44 +1184,138 @@ int fom_select_reflections(RefList *list1, RefList *list2, copy_data(refl2_acc, refl2); set_intensity(refl2_acc, val2); - ncom++; + rej.common++; } } - if ( nrej > 0 ) { - STATUS("Discarded %i reflection pairs because either or both" - " versions had I/sigma(I) < %f.\n", nrej, sigma_cutoff); - } + *plist1_acc = list1_acc; + *plist2_acc = list2_acc; + return rej; +} - if ( ignore_negs && (nneg > 0) ) { - STATUS("Discarded %i reflection pairs because either or both" - " versions had negative intensities.\n", nneg); - } - if ( zero_negs && (nneg > 0) ) { - STATUS("For %i reflection pairs, either or both versions had" - " negative intensities which were set to zero.\n", nneg); - } +struct fom_rejections fom_select_reflections(RefList *raw_list, + RefList **plist_acc, + UnitCell *cell, SymOpList *sym, + double rmin_fix, double rmax_fix, + double sigma_cutoff, int ignore_negs, + int zero_negs, int mul_cutoff) +{ + RefList *list; + Reflection *refl; + RefListIterator *iter; + struct fom_rejections rej; - if ( nmul > 0 ) { - STATUS("%i reflection pairs rejected because either or both" - " versions had too few measurements.\n", nmul); - } + *plist_acc = NULL; + + rej.common = 0; + rej.low_snr = 0; + rej.negative_deleted = 0; + rej.negative_zeroed = 0; + rej.few_measurements = 0; + rej.outside_resolution_range = 0; + rej.no_bijvoet = 0; + rej.centric = 0; + rej.nan_inf_value = 0; + + list = reflist_new(); + if ( list == NULL ) return rej; + + for ( refl = first_refl(raw_list, &iter); + refl != NULL; + refl = next_refl(refl, iter) ) { + + signed int h, k, l; + double val, sig; + int ig = 0; + Reflection *new; + + get_indices(refl, &h, &k, &l); + + val = get_intensity(refl); + sig = get_esd_intensity(refl); + + if ( !isfinite(val) || !isfinite(sig) ) { + rej.nan_inf_value++; + continue; + } + + if ( val < sigma_cutoff * sig ) { + rej.low_snr++; + ig = 1; + } + + if ( ignore_negs && (val < 0.0) ) { + rej.negative_deleted++; + ig = 1; + } + + if ( zero_negs && (val < 0.0) ) { + set_intensity(refl, 0.0); + rej.negative_zeroed++; + } + + if ( rmin_fix > 0.0 ) { + double res = 2.0*resolution(cell, h, k, l); + if ( res < rmin_fix ) { + rej.outside_resolution_range++; + continue; + } + } - if ( nres > 0 ) { - STATUS("%i reflection pairs rejected because either or both" - " versions were outside the resolution range.\n", nres); + if ( rmax_fix > 0.0 ) { + double res = 2.0*resolution(cell, h, k, l); + if ( res > rmax_fix ) { + rej.outside_resolution_range++; + continue; + } + } + + if ( ig ) continue; + + new = add_refl(list, h, k, l); + copy_data(new, refl); } - if ( nbij > 0 ) { - STATUS("%i reflection pairs rejected because either or both" - " versions did not have Bijvoet partners.\n", nres); + *plist_acc = list; + return rej; +} + + +int fom_overall_num_reflections(struct fom_context *fctx) +{ + int i; + long int n = 0; + + for ( i=0; i<fctx->nshells; i++ ) { + n += fctx->cts[i]; } + return n; +} + - if ( ncen > 0 ) { - STATUS("%i reflection pairs rejected because they were" - " centric.\n", ncen); +int fom_shell_num_reflections(struct fom_context *fctx, int i) +{ + return fctx->cts[i]; +} + + +int fom_overall_num_possible(struct fom_context *fctx) +{ + int i; + long int n = 0; + + assert(fctx->fom == FOM_COMPLETENESS); + + for ( i=0; i<fctx->nshells; i++ ) { + n += fctx->possible[i]; } + return n; +} + - return ncom; +int fom_shell_num_possible(struct fom_context *fctx, int i) +{ + assert(fctx->fom == FOM_COMPLETENESS); + return fctx->possible[i]; } diff --git a/libcrystfel/src/fom.h b/libcrystfel/src/fom.h index d3373044..c9c8dbfe 100644 --- a/libcrystfel/src/fom.h +++ b/libcrystfel/src/fom.h @@ -38,6 +38,19 @@ #include <reflist.h> #include <symmetry.h> +struct fom_rejections +{ + int common; + int low_snr; + int negative_deleted; + int negative_zeroed; + int few_measurements; + int outside_resolution_range; + int no_bijvoet; + int centric; + int nan_inf_value; +}; + enum fom_type { FOM_R1I, @@ -51,7 +64,13 @@ enum fom_type FOM_RANO, FOM_RANORSPLIT, FOM_D1SIG, - FOM_D2SIG + FOM_D2SIG, + FOM_NUM_MEASUREMENTS, + FOM_REDUNDANCY, + FOM_SNR, + FOM_MEAN_INTENSITY, + FOM_STDDEV_INTENSITY, + FOM_COMPLETENESS, }; struct fom_shells @@ -61,51 +80,53 @@ struct fom_shells double *rmaxs; }; -struct fom_context -{ - enum fom_type fom; - int nshells; - int *cts; - - /* For R-factors */ - double *num; - double *den; +struct fom_context; + +extern struct fom_rejections fom_select_reflection_pairs(RefList *list1, + RefList *list2, + RefList **plist1_acc, + RefList **plist2_acc, + UnitCell *cell, + SymOpList *sym, + int anom, + double rmin_fix, + double rmax_fix, + double sigma_cutoff, + int ignore_negs, + int zero_negs, + int mul_cutoff); + +extern struct fom_rejections fom_select_reflections(RefList *list, + RefList **plist_acc, + UnitCell *cell, + SymOpList *sym, + double rmin_fix, + double rmax_fix, + double sigma_cutoff, + int ignore_negs, + int zero_negs, + int mul_cutoff); - /* For "double" R-factors */ - double *num2; - double *den2; - - /* For CCs */ - double **vec1; - double **vec2; - int *n; - int nmax; - - /* For "counting" things e.g. d1sig or d2sig */ - int *n_within; -}; - -extern int fom_select_reflections(RefList *list1, RefList *list2, - RefList *list1_acc, RefList *list2_acc, - UnitCell *cell, SymOpList *sym, - int anom, double rmin_fix, double rmax_fix, - double sigma_cutoff, int ignore_negs, - int zero_negs, int mul_cutoff); extern struct fom_context *fom_calculate(RefList *list1, RefList *list2, UnitCell *cell, struct fom_shells *shells, enum fom_type fom, int noscale, - SymOpList *sym); + const SymOpList *sym); extern struct fom_shells *fom_make_resolution_shells(double rmin, double rmax, int nshells); -extern double fom_shell_label(struct fom_shells *s, int i); +extern double fom_shell_centre(struct fom_shells *s, int i); + +extern double fom_overall_value(struct fom_context *fctx); +extern double fom_shell_value(struct fom_context *fctx, int i); -extern double fom_shell(struct fom_context *fctx, int i); +extern int fom_overall_num_reflections(struct fom_context *fctx); +extern int fom_shell_num_reflections(struct fom_context *fctx, int i); -extern double fom_overall(struct fom_context *fctx); +extern int fom_overall_num_possible(struct fom_context *fctx); +extern int fom_shell_num_possible(struct fom_context *fctx, int i); extern enum fom_type fom_type_from_string(const char *s); diff --git a/src/check_hkl.c b/src/check_hkl.c index f743d575..d95fc51f 100644 --- a/src/check_hkl.c +++ b/src/check_hkl.c @@ -45,6 +45,7 @@ #include <reflist.h> #include <reflist-utils.h> #include <cell-utils.h> +#include <fom.h> #include "version.h" @@ -390,67 +391,16 @@ static void plot_shells(RefList *list, UnitCell *cell, const SymOpList *sym, double rmin_fix, double rmax_fix, int nshells, const char *shell_file) { - unsigned long *possible; - unsigned long *measurements; - unsigned long *measured; - unsigned long *snr_measured; - double total_vol, vol_per_shell; - double *rmins; - double *rmaxs; - double *snr; - double *mean; - double *var; double rmin, rmax; - signed int h, k, l; int i; FILE *fh; - double snr_total = 0; - unsigned long nrefl = 0; - unsigned long nmeastot = 0; - unsigned long nout = 0; - unsigned long possible_tot = 0; - unsigned long nsilly = 0; - Reflection *refl; - RefListIterator *iter; - RefList *counted; - int hmax, kmax, lmax; - double ax, ay, az; - double bx, by, bz; - double cx, cy, cz; - - possible = malloc(nshells*sizeof(unsigned long)); - measurements = malloc(nshells*sizeof(unsigned long)); - measured = malloc(nshells*sizeof(unsigned long)); - snr_measured = malloc(nshells*sizeof(unsigned long)); - if ( (possible == NULL) || (measurements == NULL) - || (measured == NULL) || (snr_measured == NULL) ) { - ERROR("Couldn't allocate memory.\n"); - free(possible); - free(measurements); - free(measured); - free(snr_measured); - return; - } - - rmins = malloc(nshells*sizeof(double)); - rmaxs = malloc(nshells*sizeof(double)); - snr = malloc(nshells*sizeof(double)); - mean = malloc(nshells*sizeof(double)); - var = malloc(nshells*sizeof(double)); - if ( (rmins == NULL) || (rmaxs == NULL) || (snr == NULL) - || (mean == NULL) || (var == NULL) ) { - ERROR("Couldn't allocate memory.\n"); - free(possible); - free(measurements); - free(measured); - free(snr_measured); - free(rmins); - free(rmaxs); - free(snr); - free(mean); - free(var); - return; - } + struct fom_shells *shells; + struct fom_context *nmeas_ctx; + struct fom_context *red_ctx; + struct fom_context *snr_ctx; + struct fom_context *mean_ctx; + struct fom_context *stddev_ctx; + struct fom_context *compl_ctx; fh = fopen(shell_file, "w"); if ( fh == NULL ) { @@ -458,224 +408,72 @@ static void plot_shells(RefList *list, UnitCell *cell, const SymOpList *sym, return; } - for ( i=0; i<nshells; i++ ) { - possible[i] = 0; - measured[i] = 0; - snr_measured[i] = 0; - measurements[i] = 0; - snr[i] = 0; - var[i] = 0; - mean[i] = 0; - } - resolution_limits(list, cell, &rmin, &rmax); STATUS("1/d goes from %f to %f nm^-1\n", rmin/1e9, rmax/1e9); - /* Widen the range just a little bit */ - rmin -= 0.001e9; - rmax += 0.001e9; - /* Fixed resolution shells if needed */ if ( rmin_fix > 0.0 ) rmin = rmin_fix; if ( rmax_fix > 0.0 ) rmax = rmax_fix; - total_vol = pow(rmax, 3.0) - pow(rmin, 3.0); - vol_per_shell = total_vol / nshells; - rmins[0] = rmin; - for ( i=1; i<nshells; i++ ) { - - double r; - - r = vol_per_shell + pow(rmins[i-1], 3.0); - r = pow(r, 1.0/3.0); - - /* Shells of constant volume */ - rmaxs[i-1] = r; - rmins[i] = r; - } - rmaxs[nshells-1] = rmax; - - /* Count the number of reflections possible in each shell */ - counted = reflist_new(); - cell_get_cartesian(cell, &ax, &ay, &az, - &bx, &by, &bz, - &cx, &cy, &cz); - hmax = rmax * modulus(ax, ay, az); - kmax = rmax * modulus(bx, by, bz); - lmax = rmax * modulus(cx, cy, cz); - for ( h=-hmax; h<=hmax; h++ ) { - for ( k=-kmax; k<=kmax; k++ ) { - for ( l=-lmax; l<=lmax; l++ ) { - - double d; - signed int hs, ks, ls; - int bin; - - get_asymm(sym, h, k, l, &hs, &ks, &ls); - d = 2.0 * resolution(cell, hs, ks, ls); - - if ( forbidden_reflection(cell, h, k, l) ) continue; - - bin = -1; - for ( i=0; i<nshells; i++ ) { - if ( (d>rmins[i]) && (d<=rmaxs[i]) ) { - bin = i; - break; - } - } - if ( bin == -1 ) continue; - - if ( find_refl(counted, hs, ks, ls) != NULL ) continue; - add_refl(counted, hs, ks, ls); - - possible[bin]++; - possible_tot++; - - } - } - } - reflist_free(counted); - - /* Calculate means */ - for ( refl = first_refl(list, &iter); - refl != NULL; - refl = next_refl(refl, iter) ) - { - signed int h, k, l, hs, ks, ls; - double d, val, esd; - int bin; - int j; - - get_indices(refl, &h, &k, &l); - if ( forbidden_reflection(cell, h, k, l) ) continue; - - get_asymm(sym, h, k, l, &hs, &ks, &ls); - d = resolution(cell, hs, ks, ls) * 2.0; - val = get_intensity(refl); - esd = get_esd_intensity(refl); - - bin = -1; - for ( j=0; j<nshells; j++ ) { - if ( (d>rmins[j]) && (d<=rmaxs[j]) ) { - bin = j; - break; - } - } - if ( bin == -1 ) continue; - - measured[bin]++; - mean[bin] += get_intensity(refl); - - if ( !isfinite(val/esd) ) nsilly++; - - } - - for ( i=0; i<nshells; i++ ) { - mean[i] /= (double)measured[i]; - } - - /* Characterise the data set */ - for ( refl = first_refl(list, &iter); - refl != NULL; - refl = next_refl(refl, iter) ) - { - signed int h, k, l; - double d; - int bin; - int j; - double val, esd; - - get_indices(refl, &h, &k, &l); - if ( forbidden_reflection(cell, h, k, l) ) continue; - - d = resolution(cell, h, k, l) * 2.0; - val = get_intensity(refl); - esd = get_esd_intensity(refl); - - bin = -1; - for ( j=0; j<nshells; j++ ) { - if ( (d>rmins[j]) && (d<=rmaxs[j]) ) { - bin = j; - break; - } - } - if ( bin == -1 ) { - nout++; - continue; - } - - /* measured[bin] was done earlier */ - measurements[bin] += get_redundancy(refl); - - if ( isfinite(val/esd) ) { - snr[bin] += val / esd; - snr_total += val / esd; - snr_measured[bin]++; - } else { - nsilly++; - } - - nrefl++; - nmeastot += get_redundancy(refl); - - var[bin] += pow(val-mean[bin], 2.0); - - } + shells = fom_make_resolution_shells(rmin, rmax, nshells); STATUS("Overall values within specified resolution range:\n"); - STATUS("%li measurements in total.\n", nmeastot); - STATUS("%li reflections in total.\n", nrefl); - STATUS("%li reflections possible.\n", possible_tot); - STATUS("Overall <snr> = %f\n", snr_total/(double)nrefl); - STATUS("Overall redundancy = %f measurements/unique reflection\n", - nmeastot/(double)nrefl); - STATUS("Overall completeness = %f %%\n", 100.0*nrefl/(double)possible_tot); - if ( nout ) { - STATUS("WARNING; %li reflections outside resolution range.\n", - nout); - } - - if ( nsilly ) { - STATUS("WARNING; %li reflections had infinite or invalid values" - " of I/sigma(I).\n", nsilly); - } + nmeas_ctx = fom_calculate(list, NULL, cell, shells, + FOM_NUM_MEASUREMENTS, 0, sym); + red_ctx = fom_calculate(list, NULL, cell, shells, + FOM_REDUNDANCY, 0, sym); + snr_ctx = fom_calculate(list, NULL, cell, shells, + FOM_SNR, 0, sym); + mean_ctx = fom_calculate(list, NULL, cell, shells, + FOM_MEAN_INTENSITY, 0, sym); + stddev_ctx = fom_calculate(list, NULL, cell, shells, + FOM_STDDEV_INTENSITY, 0, sym); + compl_ctx = fom_calculate(list, NULL, cell, shells, + FOM_COMPLETENESS, 0, sym); + + STATUS("%.0f measurements in total.\n", + fom_overall_value(nmeas_ctx)); + STATUS("%li reflections in total.\n", + fom_overall_num_reflections(compl_ctx)); + STATUS("%li reflections possible.\n", + fom_overall_num_possible(compl_ctx)); + STATUS("Overall <snr> = %f\n", fom_overall_value(snr_ctx)); + STATUS("Overall redundancy = %f measurements/unique reflection\n", + fom_overall_value(red_ctx)); + STATUS("Overall completeness = %f %%\n", + 100.0*fom_overall_value(compl_ctx)); fprintf(fh, "Center 1/nm # refs Possible Compl " "Meas Red SNR Std dev Mean d(A) " "Min 1/nm Max 1/nm\n"); for ( i=0; i<nshells; i++ ) { - double cen; - cen = rmins[i] + (rmaxs[i] - rmins[i])/2.0; - fprintf(fh, "%10.3f %8li %8li %6.2f %10li %5.1f" + long int measured, possible; + + measured = fom_shell_num_reflections(compl_ctx, i); + possible = fom_shell_num_possible(compl_ctx, i); + + fprintf(fh, "%10.3f %8li %8li %6.2f %10.0f %5.1f" " %5.2f %10.2f %10.2f %8.2f %10.3f %10.3f\n", - cen*1.0e-9, - measured[i], - possible[i], - 100.0*(double)measured[i]/possible[i], - measurements[i], - (double)measurements[i]/measured[i], - snr[i]/(double)snr_measured[i], - sqrt(var[i]/measured[i]), - mean[i], (1.0/cen)*1e10, - rmins[i]*1.0e-9, rmaxs[i]*1.0e-9); + fom_shell_centre(shells, i)*1.0e-9, + measured, + possible, + 100.0*fom_shell_value(compl_ctx, i), + fom_shell_value(nmeas_ctx, i), + fom_shell_value(red_ctx, i), + fom_shell_value(snr_ctx, i), + fom_shell_value(stddev_ctx, i), + fom_shell_value(mean_ctx, i), + (1.0/fom_shell_centre(shells, i))*1e10, + shells->rmins[i]*1.0e-9, + shells->rmaxs[i]*1.0e-9); } fclose(fh); STATUS("Resolution shell information written to %s.\n", shell_file); - - free(possible); - free(measurements); - free(measured); - free(snr_measured); - free(rmins); - free(rmaxs); - free(snr); - free(mean); - free(var); } @@ -711,10 +509,7 @@ int main(int argc, char *argv[]) SymOpList *sym; RefList *raw_list; RefList *list; - Reflection *refl; - RefListIterator *iter; char *cellfile = NULL; - int rej = 0; float rmin_fix = -1.0; float rmax_fix = -1.0; float sigma_cutoff = -INFINITY; @@ -725,9 +520,8 @@ int main(int argc, char *argv[]) int ltest = 0; int ignorenegs = 0; int zeronegs = 0; - int nneg = 0; - int nres = 0; float highres, lowres; + struct fom_rejections rej; /* Long options */ const struct option longopts[] = { @@ -905,74 +699,32 @@ int main(int argc, char *argv[]) } /* Reject some reflections */ - list = reflist_new(); - for ( refl = first_refl(raw_list, &iter); - refl != NULL; - refl = next_refl(refl, iter) ) { - - signed int h, k, l; - double val, sig; - int ig = 0; - Reflection *new; - - get_indices(refl, &h, &k, &l); + rej = fom_select_reflections(raw_list, &list, cell, sym, + rmin_fix, rmax_fix, sigma_cutoff, + ignorenegs, zeronegs, 0); - val = get_intensity(refl); - sig = get_esd_intensity(refl); - - if ( val < sigma_cutoff * sig ) { - rej++; - ig = 1; - } - - if ( ignorenegs && (val < 0.0) ) { - nneg++; - ig = 1; - } - - if ( zeronegs && (val < 0.0) ) { - set_intensity(refl, 0.0); - nneg++; - } - - if ( rmin_fix > 0.0 ) { - double res = 2.0*resolution(cell, h, k, l); - if ( res < rmin_fix ) { - nres++; - continue; - } - } - - if ( rmax_fix > 0.0 ) { - double res = 2.0*resolution(cell, h, k, l); - if ( res > rmax_fix ) { - nres++; - continue; - } - } - - if ( ig ) continue; - - new = add_refl(list, h, k, l); - copy_data(new, refl); - - } STATUS("Discarded %i reflections (out of %i) with I/sigma(I) < %f\n", - rej, num_reflections(raw_list), sigma_cutoff); + rej.low_snr, num_reflections(raw_list), sigma_cutoff); reflist_free(raw_list); - if ( ignorenegs && (nneg > 0) ) { + if ( rej.negative_deleted > 0 ) { STATUS("Discarded %i reflections because they had negative " - "intensities.\n", nneg); + "intensities.\n", rej.negative_deleted); } - if ( zeronegs && (nneg > 0) ) { - STATUS("Set %i negative intensities to zero\n", nneg); + if ( rej.negative_zeroed > 0 ) { + STATUS("Set %i negative intensities to zero\n", + rej.negative_zeroed); } - if ( nres > 0 ) { + if ( rej.outside_resolution_range > 0 ) { STATUS("%i reflections rejected because they were outside the " - "resolution range.\n", nres); + "resolution range.\n", rej.outside_resolution_range); + } + + if ( rej.nan_inf_value ) { + STATUS("WARNING: %i reflections had infinite or invalid values" + " of I or sigma(I).\n", rej.nan_inf_value); } if ( wilson ) { diff --git a/src/compare_hkl.c b/src/compare_hkl.c index 2cd7207b..cb5bee38 100644 --- a/src/compare_hkl.c +++ b/src/compare_hkl.c @@ -106,55 +106,58 @@ static void do_fom(RefList *list1, RefList *list2, UnitCell *cell, switch ( fom ) { case FOM_R1I : - STATUS("Overall R1(I) = %.2f %%\n", 100.0*fom_overall(fctx)); + STATUS("Overall R1(I) = %.2f %%\n", 100.0*fom_overall_value(fctx)); break; case FOM_R1F : - STATUS("Overall R1(F) = %.2f %%\n", 100.0*fom_overall(fctx)); + STATUS("Overall R1(F) = %.2f %%\n", 100.0*fom_overall_value(fctx)); break; case FOM_R2 : - STATUS("Overall R(2) = %.2f %%\n", 100.0*fom_overall(fctx)); + STATUS("Overall R(2) = %.2f %%\n", 100.0*fom_overall_value(fctx)); break; case FOM_RSPLIT : - STATUS("Overall Rsplit = %.2f %%\n", 100.0*fom_overall(fctx)); + STATUS("Overall Rsplit = %.2f %%\n", 100.0*fom_overall_value(fctx)); break; case FOM_CC : - STATUS("Overall CC = %.7f\n", fom_overall(fctx)); + STATUS("Overall CC = %.7f\n", fom_overall_value(fctx)); break; case FOM_CCSTAR : - STATUS("Overall CC* = %.7f\n", fom_overall(fctx)); + STATUS("Overall CC* = %.7f\n", fom_overall_value(fctx)); break; case FOM_CCANO : - STATUS("Overall CCano = %.7f\n", fom_overall(fctx)); + STATUS("Overall CCano = %.7f\n", fom_overall_value(fctx)); break; case FOM_CRDANO : - STATUS("Overall CRDano = %.7f\n", fom_overall(fctx)); + STATUS("Overall CRDano = %.7f\n", fom_overall_value(fctx)); break; case FOM_RANO : - STATUS("Overall Rano = %.2f %%\n", 100.0*fom_overall(fctx)); + STATUS("Overall Rano = %.2f %%\n", 100.0*fom_overall_value(fctx)); break; case FOM_RANORSPLIT : - STATUS("Overall Rano/Rsplit = %.7f\n", fom_overall(fctx)); + STATUS("Overall Rano/Rsplit = %.7f\n", fom_overall_value(fctx)); break; case FOM_D1SIG : STATUS("Fraction of differences less than 1 sigma = %.7f %%\n", - 100.0*fom_overall(fctx)); + 100.0*fom_overall_value(fctx)); break; case FOM_D2SIG : STATUS("Fraction of differences less than 2 sigma = %.7f %%\n", - 100.0*fom_overall(fctx)); + 100.0*fom_overall_value(fctx)); break; + default : + ERROR("Unhandled figure of merit type\n"); + break; } @@ -217,14 +220,17 @@ static void do_fom(RefList *list1, RefList *list2, UnitCell *cell, fprintf(fh, "%s D<2sigma/%% nref%s\n", t1, t2); break; + default : + break; + } for ( i=0; i<nshells; i++ ) { double r, cen; - cen = fom_shell_label(shells, i); - r = fom_shell(fctx, i); + cen = fom_shell_centre(shells, i); + r = fom_shell_value(fctx, i); switch ( fom ) { @@ -235,7 +241,8 @@ static void do_fom(RefList *list1, RefList *list2, UnitCell *cell, case FOM_RANO : fprintf(fh, "%10.3f %10.2f %10i %10.2f " "%10.3f %10.3f\n", - cen*1.0e-9, r*100.0, fctx->cts[i], + cen*1.0e-9, r*100.0, + fom_shell_num_reflections(fctx, i), (1.0/cen)*1e10, shells->rmins[i]*1.0e-9, shells->rmaxs[i]*1.0e-9); @@ -247,7 +254,9 @@ static void do_fom(RefList *list1, RefList *list2, UnitCell *cell, case FOM_CRDANO : fprintf(fh, "%10.3f %10.7f %10i %10.2f " "%10.3f %10.3f\n", - cen*1.0e-9, r, fctx->cts[i], (1.0/cen)*1e10, + cen*1.0e-9, r, + fom_shell_num_reflections(fctx, i), + (1.0/cen)*1e10, shells->rmins[i]*1.0e-9, shells->rmaxs[i]*1.0e-9); break; @@ -255,7 +264,9 @@ static void do_fom(RefList *list1, RefList *list2, UnitCell *cell, case FOM_RANORSPLIT : fprintf(fh, "%10.3f %10.7f %10i %10.2f " "%10.3f %10.3f\n", - cen*1.0e-9, r, fctx->cts[i], (1.0/cen)*1e10, + cen*1.0e-9, r, + fom_shell_num_reflections(fctx, i), + (1.0/cen)*1e10, shells->rmins[i]*1.0e-9, shells->rmaxs[i]*1.0e-9); break; @@ -264,12 +275,16 @@ static void do_fom(RefList *list1, RefList *list2, UnitCell *cell, case FOM_D2SIG : fprintf(fh, "%10.3f %10.2f %10i %10.2f " "%10.3f %10.3f\n", - cen*1.0e-9, r*100.0, fctx->cts[i], + cen*1.0e-9, r*100.0, + fom_shell_num_reflections(fctx, i), (1.0/cen)*1e10, shells->rmins[i]*1.0e-9, shells->rmaxs[i]*1.0e-9); break; + default : + break; + } } @@ -311,7 +326,6 @@ int main(int argc, char *argv[]) char *sym_str_fromfile1 = NULL; char *sym_str_fromfile2 = NULL; SymOpList *sym; - int ncom; RefList *list1_acc; RefList *list2_acc; RefList *list1; @@ -332,6 +346,7 @@ int main(int argc, char *argv[]) float highres, lowres; int mul_cutoff = 0; int anom; + struct fom_rejections rej; /* Long options */ const struct option longopts[] = { @@ -478,7 +493,7 @@ int main(int argc, char *argv[]) " intensities.\n"); ERROR("Please try again with --ignore-negs or" " --zero-negs.\n"); - exit(1); + return 1; case FOM_R2 : case FOM_R1I : @@ -492,6 +507,10 @@ int main(int argc, char *argv[]) case FOM_D1SIG : case FOM_D2SIG : break; + + default : + ERROR("Unhandled figure of merit!\n"); + return 1; } } @@ -579,6 +598,10 @@ int main(int argc, char *argv[]) " try again using a non-centrosymmetric point" " group for '-y'.\n"); return 1; + + default : + ERROR("Unhandled figure of merit type!\n"); + return 1; } } @@ -622,21 +645,60 @@ int main(int argc, char *argv[]) reflist_free(list1_raw); reflist_free(list2_raw); - list1_acc = reflist_new(); - list2_acc = reflist_new(); anom = ( (fom == FOM_CCANO) || (fom == FOM_CRDANO) || (fom == FOM_RANO) || (fom == FOM_RANORSPLIT) ); - ncom = fom_select_reflections(list1, list2, list1_acc, list2_acc, - cell, sym, - anom, rmin_fix, rmax_fix, sigma_cutoff, - config_ignorenegs, config_zeronegs, - mul_cutoff); + rej = fom_select_reflection_pairs(list1, list2, &list1_acc, &list2_acc, + cell, sym, + anom, rmin_fix, rmax_fix, sigma_cutoff, + config_ignorenegs, config_zeronegs, + mul_cutoff); reflist_free(list1); reflist_free(list2); gsl_set_error_handler_off(); - STATUS("%i reflection pairs accepted.\n", ncom); + if ( rej.low_snr > 0 ) { + STATUS("Discarded %i reflection pairs because either or both" + " versions had I/sigma(I) < %f.\n", + rej.low_snr, sigma_cutoff); + } + + if ( rej.negative_deleted > 0 ) { + STATUS("Discarded %i reflection pairs because either or both" + " versions had negative intensities.\n", + rej.negative_deleted); + } + + if ( rej.negative_zeroed > 0 ) { + STATUS("For %i reflection pairs, either or both versions had" + " negative intensities which were set to zero.\n", + rej.negative_zeroed); + } + + if ( rej.few_measurements > 0 ) { + STATUS("%i reflection pairs rejected because either or both" + " versions had too few measurements.\n", + rej.few_measurements); + } + + if ( rej.outside_resolution_range > 0 ) { + STATUS("%i reflection pairs rejected because either or both" + " versions were outside the resolution range.\n", + rej.outside_resolution_range); + } + + if ( rej.no_bijvoet > 0 ) { + STATUS("%i reflection pairs rejected because either or both" + " versions did not have Bijvoet partners.\n", + rej.no_bijvoet); + } + + if ( rej.centric > 0 ) { + STATUS("%i reflection pairs rejected because they were" + " centric.\n", rej.centric); + } + + STATUS("%i reflection pairs accepted.\n", rej.common); resolution_limits(list1_acc, cell, &rmin, &rmax); resolution_limits(list2_acc, cell, &rmin, &rmax); |