From d055c5c8f2e489bae2d34d567d6e2457910b4c51 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Thu, 21 Nov 2013 22:03:24 -0800 Subject: compare_hkl: Add Rano and Rano/Rsplit --- src/compare_hkl.c | 153 ++++++++++++++++++++++++++++++++++++++++++++---------- 1 file changed, 127 insertions(+), 26 deletions(-) (limited to 'src/compare_hkl.c') diff --git a/src/compare_hkl.c b/src/compare_hkl.c index 66b58adc..88e162b5 100644 --- a/src/compare_hkl.c +++ b/src/compare_hkl.c @@ -3,11 +3,11 @@ * * Compare reflection lists * - * Copyright © 2012 Deutsches Elektronen-Synchrotron DESY, - * a research centre of the Helmholtz Association. + * Copyright © 2012-2013 Deutsches Elektronen-Synchrotron DESY, + * a research centre of the Helmholtz Association. * * Authors: - * 2010-2012 Thomas White + * 2010-2013 Thomas White * 2013 Lorenzo Galli * * This file is part of CrystFEL. @@ -58,7 +58,9 @@ enum fom FOM_CC, FOM_CCSTAR, FOM_CCANO, - FOM_CRDANO + FOM_CRDANO, + FOM_RANO, + FOM_RANORSPLIT }; @@ -72,6 +74,8 @@ static enum fom get_fom(const char *s) 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; ERROR("Unknown figure of merit '%s'.\n", s); exit(1); @@ -88,7 +92,7 @@ static void show_help(const char *s) " -p, --pdb= PDB file to use.\n" " --fom= Calculate this figure of merit Choose from:\n" " R1I, R1F, R2, Rsplit, CC, CCstar,\n" -" CCano, CRDano.\n" +" CCano, CRDano, Rano, Rano/Rsplit\n" " --nshells= Use resolution shells.\n" " -u Force scale factor to 1.\n" " --shell-file= Write resolution shells to .\n" @@ -116,6 +120,10 @@ struct fom_context double *num; double *den; + /* For "double" R-factors */ + double *num2; + double *den2; + /* For CCs */ double **vec1; double **vec2; @@ -137,10 +145,21 @@ static struct fom_context *init_fom(enum fom fom, int nmax, int nshells) 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; inum2[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; @@ -182,6 +201,7 @@ static void add_to_fom(struct fom_context *fctx, double i1, double i2, double i1bij, double i2bij, int bin) { double f1, f2; + double im, imbij; /* Negative intensities have already been weeded out. */ f1 = sqrt(i1); @@ -225,6 +245,18 @@ static void add_to_fom(struct fom_context *fctx, double i1, double i2, 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; + } } @@ -234,6 +266,8 @@ 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; @@ -250,12 +284,28 @@ static double fom_overall(struct fom_context *fctx) case FOM_R1F : case FOM_R2 : case FOM_RSPLIT : + case FOM_RANO : + overall_num = 0.0; + overall_den = 0.0; + for ( i=0; inshells; 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; inshells; i++ ) { overall_num += fctx->num[i]; overall_den += fctx->den[i]; } + overall_num2 = 0.0; + overall_den2 = 0.0; + for ( i=0; inshells; i++ ) { + overall_num2 += fctx->num2[i]; + overall_den2 += fctx->den2[i]; + } break; case FOM_CC : @@ -326,6 +376,13 @@ static double fom_overall(struct fom_context *fctx) 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)); + } ERROR("This point is never reached.\n"); @@ -364,6 +421,12 @@ static double fom_shell(struct fom_context *fctx, int i) 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)); @@ -618,7 +681,9 @@ static void do_fom(RefList *list1, RefList *list2, UnitCell *cell, i1 = get_intensity(refl1); i2 = scale * get_intensity(refl2); - if ( fom == FOM_CCANO || fom == FOM_CRDANO ) { + if ( (fom == FOM_CCANO) || (fom == FOM_CRDANO) + || (fom == FOM_RANO) || (fom == FOM_RANORSPLIT) ) + { Reflection *refl1_bij = NULL; Reflection *refl2_bij = NULL; @@ -691,6 +756,14 @@ static void do_fom(RefList *list1, RefList *list2, UnitCell *cell, STATUS("Overall CRDano = %.7f\n", fom_overall(fctx)); break; + case FOM_RANO : + STATUS("Overall Rano = %.2f %%\n", 100.0*fom_overall(fctx)); + break; + + case FOM_RANORSPLIT : + STATUS("Overall Rano/Rsplit = %.7f\n", fom_overall(fctx)); + break; + } fh = fopen(filename, "w"); @@ -741,6 +814,14 @@ static void do_fom(RefList *list1, RefList *list2, UnitCell *cell, fprintf(fh, "%s CRDano nref%s\n", t1, t2); break; + case FOM_RANO : + fprintf(fh, "%s Rano/%% nref%s\n", t1, t2); + break; + + case FOM_RANORSPLIT : + fprintf(fh, "%s Rano/Rsplit nref%s\n", t1, t2); + break; + } for ( i=0; i