diff options
-rw-r--r-- | libcrystfel/CMakeLists.txt | 2 | ||||
-rw-r--r-- | libcrystfel/meson.build | 4 | ||||
-rw-r--r-- | libcrystfel/src/fom.c | 789 | ||||
-rw-r--r-- | libcrystfel/src/fom.h | 109 | ||||
-rw-r--r-- | src/compare_hkl.c | 921 |
5 files changed, 970 insertions, 855 deletions
diff --git a/libcrystfel/CMakeLists.txt b/libcrystfel/CMakeLists.txt index 3ab45b30..a199a68c 100644 --- a/libcrystfel/CMakeLists.txt +++ b/libcrystfel/CMakeLists.txt @@ -59,6 +59,7 @@ set(LIBCRYSTFEL_SOURCES src/detgeom.c src/image-cbf.c src/image-hdf5.c + src/fom.c ${BISON_symopp_OUTPUTS} ${FLEX_symopl_OUTPUTS} src/indexers/dirax.c @@ -97,6 +98,7 @@ set(LIBCRYSTFEL_HEADERS src/colscale.h src/detgeom.h src/image-msgpack.h + src/fom.h ) if (MSGPACK_FOUND) diff --git a/libcrystfel/meson.build b/libcrystfel/meson.build index 244fc4cc..48822d60 100644 --- a/libcrystfel/meson.build +++ b/libcrystfel/meson.build @@ -94,6 +94,7 @@ libcrystfel_sources = ['src/image.c', 'src/datatemplate.c', 'src/colscale.c', 'src/detgeom.c', + 'src/fom.c', 'src/image-cbf.c', 'src/image-hdf5.c', 'src/indexers/dirax.c', @@ -154,7 +155,8 @@ install_headers(['src/reflist.h', 'src/datatemplate.h', 'src/colscale.h', 'src/detgeom.h', - 'src/image-msgpack.h'], + 'src/image-msgpack.h', + 'src/fom.h'], subdir: 'crystfel') # API documentation (Doxygen) diff --git a/libcrystfel/src/fom.c b/libcrystfel/src/fom.c new file mode 100644 index 00000000..1296b53f --- /dev/null +++ b/libcrystfel/src/fom.c @@ -0,0 +1,789 @@ +/* + * fom.c + * + * Figure of merit calculation + * + * Copyright © 2012-2021 Deutsches Elektronen-Synchrotron DESY, + * a research centre of the Helmholtz Association. + * + * Authors: + * 2010-2021 Thomas White <taw@physics.org> + * 2013 Lorenzo Galli <lorenzo.galli@desy.de> + * + * This file is part of CrystFEL. + * + * CrystFEL is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * CrystFEL is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with CrystFEL. If not, see <http://www.gnu.org/licenses/>. + * + */ + +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif + +#include <gsl/gsl_errno.h> +#include <gsl/gsl_statistics.h> +#include <gsl/gsl_fit.h> +#include <assert.h> + +#include "utils.h" +#include "fom.h" +#include "cell.h" +#include "cell-utils.h" +#include "reflist.h" +#include "reflist-utils.h" + +enum fom get_fom(const char *s) +{ + if ( strcasecmp(s, "r1i") == 0 ) return FOM_R1I; + if ( strcasecmp(s, "r1f") == 0 ) return FOM_R1F; + if ( strcasecmp(s, "r2") == 0 ) return FOM_R2; + if ( strcasecmp(s, "rsplit") == 0 ) return FOM_RSPLIT; + if ( strcasecmp(s, "cc") == 0 ) return FOM_CC; + if ( strcasecmp(s, "ccstar") == 0 ) return FOM_CCSTAR; + if ( strcasecmp(s, "ccano") == 0 ) return FOM_CCANO; + if ( strcasecmp(s, "crdano") == 0 ) return FOM_CRDANO; + if ( strcasecmp(s, "rano") == 0 ) return FOM_RANO; + if ( strcasecmp(s, "rano/rsplit") == 0 ) return FOM_RANORSPLIT; + if ( strcasecmp(s, "d1sig") == 0 ) return FOM_D1SIG; + if ( strcasecmp(s, "d2sig") == 0 ) return FOM_D2SIG; + + ERROR("Unknown figure of merit '%s'.\n", s); + exit(1); +} + + +static struct fom_context *init_fom(enum fom fom, int nmax, int nshells) +{ + struct fom_context *fctx; + int i; + + fctx = malloc(sizeof(struct fom_context)); + if ( fctx == NULL ) return NULL; + + fctx->fom = fom; + fctx->nshells = nshells; + fctx->cts = malloc(nshells*sizeof(int)); + for ( i=0; i<nshells; i++ ) { + fctx->cts[i] = 0; + } + + switch ( fctx->fom ) { + + case FOM_RANORSPLIT : + fctx->num2 = malloc(nshells*sizeof(double)); + fctx->den2 = malloc(nshells*sizeof(double)); + if ( (fctx->num2 == NULL) || (fctx->den2 == NULL) ) return NULL; + for ( i=0; i<nshells; i++ ) { + fctx->num2[i] = 0.0; + fctx->den2[i] = 0.0; + } + /* Intentional fall-through (no break) */ + + case FOM_R1I : + case FOM_R1F : + case FOM_R2 : + case FOM_RSPLIT : + case FOM_RANO : + fctx->num = malloc(nshells*sizeof(double)); + fctx->den = malloc(nshells*sizeof(double)); + if ( (fctx->num == NULL) || (fctx->den == NULL) ) return NULL; + for ( i=0; i<nshells; i++ ) { + fctx->num[i] = 0.0; + fctx->den[i] = 0.0; + } + break; + + case FOM_CC : + case FOM_CCSTAR : + case FOM_CCANO : + case FOM_CRDANO : + fctx->vec1 = malloc(nshells*sizeof(double *)); + fctx->vec2 = malloc(nshells*sizeof(double *)); + if ( (fctx->vec1 == NULL) || (fctx->vec2 == NULL) ) return NULL; + for ( i=0; i<nshells; i++ ) { + fctx->vec1[i] = malloc(nmax*sizeof(double)); + if ( fctx->vec1[i] == NULL ) return NULL; + fctx->vec2[i] = malloc(nmax*sizeof(double)); + if ( fctx->vec2[i] == NULL ) return NULL; + fctx->n = malloc(nshells*sizeof(int)); + if ( fctx->n == NULL ) return NULL; + } + for ( i=0; i<nshells; i++ ) { + fctx->n[i] = 0; + } + fctx->nmax = nmax; + break; + + case FOM_D1SIG : + case FOM_D2SIG : + fctx->n_within = malloc(nshells*sizeof(int)); + if ( fctx->n_within == NULL ) return NULL; + for ( i=0; i<nshells; i++ ) { + fctx->n_within[i] = 0; + } + break; + + } + + return fctx; +} + + +static void add_to_fom(struct fom_context *fctx, double i1, double i2, + double i1bij, double i2bij, double sig1, double sig2, + int bin) +{ + double f1, f2; + double im, imbij; + + fctx->cts[bin]++; + + /* Negative intensities have already been weeded out. */ + f1 = sqrt(i1); + f2 = sqrt(i2); + + switch ( fctx->fom ) { + + case FOM_R1I : + fctx->num[bin] += fabs(i1 - i2); + fctx->den[bin] += i1; + break; + + case FOM_R1F : + fctx->num[bin] += fabs(f1 - f2); + fctx->den[bin] += f1; + break; + + case FOM_R2 : + fctx->num[bin] += pow(i1 - i2, 2.0); + fctx->den[bin] += pow(i1, 2.0); + break; + + case FOM_RSPLIT : + fctx->num[bin] += fabs(i1 - i2); + fctx->den[bin] += i1 + i2; + break; + + case FOM_CC : + case FOM_CCSTAR : + assert(fctx->n[bin] < fctx->nmax); + fctx->vec1[bin][fctx->n[bin]] = i1; + fctx->vec2[bin][fctx->n[bin]] = i2; + fctx->n[bin]++; + break; + + case FOM_CCANO : + case FOM_CRDANO : + assert(fctx->n[bin] < fctx->nmax); + fctx->vec1[bin][fctx->n[bin]] = i1 - i1bij; + fctx->vec2[bin][fctx->n[bin]] = i2 - i2bij; + fctx->n[bin]++; + break; + + case FOM_RANORSPLIT : + fctx->num2[bin] += fabs(i1 - i2); + fctx->den2[bin] += i1 + i2; + /* Intentional fall-through (no break) */ + + case FOM_RANO : + im = (i1 + i2)/2.0; + imbij = (i1bij + i2bij)/2.0; + fctx->num[bin] += fabs(im - imbij); + fctx->den[bin] += im + imbij; + break; + + case FOM_D1SIG : + if ( fabs(i1-i2) < sqrt(sig1*sig1 + sig2*sig2) ) { + fctx->n_within[bin]++; + } + break; + + case FOM_D2SIG : + if ( fabs(i1-i2) < 2.0*sqrt(sig1*sig1 + sig2*sig2) ) { + fctx->n_within[bin]++; + } + break; + + } +} + + +double fom_overall(struct fom_context *fctx) +{ + double overall_num = INFINITY; + double overall_den = 0.0; + double overall_num2 = INFINITY; + double overall_den2 = 0.0; + int i; + double *overall_vec1; + double *overall_vec2; + int overall_n; + double *overall_along_diagonal; + double *overall_perpend_diagonal; + double variance_signal; + double variance_error; + double cc = INFINITY; + + switch ( fctx->fom ) { + + case FOM_R1I : + case FOM_R1F : + case FOM_R2 : + case FOM_RSPLIT : + case FOM_RANO : + overall_num = 0.0; + overall_den = 0.0; + for ( i=0; i<fctx->nshells; i++ ) { + overall_num += fctx->num[i]; + overall_den += fctx->den[i]; + } + break; + + case FOM_RANORSPLIT : + overall_num = 0.0; + overall_den = 0.0; + for ( i=0; i<fctx->nshells; i++ ) { + overall_num += fctx->num[i]; + overall_den += fctx->den[i]; + } + overall_num2 = 0.0; + overall_den2 = 0.0; + for ( i=0; i<fctx->nshells; i++ ) { + overall_num2 += fctx->num2[i]; + overall_den2 += fctx->den2[i]; + } + break; + + case FOM_CC : + case FOM_CCSTAR : + case FOM_CCANO : + overall_vec1 = malloc(fctx->nmax*sizeof(double)); + overall_vec2 = malloc(fctx->nmax*sizeof(double)); + overall_n = 0; + for ( i=0; i<fctx->nshells; i++ ) { + int j; + for ( j=0; j<fctx->n[i]; j++ ) { + overall_vec1[overall_n] = fctx->vec1[i][j]; + overall_vec2[overall_n] = fctx->vec2[i][j]; + overall_n++; + } + } + cc = gsl_stats_correlation(overall_vec1, 1, overall_vec2, 1, + overall_n); + free(overall_vec1); + free(overall_vec2); + break; + + case FOM_CRDANO : + overall_along_diagonal = malloc(fctx->nmax*sizeof(double)); + overall_perpend_diagonal = malloc(fctx->nmax*sizeof(double)); + overall_n = 0; + for ( i=0; i<fctx->nshells; i++ ) { + int j; + for ( j=0; j<fctx->n[i]; j++ ) { + overall_along_diagonal[overall_n] = + ( fctx->vec1[i][j] + fctx->vec2[i][j] ) + / sqrt(2.0); + overall_perpend_diagonal[overall_n] = + ( fctx->vec1[i][j] - fctx->vec2[i][j] ) + / sqrt(2.0); + overall_n++; + } + } + variance_signal = gsl_stats_variance_m(overall_along_diagonal, + 1, overall_n, 0.0); + variance_error = gsl_stats_variance_m(overall_perpend_diagonal, + 1, overall_n, 0.0); + cc = sqrt(variance_signal / variance_error ); + + free(overall_along_diagonal); + free(overall_perpend_diagonal); + break; + + case FOM_D1SIG : + case FOM_D2SIG : + overall_num = 0.0; + overall_den = 0.0; + for ( i=0; i<fctx->nshells; i++ ) { + overall_num += fctx->n_within[i]; + overall_den += fctx->cts[i]; + } + break; + + } + + switch ( fctx->fom ) { + + case FOM_R1I : + case FOM_R1F : + return overall_num/overall_den; + + case FOM_R2 : + return sqrt(overall_num/overall_den); + + case FOM_RSPLIT : + return 2.0*(overall_num/overall_den) / sqrt(2.0); + + case FOM_CC : + case FOM_CCANO : + case FOM_CRDANO : + return cc; + + case FOM_CCSTAR : + return sqrt((2.0*cc)/(1.0+cc)); + + case FOM_RANO : + return 2.0*(overall_num/overall_den); + + case FOM_RANORSPLIT : + return (2.0*(overall_num/overall_den)) / + (2.0*(overall_num2/overall_den2) / sqrt(2.0)); + + case FOM_D1SIG : + case FOM_D2SIG : + return overall_num/overall_den; + + } + + ERROR("This point is never reached.\n"); + abort(); +} + + +double fom_shell(struct fom_context *fctx, int i) +{ + double cc; + int j; + double variance_signal; + double variance_error; + double *along_diagonal; + double *perpend_diagonal; + + switch ( fctx->fom ) { + + case FOM_R1I : + case FOM_R1F : + return fctx->num[i]/fctx->den[i]; + + case FOM_R2 : + return sqrt(fctx->num[i]/fctx->den[i]); + + case FOM_RSPLIT : + return 2.0*(fctx->num[i]/fctx->den[i]) / sqrt(2.0); + + case FOM_CC : + case FOM_CCANO : + return gsl_stats_correlation(fctx->vec1[i], 1, fctx->vec2[i], 1, + fctx->n[i]); + + case FOM_CCSTAR : + cc = gsl_stats_correlation(fctx->vec1[i], 1, fctx->vec2[i], 1, + fctx->n[i]); + return sqrt((2.0*cc)/(1.0+cc)); + + case FOM_RANO : + return 2.0 * fctx->num[i]/fctx->den[i]; + + case FOM_RANORSPLIT : + return (2.0*fctx->num[i]/fctx->den[i]) / + (2.0*(fctx->num2[i]/fctx->den2[i]) / sqrt(2.0)); + + case FOM_CRDANO : + along_diagonal = malloc(fctx->n[i] * sizeof(double)); + perpend_diagonal = malloc(fctx->n[i] * sizeof(double)); + for ( j=0; j<fctx->n[i]; j++ ) { + along_diagonal[j] = ( fctx->vec1[i][j] + + fctx->vec2[i][j] ) / sqrt(2.0); + perpend_diagonal[j] = ( fctx->vec1[i][j] - + fctx->vec2[i][j] ) / sqrt(2.0); + } + variance_signal = gsl_stats_variance_m(along_diagonal, 1, + fctx->n[i], 0.0); + variance_error = gsl_stats_variance_m(perpend_diagonal, 1, + fctx->n[i], 0.0); + free(along_diagonal); + free(perpend_diagonal); + return sqrt(variance_signal / variance_error); + + case FOM_D1SIG : + case FOM_D2SIG : + return (double)fctx->n_within[i] / fctx->cts[i]; + + } + + ERROR("This point is never reached.\n"); + abort(); +} + + +struct shells *make_intensity_shells(double min_I, double max_I, int nshells) +{ + struct shells *s; + int i; + + if ( min_I >= max_I ) { + ERROR("Invalid intensity range.\n"); + return NULL; + } + + /* Adjust minimum and maximum intensities to get the most densely + * populated part of the reflections */ + max_I = min_I + (max_I-min_I)/5000.0; + + s = malloc(sizeof(struct shells)); + if ( s == NULL ) return NULL; + + s->rmins = malloc(nshells*sizeof(double)); + s->rmaxs = malloc(nshells*sizeof(double)); + + if ( (s->rmins==NULL) || (s->rmaxs==NULL) ) { + ERROR("Couldn't allocate memory for shells.\n"); + free(s); + return NULL; + } + + s->config_intshells = 1; + s->nshells = nshells; + + for ( i=0; i<nshells; i++ ) { + + s->rmins[i] = min_I + i*(max_I - min_I)/nshells;; + s->rmaxs[i] = min_I + (i+1)*(max_I - min_I)/nshells;; + + } + + return s; +} + + +struct shells *make_resolution_shells(double rmin, double rmax, int nshells) +{ + struct shells *s; + double total_vol, vol_per_shell; + int i; + + s = malloc(sizeof(struct shells)); + if ( s == NULL ) return NULL; + + s->rmins = malloc(nshells*sizeof(double)); + s->rmaxs = malloc(nshells*sizeof(double)); + + if ( (s->rmins==NULL) || (s->rmaxs==NULL) ) { + ERROR("Couldn't allocate memory for resolution shells.\n"); + free(s); + return NULL; + } + + s->config_intshells = 0; + s->nshells = nshells; + + total_vol = pow(rmax, 3.0) - pow(rmin, 3.0); + vol_per_shell = total_vol / nshells; + s->rmins[0] = rmin; + for ( i=1; i<nshells; i++ ) { + + double r; + + r = vol_per_shell + pow(s->rmins[i-1], 3.0); + r = pow(r, 1.0/3.0); + + /* Shells of constant volume */ + s->rmaxs[i-1] = r; + s->rmins[i] = r; + + } + s->rmaxs[nshells-1] = rmax; + + return s; +} + + +double shell_label(struct shells *s, int i) +{ + if ( s->config_intshells ) { + return (i+0.5) / s->nshells; + } else { + return s->rmins[i] + (s->rmaxs[i] - s->rmins[i])/2.0; + } +} + + +static int get_bin(struct shells *s, Reflection *refl, UnitCell *cell) +{ + if ( s->config_intshells ) { + + double intensity; + int bin, j; + + intensity = get_intensity(refl); + + bin = -1; + for ( j=0; j<s->nshells; j++ ) { + if ( (intensity>s->rmins[j]) + && (intensity<=s->rmaxs[j]) ) + { + bin = j; + break; + } + } + + return bin; + + } else { + + double d; + int bin, j; + signed int h, k, l; + + get_indices(refl, &h, &k, &l); + d = 2.0 * resolution(cell, h, k, l); + + bin = -1; + for ( j=0; j<s->nshells; j++ ) { + if ( (d>s->rmins[j]) && (d<=s->rmaxs[j]) ) { + bin = j; + break; + } + } + + /* Allow for slight rounding errors */ + if ( (bin == -1) && (d <= s->rmins[0]) ) bin = 0; + if ( (bin == -1) && (d >= s->rmaxs[s->nshells-1]) ) bin = 0; + assert(bin != -1); + + return bin; + + } +} + + +static int wilson_scale(RefList *list1, RefList *list2, UnitCell *cell) +{ + Reflection *refl1; + Reflection *refl2; + RefListIterator *iter; + int max_n = 256; + int n = 0; + double *x; + double *y; + int r; + double G, B; + double c0, c1, cov00, cov01, cov11, chisq; + + x = malloc(max_n*sizeof(double)); + y = malloc(max_n*sizeof(double)); + if ( (x==NULL) || (y==NULL) ) { + ERROR("Failed to allocate memory for scaling.\n"); + return 1; + } + + for ( refl1 = first_refl(list1, &iter); + refl1 != NULL; + refl1 = next_refl(refl1, iter) ) + { + signed int h, k, l; + double Ih1, Ih2; + double res; + + get_indices(refl1, &h, &k, &l); + res = resolution(cell, h, k, l); + + refl2 = find_refl(list2, h, k, l); + assert(refl2 != NULL); + + Ih1 = get_intensity(refl1); + Ih2 = get_intensity(refl2); + + if ( (Ih1 <= 0.0) || (Ih2 <= 0.0) ) continue; + if ( isnan(Ih1) || isinf(Ih1) ) continue; + if ( isnan(Ih2) || isinf(Ih2) ) continue; + + if ( n == max_n ) { + max_n *= 2; + x = realloc(x, max_n*sizeof(double)); + y = realloc(y, max_n*sizeof(double)); + if ( (x==NULL) || (y==NULL) ) { + ERROR("Failed to allocate memory for scaling.\n"); + return 1; + } + } + + x[n] = res*res; + y[n] = log(Ih1/Ih2); + n++; + + } + + if ( n < 2 ) { + ERROR("Not enough reflections for scaling\n"); + return 1; + } + + r = gsl_fit_linear(x, 1, y, 1, n, &c0, &c1, + &cov00, &cov01, &cov11, &chisq); + + if ( r ) { + ERROR("Scaling failed.\n"); + return 1; + } + + G = exp(c0); + B = c1/2.0; + + STATUS("Relative scale factor = %f, relative B factor = %f A^2\n", + G, B*1e20); + STATUS("A scale factor greater than 1 means that the second reflection " + "list is weaker than the first.\n"); + STATUS("A positive relative B factor means that the second reflection " + "list falls off with resolution more quickly than the first.\n"); + + free(x); + free(y); + + /* Apply the scaling factor */ + for ( refl2 = first_refl(list2, &iter); + refl2 != NULL; + refl2 = next_refl(refl2, iter) ) + { + signed int h, k, l; + double res; + double corr; + + get_indices(refl2, &h, &k, &l); + res = resolution(cell, h, k, l); + + corr = G * exp(2.0*B*res*res); + set_intensity(refl2, get_intensity(refl2)*corr); + set_esd_intensity(refl2, get_esd_intensity(refl2)*corr); + + } + return 0; +} + + + +struct fom_context *fom_calculate(RefList *list1, RefList *list2, UnitCell *cell, + struct shells *shells, enum fom fom, + int noscale, SymOpList *sym) +{ + Reflection *refl1; + RefListIterator *iter; + struct fom_context *fctx; + int n_out; + + fctx = init_fom(fom, num_reflections(list1), shells->nshells); + + if ( fctx == NULL ) { + ERROR("Couldn't allocate memory for resolution shells.\n"); + return NULL; + } + + 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); + } + + 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; + + get_indices(refl1, &h, &k, &l); + + refl2 = find_refl(list2, h, k, l); + if ( refl2 == NULL ) continue; + + bin = get_bin(shells, refl1, cell); + if ( bin == -1 ) { + n_out++; + continue; + } + + i1 = get_intensity(refl1); + i2 = get_intensity(refl2); + sig1 = get_esd_intensity(refl1); + sig2 = get_esd_intensity(refl2); + + 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, + &hb, &kb, &lb) ) + { + refl1_bij = find_refl(list1, hb, kb, lb); + } + + if ( find_equiv_in_list(list2, -h, -k, -l, sym, + &hb, &kb, &lb) ) + { + refl2_bij = find_refl(list2, hb, kb, lb); + } + + /* Each reflection must only be counted once, whether + * we are visiting it now as "normal" or "bij" */ + if ( get_flag(refl1) ) continue; + assert(!get_flag(refl2)); + set_flag(refl1, 1); + set_flag(refl1_bij, 1); + set_flag(refl2, 1); + set_flag(refl2_bij, 1); + + 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); + } + if ( n_out) { + ERROR("WARNING: %i reflection pairs outside range.\n", n_out); + } + + return fctx; +} + diff --git a/libcrystfel/src/fom.h b/libcrystfel/src/fom.h new file mode 100644 index 00000000..968278b4 --- /dev/null +++ b/libcrystfel/src/fom.h @@ -0,0 +1,109 @@ +/* + * fom.h + * + * Figure of merit calculation + * + * Copyright © 2012-2021 Deutsches Elektronen-Synchrotron DESY, + * a research centre of the Helmholtz Association. + * + * Authors: + * 2010-2021 Thomas White <taw@physics.org> + * 2013 Lorenzo Galli <lorenzo.galli@desy.de> + * + * This file is part of CrystFEL. + * + * CrystFEL is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * CrystFEL is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with CrystFEL. If not, see <http://www.gnu.org/licenses/>. + * + */ + +#ifndef FOM_H +#define FOM_H + +/** + * \file fom.h + * Figure of merit calculation + */ + +#include <reflist.h> +#include <symmetry.h> + +enum fom +{ + FOM_R1I, + FOM_R1F, + FOM_R2, + FOM_RSPLIT, + FOM_CC, + FOM_CCSTAR, + FOM_CCANO, + FOM_CRDANO, + FOM_RANO, + FOM_RANORSPLIT, + FOM_D1SIG, + FOM_D2SIG +}; + +struct shells +{ + int config_intshells; + int nshells; + double *rmins; + double *rmaxs; +}; + +struct fom_context +{ + enum fom 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; +}; + + +extern struct fom_context *fom_calculate(RefList *list1, RefList *list2, + UnitCell *cell, struct shells *shells, + enum fom fom, int noscale, + SymOpList *sym); + +extern struct shells *make_resolution_shells(double rmin, double rmax, + int nshells); + +extern struct shells *make_intensity_shells(double min_I, double max_I, + int nshells); + +extern double shell_label(struct shells *s, int i); + +extern double fom_shell(struct fom_context *fctx, int i); + +extern double fom_overall(struct fom_context *fctx); + +extern enum fom get_fom(const char *s); + +#endif /* FOM */ diff --git a/src/compare_hkl.c b/src/compare_hkl.c index 2a89596e..d4efcd05 100644 --- a/src/compare_hkl.c +++ b/src/compare_hkl.c @@ -39,55 +39,16 @@ #include <unistd.h> #include <getopt.h> #include <assert.h> -#include <gsl/gsl_errno.h> -#include <gsl/gsl_statistics.h> -#include <gsl/gsl_fit.h> #include <utils.h> #include <symmetry.h> #include <reflist-utils.h> #include <cell-utils.h> +#include <fom.h> #include "version.h" -enum fom -{ - FOM_R1I, - FOM_R1F, - FOM_R2, - FOM_RSPLIT, - FOM_CC, - FOM_CCSTAR, - FOM_CCANO, - FOM_CRDANO, - FOM_RANO, - FOM_RANORSPLIT, - FOM_D1SIG, - FOM_D2SIG -}; - - -static enum fom get_fom(const char *s) -{ - if ( strcasecmp(s, "r1i") == 0 ) return FOM_R1I; - if ( strcasecmp(s, "r1f") == 0 ) return FOM_R1F; - if ( strcasecmp(s, "r2") == 0 ) return FOM_R2; - if ( strcasecmp(s, "rsplit") == 0 ) return FOM_RSPLIT; - if ( strcasecmp(s, "cc") == 0 ) return FOM_CC; - if ( strcasecmp(s, "ccstar") == 0 ) return FOM_CCSTAR; - if ( strcasecmp(s, "ccano") == 0 ) return FOM_CCANO; - if ( strcasecmp(s, "crdano") == 0 ) return FOM_CRDANO; - if ( strcasecmp(s, "rano") == 0 ) return FOM_RANO; - if ( strcasecmp(s, "rano/rsplit") == 0 ) return FOM_RANORSPLIT; - if ( strcasecmp(s, "d1sig") == 0 ) return FOM_D1SIG; - if ( strcasecmp(s, "d2sig") == 0 ) return FOM_D2SIG; - - ERROR("Unknown figure of merit '%s'.\n", s); - exit(1); -} - - static void show_help(const char *s) { printf("Syntax: %s [options] <file1.hkl> <file2.hkl>\n\n", s); @@ -121,684 +82,23 @@ static void show_help(const char *s) } -struct fom_context -{ - enum fom 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; -}; - - -static struct fom_context *init_fom(enum fom fom, int nmax, int nshells) -{ - struct fom_context *fctx; - int i; - - fctx = malloc(sizeof(struct fom_context)); - if ( fctx == NULL ) return NULL; - - fctx->fom = fom; - fctx->nshells = nshells; - fctx->cts = malloc(nshells*sizeof(int)); - for ( i=0; i<nshells; i++ ) { - fctx->cts[i] = 0; - } - - switch ( fctx->fom ) { - - case FOM_RANORSPLIT : - fctx->num2 = malloc(nshells*sizeof(double)); - fctx->den2 = malloc(nshells*sizeof(double)); - if ( (fctx->num2 == NULL) || (fctx->den2 == NULL) ) return NULL; - for ( i=0; i<nshells; i++ ) { - fctx->num2[i] = 0.0; - fctx->den2[i] = 0.0; - } - /* Intentional fall-through (no break) */ - - case FOM_R1I : - case FOM_R1F : - case FOM_R2 : - case FOM_RSPLIT : - case FOM_RANO : - fctx->num = malloc(nshells*sizeof(double)); - fctx->den = malloc(nshells*sizeof(double)); - if ( (fctx->num == NULL) || (fctx->den == NULL) ) return NULL; - for ( i=0; i<nshells; i++ ) { - fctx->num[i] = 0.0; - fctx->den[i] = 0.0; - } - break; - - case FOM_CC : - case FOM_CCSTAR : - case FOM_CCANO : - case FOM_CRDANO : - fctx->vec1 = malloc(nshells*sizeof(double *)); - fctx->vec2 = malloc(nshells*sizeof(double *)); - if ( (fctx->vec1 == NULL) || (fctx->vec2 == NULL) ) return NULL; - for ( i=0; i<nshells; i++ ) { - fctx->vec1[i] = malloc(nmax*sizeof(double)); - if ( fctx->vec1[i] == NULL ) return NULL; - fctx->vec2[i] = malloc(nmax*sizeof(double)); - if ( fctx->vec2[i] == NULL ) return NULL; - fctx->n = malloc(nshells*sizeof(int)); - if ( fctx->n == NULL ) return NULL; - } - for ( i=0; i<nshells; i++ ) { - fctx->n[i] = 0; - } - fctx->nmax = nmax; - break; - - case FOM_D1SIG : - case FOM_D2SIG : - fctx->n_within = malloc(nshells*sizeof(int)); - if ( fctx->n_within == NULL ) return NULL; - for ( i=0; i<nshells; i++ ) { - fctx->n_within[i] = 0; - } - break; - - } - - return fctx; -} - - -static void add_to_fom(struct fom_context *fctx, double i1, double i2, - double i1bij, double i2bij, double sig1, double sig2, - int bin) -{ - double f1, f2; - double im, imbij; - - fctx->cts[bin]++; - - /* Negative intensities have already been weeded out. */ - f1 = sqrt(i1); - f2 = sqrt(i2); - - switch ( fctx->fom ) { - - case FOM_R1I : - fctx->num[bin] += fabs(i1 - i2); - fctx->den[bin] += i1; - break; - - case FOM_R1F : - fctx->num[bin] += fabs(f1 - f2); - fctx->den[bin] += f1; - break; - - case FOM_R2 : - fctx->num[bin] += pow(i1 - i2, 2.0); - fctx->den[bin] += pow(i1, 2.0); - break; - - case FOM_RSPLIT : - fctx->num[bin] += fabs(i1 - i2); - fctx->den[bin] += i1 + i2; - break; - - case FOM_CC : - case FOM_CCSTAR : - assert(fctx->n[bin] < fctx->nmax); - fctx->vec1[bin][fctx->n[bin]] = i1; - fctx->vec2[bin][fctx->n[bin]] = i2; - fctx->n[bin]++; - break; - - case FOM_CCANO : - case FOM_CRDANO : - assert(fctx->n[bin] < fctx->nmax); - fctx->vec1[bin][fctx->n[bin]] = i1 - i1bij; - fctx->vec2[bin][fctx->n[bin]] = i2 - i2bij; - fctx->n[bin]++; - break; - - case FOM_RANORSPLIT : - fctx->num2[bin] += fabs(i1 - i2); - fctx->den2[bin] += i1 + i2; - /* Intentional fall-through (no break) */ - - case FOM_RANO : - im = (i1 + i2)/2.0; - imbij = (i1bij + i2bij)/2.0; - fctx->num[bin] += fabs(im - imbij); - fctx->den[bin] += im + imbij; - break; - - case FOM_D1SIG : - if ( fabs(i1-i2) < sqrt(sig1*sig1 + sig2*sig2) ) { - fctx->n_within[bin]++; - } - break; - - case FOM_D2SIG : - if ( fabs(i1-i2) < 2.0*sqrt(sig1*sig1 + sig2*sig2) ) { - fctx->n_within[bin]++; - } - break; - - } -} - - -static double fom_overall(struct fom_context *fctx) -{ - double overall_num = INFINITY; - double overall_den = 0.0; - double overall_num2 = INFINITY; - double overall_den2 = 0.0; - int i; - double *overall_vec1; - double *overall_vec2; - int overall_n; - double *overall_along_diagonal; - double *overall_perpend_diagonal; - double variance_signal; - double variance_error; - double cc = INFINITY; - - switch ( fctx->fom ) { - - case FOM_R1I : - case FOM_R1F : - case FOM_R2 : - case FOM_RSPLIT : - case FOM_RANO : - overall_num = 0.0; - overall_den = 0.0; - for ( i=0; i<fctx->nshells; i++ ) { - overall_num += fctx->num[i]; - overall_den += fctx->den[i]; - } - break; - - case FOM_RANORSPLIT : - overall_num = 0.0; - overall_den = 0.0; - for ( i=0; i<fctx->nshells; i++ ) { - overall_num += fctx->num[i]; - overall_den += fctx->den[i]; - } - overall_num2 = 0.0; - overall_den2 = 0.0; - for ( i=0; i<fctx->nshells; i++ ) { - overall_num2 += fctx->num2[i]; - overall_den2 += fctx->den2[i]; - } - break; - - case FOM_CC : - case FOM_CCSTAR : - case FOM_CCANO : - overall_vec1 = malloc(fctx->nmax*sizeof(double)); - overall_vec2 = malloc(fctx->nmax*sizeof(double)); - overall_n = 0; - for ( i=0; i<fctx->nshells; i++ ) { - int j; - for ( j=0; j<fctx->n[i]; j++ ) { - overall_vec1[overall_n] = fctx->vec1[i][j]; - overall_vec2[overall_n] = fctx->vec2[i][j]; - overall_n++; - } - } - cc = gsl_stats_correlation(overall_vec1, 1, overall_vec2, 1, - overall_n); - free(overall_vec1); - free(overall_vec2); - break; - - case FOM_CRDANO : - overall_along_diagonal = malloc(fctx->nmax*sizeof(double)); - overall_perpend_diagonal = malloc(fctx->nmax*sizeof(double)); - overall_n = 0; - for ( i=0; i<fctx->nshells; i++ ) { - int j; - for ( j=0; j<fctx->n[i]; j++ ) { - overall_along_diagonal[overall_n] = - ( fctx->vec1[i][j] + fctx->vec2[i][j] ) - / sqrt(2.0); - overall_perpend_diagonal[overall_n] = - ( fctx->vec1[i][j] - fctx->vec2[i][j] ) - / sqrt(2.0); - overall_n++; - } - } - variance_signal = gsl_stats_variance_m(overall_along_diagonal, - 1, overall_n, 0.0); - variance_error = gsl_stats_variance_m(overall_perpend_diagonal, - 1, overall_n, 0.0); - cc = sqrt(variance_signal / variance_error ); - - free(overall_along_diagonal); - free(overall_perpend_diagonal); - break; - - case FOM_D1SIG : - case FOM_D2SIG : - overall_num = 0.0; - overall_den = 0.0; - for ( i=0; i<fctx->nshells; i++ ) { - overall_num += fctx->n_within[i]; - overall_den += fctx->cts[i]; - } - break; - - } - - switch ( fctx->fom ) { - - case FOM_R1I : - case FOM_R1F : - return overall_num/overall_den; - - case FOM_R2 : - return sqrt(overall_num/overall_den); - - case FOM_RSPLIT : - return 2.0*(overall_num/overall_den) / sqrt(2.0); - - case FOM_CC : - case FOM_CCANO : - case FOM_CRDANO : - return cc; - - case FOM_CCSTAR : - return sqrt((2.0*cc)/(1.0+cc)); - - case FOM_RANO : - return 2.0*(overall_num/overall_den); - - case FOM_RANORSPLIT : - return (2.0*(overall_num/overall_den)) / - (2.0*(overall_num2/overall_den2) / sqrt(2.0)); - - case FOM_D1SIG : - case FOM_D2SIG : - return overall_num/overall_den; - - } - - ERROR("This point is never reached.\n"); - abort(); -} - - -static double fom_shell(struct fom_context *fctx, int i) -{ - double cc; - int j; - double variance_signal; - double variance_error; - double *along_diagonal; - double *perpend_diagonal; - - switch ( fctx->fom ) { - - case FOM_R1I : - case FOM_R1F : - return fctx->num[i]/fctx->den[i]; - - case FOM_R2 : - return sqrt(fctx->num[i]/fctx->den[i]); - - case FOM_RSPLIT : - return 2.0*(fctx->num[i]/fctx->den[i]) / sqrt(2.0); - - case FOM_CC : - case FOM_CCANO : - return gsl_stats_correlation(fctx->vec1[i], 1, fctx->vec2[i], 1, - fctx->n[i]); - - case FOM_CCSTAR : - cc = gsl_stats_correlation(fctx->vec1[i], 1, fctx->vec2[i], 1, - fctx->n[i]); - return sqrt((2.0*cc)/(1.0+cc)); - - case FOM_RANO : - return 2.0 * fctx->num[i]/fctx->den[i]; - - case FOM_RANORSPLIT : - return (2.0*fctx->num[i]/fctx->den[i]) / - (2.0*(fctx->num2[i]/fctx->den2[i]) / sqrt(2.0)); - - case FOM_CRDANO : - along_diagonal = malloc(fctx->n[i] * sizeof(double)); - perpend_diagonal = malloc(fctx->n[i] * sizeof(double)); - for ( j=0; j<fctx->n[i]; j++ ) { - along_diagonal[j] = ( fctx->vec1[i][j] + - fctx->vec2[i][j] ) / sqrt(2.0); - perpend_diagonal[j] = ( fctx->vec1[i][j] - - fctx->vec2[i][j] ) / sqrt(2.0); - } - variance_signal = gsl_stats_variance_m(along_diagonal, 1, - fctx->n[i], 0.0); - variance_error = gsl_stats_variance_m(perpend_diagonal, 1, - fctx->n[i], 0.0); - free(along_diagonal); - free(perpend_diagonal); - return sqrt(variance_signal / variance_error); - - case FOM_D1SIG : - case FOM_D2SIG : - return (double)fctx->n_within[i] / fctx->cts[i]; - - } - - ERROR("This point is never reached.\n"); - abort(); -} - - -struct shells -{ - int config_intshells; - int nshells; - double *rmins; - double *rmaxs; -}; - - -static struct shells *set_intensity_shells(double min_I, double max_I, - int nshells) -{ - struct shells *s; - int i; - - if ( min_I >= max_I ) { - ERROR("Invalid intensity range.\n"); - return NULL; - } - - /* Adjust minimum and maximum intensities to get the most densely - * populated part of the reflections */ - max_I = min_I + (max_I-min_I)/5000.0; - - s = malloc(sizeof(struct shells)); - if ( s == NULL ) return NULL; - - s->rmins = malloc(nshells*sizeof(double)); - s->rmaxs = malloc(nshells*sizeof(double)); - - if ( (s->rmins==NULL) || (s->rmaxs==NULL) ) { - ERROR("Couldn't allocate memory for shells.\n"); - free(s); - return NULL; - } - - s->config_intshells = 1; - s->nshells = nshells; - - for ( i=0; i<nshells; i++ ) { - - s->rmins[i] = min_I + i*(max_I - min_I)/nshells;; - s->rmaxs[i] = min_I + (i+1)*(max_I - min_I)/nshells;; - - } - - return s; -} - - -static struct shells *set_resolution_shells(double rmin, double rmax, - int nshells) -{ - struct shells *s; - double total_vol, vol_per_shell; - int i; - - s = malloc(sizeof(struct shells)); - if ( s == NULL ) return NULL; - - s->rmins = malloc(nshells*sizeof(double)); - s->rmaxs = malloc(nshells*sizeof(double)); - - if ( (s->rmins==NULL) || (s->rmaxs==NULL) ) { - ERROR("Couldn't allocate memory for resolution shells.\n"); - free(s); - return NULL; - } - - s->config_intshells = 0; - s->nshells = nshells; - - total_vol = pow(rmax, 3.0) - pow(rmin, 3.0); - vol_per_shell = total_vol / nshells; - s->rmins[0] = rmin; - for ( i=1; i<nshells; i++ ) { - - double r; - - r = vol_per_shell + pow(s->rmins[i-1], 3.0); - r = pow(r, 1.0/3.0); - - /* Shells of constant volume */ - s->rmaxs[i-1] = r; - s->rmins[i] = r; - - } - s->rmaxs[nshells-1] = rmax; - - return s; -} - - -static double shell_label(struct shells *s, int i) -{ - if ( s->config_intshells ) { - return (i+0.5) / s->nshells; - } else { - return s->rmins[i] + (s->rmaxs[i] - s->rmins[i])/2.0; - } -} - - -static int get_bin(struct shells *s, Reflection *refl, UnitCell *cell) -{ - if ( s->config_intshells ) { - - double intensity; - int bin, j; - - intensity = get_intensity(refl); - - bin = -1; - for ( j=0; j<s->nshells; j++ ) { - if ( (intensity>s->rmins[j]) - && (intensity<=s->rmaxs[j]) ) - { - bin = j; - break; - } - } - - return bin; - - - } else { - - double d; - int bin, j; - signed int h, k, l; - - get_indices(refl, &h, &k, &l); - d = 2.0 * resolution(cell, h, k, l); - - bin = -1; - for ( j=0; j<s->nshells; j++ ) { - if ( (d>s->rmins[j]) && (d<=s->rmaxs[j]) ) { - bin = j; - break; - } - } - - /* Allow for slight rounding errors */ - if ( (bin == -1) && (d <= s->rmins[0]) ) bin = 0; - if ( (bin == -1) && (d >= s->rmaxs[s->nshells-1]) ) bin = 0; - assert(bin != -1); - - return bin; - - } -} - - -static int wilson_scale(RefList *list1, RefList *list2, UnitCell *cell) -{ - Reflection *refl1; - Reflection *refl2; - RefListIterator *iter; - int max_n = 256; - int n = 0; - double *x; - double *y; - int r; - double G, B; - double c0, c1, cov00, cov01, cov11, chisq; - - x = malloc(max_n*sizeof(double)); - y = malloc(max_n*sizeof(double)); - if ( (x==NULL) || (y==NULL) ) { - ERROR("Failed to allocate memory for scaling.\n"); - return 1; - } - - for ( refl1 = first_refl(list1, &iter); - refl1 != NULL; - refl1 = next_refl(refl1, iter) ) - { - signed int h, k, l; - double Ih1, Ih2; - double res; - - get_indices(refl1, &h, &k, &l); - res = resolution(cell, h, k, l); - - refl2 = find_refl(list2, h, k, l); - assert(refl2 != NULL); - - Ih1 = get_intensity(refl1); - Ih2 = get_intensity(refl2); - - if ( (Ih1 <= 0.0) || (Ih2 <= 0.0) ) continue; - if ( isnan(Ih1) || isinf(Ih1) ) continue; - if ( isnan(Ih2) || isinf(Ih2) ) continue; - - if ( n == max_n ) { - max_n *= 2; - x = realloc(x, max_n*sizeof(double)); - y = realloc(y, max_n*sizeof(double)); - if ( (x==NULL) || (y==NULL) ) { - ERROR("Failed to allocate memory for scaling.\n"); - return 1; - } - } - - x[n] = res*res; - y[n] = log(Ih1/Ih2); - n++; - - } - - if ( n < 2 ) { - ERROR("Not enough reflections for scaling\n"); - return 1; - } - - r = gsl_fit_linear(x, 1, y, 1, n, &c0, &c1, - &cov00, &cov01, &cov11, &chisq); - - if ( r ) { - ERROR("Scaling failed.\n"); - return 1; - } - - G = exp(c0); - B = c1/2.0; - - STATUS("Relative scale factor = %f, relative B factor = %f A^2\n", - G, B*1e20); - STATUS("A scale factor greater than 1 means that the second reflection " - "list is weaker than the first.\n"); - STATUS("A positive relative B factor means that the second reflection " - "list falls off with resolution more quickly than the first.\n"); - - free(x); - free(y); - - /* Apply the scaling factor */ - for ( refl2 = first_refl(list2, &iter); - refl2 != NULL; - refl2 = next_refl(refl2, iter) ) - { - signed int h, k, l; - double res; - double corr; - - get_indices(refl2, &h, &k, &l); - res = resolution(cell, h, k, l); - - corr = G * exp(2.0*B*res*res); - set_intensity(refl2, get_intensity(refl2)*corr); - set_esd_intensity(refl2, get_esd_intensity(refl2)*corr); - - } - return 0; -} - - static void do_fom(RefList *list1, RefList *list2, UnitCell *cell, double rmin, double rmax, enum fom fom, int config_unity, int nshells, const char *filename, int config_intshells, double min_I, double max_I, SymOpList *sym) { - int i; - Reflection *refl1; - RefListIterator *iter; - FILE *fh; - struct fom_context *fctx; struct shells *shells; + struct fom_context *fctx; + FILE *fh; + int i; const char *t1, *t2; - int n_out; - - fctx = init_fom(fom, num_reflections(list1), nshells); - - if ( fctx==NULL ) { - ERROR("Couldn't allocate memory for resolution shells.\n"); - return; - } - - if ( !config_unity && wilson_scale(list1, list2, cell) ) { - ERROR("Error with scaling.\n"); - return; - } /* Calculate the bins */ if ( config_intshells ) { - shells = set_intensity_shells(min_I, max_I, nshells); + shells = make_intensity_shells(min_I, max_I, nshells); } else { - shells = set_resolution_shells(rmin, rmax, nshells); + shells = make_resolution_shells(rmin, rmax, nshells); } if ( shells == NULL ) { @@ -806,95 +106,8 @@ static void do_fom(RefList *list1, RefList *list2, UnitCell *cell, return; } - 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; - - get_indices(refl1, &h, &k, &l); - - refl2 = find_refl(list2, h, k, l); - if ( refl2 == NULL ) continue; - - bin = get_bin(shells, refl1, cell); - if ( bin == -1 ) { - n_out++; - continue; - } - - i1 = get_intensity(refl1); - i2 = get_intensity(refl2); - sig1 = get_esd_intensity(refl1); - sig2 = get_esd_intensity(refl2); - - 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, - &hb, &kb, &lb) ) - { - refl1_bij = find_refl(list1, hb, kb, lb); - } - - if ( find_equiv_in_list(list2, -h, -k, -l, sym, - &hb, &kb, &lb) ) - { - refl2_bij = find_refl(list2, hb, kb, lb); - } - - /* Each reflection must only be counted once, whether - * we are visiting it now as "normal" or "bij" */ - if ( get_flag(refl1) ) continue; - assert(!get_flag(refl2)); - set_flag(refl1, 1); - set_flag(refl1_bij, 1); - set_flag(refl2, 1); - set_flag(refl2_bij, 1); - - 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); - } - if ( n_out) { - ERROR("WARNING: %i reflection pairs outside range.\n", n_out); - } + fctx = fom_calculate(list1, list2, cell, shells, fom, + config_unity, sym); switch ( fom ) { @@ -1026,69 +239,69 @@ static void do_fom(RefList *list1, RefList *list2, UnitCell *cell, switch ( fom ) { - case FOM_R1I : - case FOM_R1F : - case FOM_R2 : - case FOM_RSPLIT : - case FOM_RANO : - if ( config_intshells ) { - fprintf(fh, "%10.3f %10.2f %10i\n", - cen, r*100.0, fctx->cts[i]); - } else { - fprintf(fh, "%10.3f %10.2f %10i %10.2f " - "%10.3f %10.3f\n", - cen*1.0e-9, r*100.0, fctx->cts[i], - (1.0/cen)*1e10, - shells->rmins[i]*1.0e-9, - shells->rmaxs[i]*1.0e-9); - } - break; + case FOM_R1I : + case FOM_R1F : + case FOM_R2 : + case FOM_RSPLIT : + case FOM_RANO : + if ( config_intshells ) { + fprintf(fh, "%10.3f %10.2f %10i\n", + cen, r*100.0, fctx->cts[i]); + } else { + fprintf(fh, "%10.3f %10.2f %10i %10.2f " + "%10.3f %10.3f\n", + cen*1.0e-9, r*100.0, fctx->cts[i], + (1.0/cen)*1e10, + shells->rmins[i]*1.0e-9, + shells->rmaxs[i]*1.0e-9); + } + break; - case FOM_CC : - case FOM_CCSTAR : - case FOM_CCANO : - case FOM_CRDANO : - if ( config_intshells ) { - fprintf(fh, "%10.3f %10.7f %10i\n", - cen, r, fctx->cts[i]); - } else { - 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, - shells->rmins[i]*1.0e-9, - shells->rmaxs[i]*1.0e-9); - } - break; + case FOM_CC : + case FOM_CCSTAR : + case FOM_CCANO : + case FOM_CRDANO : + if ( config_intshells ) { + fprintf(fh, "%10.3f %10.7f %10i\n", + cen, r, fctx->cts[i]); + } else { + 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, + shells->rmins[i]*1.0e-9, + shells->rmaxs[i]*1.0e-9); + } + break; - case FOM_RANORSPLIT : - if ( config_intshells ) { - fprintf(fh, "%10.3f %10.7f %10i\n", - cen, r, fctx->cts[i]); - } else { - 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, - shells->rmins[i]*1.0e-9, - shells->rmaxs[i]*1.0e-9); - } - break; + case FOM_RANORSPLIT : + if ( config_intshells ) { + fprintf(fh, "%10.3f %10.7f %10i\n", + cen, r, fctx->cts[i]); + } else { + 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, + shells->rmins[i]*1.0e-9, + shells->rmaxs[i]*1.0e-9); + } + break; - case FOM_D1SIG : - case FOM_D2SIG : - if ( config_intshells ) { - fprintf(fh, "%10.3f %10.2f %10i\n", - cen, r*100.0, fctx->cts[i]); - } else { - fprintf(fh, "%10.3f %10.2f %10i %10.2f " - "%10.3f %10.3f\n", - cen*1.0e-9, r*100.0, fctx->cts[i], - (1.0/cen)*1e10, - shells->rmins[i]*1.0e-9, - shells->rmaxs[i]*1.0e-9); - } - break; + case FOM_D1SIG : + case FOM_D2SIG : + if ( config_intshells ) { + fprintf(fh, "%10.3f %10.2f %10i\n", + cen, r*100.0, fctx->cts[i]); + } else { + fprintf(fh, "%10.3f %10.2f %10i %10.2f " + "%10.3f %10.3f\n", + cen*1.0e-9, r*100.0, fctx->cts[i], + (1.0/cen)*1e10, + shells->rmins[i]*1.0e-9, + shells->rmaxs[i]*1.0e-9); + } + break; - } + } } |