From 491295f2d71efae499ee0b4b87275260b971e25e Mon Sep 17 00:00:00 2001 From: Lorenzo Galli Date: Fri, 21 Jun 2013 16:58:39 +0200 Subject: Added rms correlation ratio to compare_hkl (CRDano) --- src/compare_hkl.c | 86 +++++++++++++++++++++++++++++++++++++++++++++++++------ 1 file changed, 77 insertions(+), 9 deletions(-) (limited to 'src') diff --git a/src/compare_hkl.c b/src/compare_hkl.c index 687d3a31..47b52419 100644 --- a/src/compare_hkl.c +++ b/src/compare_hkl.c @@ -8,6 +8,7 @@ * * Authors: * 2010-2012 Thomas White + * 2013 Lorenzo Galli * * This file is part of CrystFEL. * @@ -56,7 +57,8 @@ enum fom FOM_RSPLIT, FOM_CC, FOM_CCSTAR, - FOM_CCANO + FOM_CCANO, + FOM_CRDANO }; @@ -69,6 +71,7 @@ static enum fom get_fom(const char *s) 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; ERROR("Unknown figure of merit '%s'.\n", s); exit(1); @@ -84,7 +87,8 @@ static void show_help(const char *s) " -y, --symmetry= The symmetry of both the input files.\n" " -p, --pdb= PDB file to use.\n" " --fom= Calculate this figure of merit Choose from:\n" -" R1I, R1F, R2, Rsplit, CC, CCstar, CCano.\n" +" R1I, R1F, R2, Rsplit, CC, CCstar,\n" +" CCano, CRDano.\n" " --nshells= Use resolution shells.\n" " -u Force scale factor to 1.\n" " --shell-file= Write resolution shells to .\n" @@ -149,6 +153,7 @@ static struct fom_context *init_fom(enum fom fom, int nmax, int nshells) 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; @@ -213,6 +218,7 @@ static void add_to_fom(struct fom_context *fctx, double i1, double i2, 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; @@ -232,6 +238,10 @@ static double fom_overall(struct fom_context *fctx) 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 ) { @@ -268,6 +278,32 @@ static double fom_overall(struct fom_context *fctx) 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; inshells; i++ ) { + int j; + for ( j=0; jn[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(overall_along_diagonal, 1, + overall_n); + variance_error = gsl_stats_variance(overall_perpend_diagonal, 1, + overall_n); + cc = sqrt(variance_signal / variance_error ); + + free(overall_along_diagonal); + free(overall_perpend_diagonal); + break; + } switch ( fctx->fom ) { @@ -284,6 +320,7 @@ static double fom_overall(struct fom_context *fctx) case FOM_CC : case FOM_CCANO : + case FOM_CRDANO : return cc; case FOM_CCSTAR : @@ -299,6 +336,11 @@ static double fom_overall(struct fom_context *fctx) static double fom_shell(struct fom_context *fctx, int i) { double cc; + int j; + double variance_signal; + double variance_error; + double along_diagonal[fctx->n[i]]; + double perpend_diagonal[fctx->n[i]]; switch ( fctx->fom ) { @@ -322,6 +364,20 @@ static double fom_shell(struct fom_context *fctx, int i) fctx->n[i]); return sqrt((2.0*cc)/(1.0+cc)); + + case FOM_CRDANO : + for ( j=0; jn[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(along_diagonal, 1, + fctx->n[i]); + variance_error = gsl_stats_variance(perpend_diagonal, 1, + fctx->n[i]); + return sqrt(variance_signal / variance_error); + } ERROR("This point is never reached.\n"); @@ -558,7 +614,7 @@ static void do_fom(RefList *list1, RefList *list2, UnitCell *cell, i1 = get_intensity(refl1); i2 = scale * get_intensity(refl2); - if ( fom == FOM_CCANO ) { + if ( fom == FOM_CCANO || fom == FOM_CRDANO ) { Reflection *refl1_bij = NULL; Reflection *refl2_bij = NULL; @@ -627,6 +683,10 @@ static void do_fom(RefList *list1, RefList *list2, UnitCell *cell, STATUS("Overall CCano = %.7f\n", fom_overall(fctx)); break; + case FOM_CRDANO : + STATUS("Overall CRDano = %.7f\n", fom_overall(fctx)); + break; + } fh = fopen(filename, "w"); @@ -673,6 +733,10 @@ static void do_fom(RefList *list1, RefList *list2, UnitCell *cell, fprintf(fh, "%s CCano nref%s\n", t1, t2); break; + case FOM_CRDANO : + fprintf(fh, "%s CRDano nref%s\n", t1, t2); + break; + } for ( i=0; i 0 ) { STATUS("%i reflection pairs rejected because either or both" - " versions did not have Bijvoet partners.\n", nbij); + " versions did not have Bijvoet partners.\n", nres); } if ( ncen > 0 ) { -- cgit v1.2.3