aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--libcrystfel/src/fom.c612
-rw-r--r--libcrystfel/src/fom.h89
-rw-r--r--src/check_hkl.c384
-rw-r--r--src/compare_hkl.c118
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);
generated by cgit v1.2.3 (git 2.25.1) at 2024-10-01 13:27:35 +0000