compare_hkl: Add d1sig and d2sig

2 files changed, 115 insertions, 24 deletions

diff --git a/doc/man/compare_hkl.1 b/doc/man/compare_hkl.1
index 270e990c..b628f0b6 100644
--- a/doc/man/compare_hkl.1
+++ b/doc/man/compare_hkl.1
@@ -70,6 +70,9 @@ Note that this figure of merit compares measurements within one data set, so it
 .IP \fBRano/Rsplit\fR
 .PD
 The ratio of Rano to Rsplit, as defined above.
+.IP "\fBd1sig\fR and \fBd2sig\fR"
+.PD
+The fraction of differences between intensities which are within one times (for \fBd1sig\fR) and two times (for \fBd2sig\fR) the mean of the corresponding sigma(I) values.
 .PP
 I1 and I2 are the intensities of the same reflection in both reflection lists.  The scale factor, k, is given by sum(I1*i2) / sum(I2^2), unless you use \fB-u\fR.
 .RE
diff --git a/src/compare_hkl.c b/src/compare_hkl.c
index d3f60feb..eed92af0 100644
--- a/src/compare_hkl.c
+++ b/src/compare_hkl.c
@@ -60,7 +60,9 @@ enum fom
 	FOM_CCANO,
 	FOM_CRDANO,
 	FOM_RANO,
-	FOM_RANORSPLIT
+	FOM_RANORSPLIT,
+	FOM_D1SIG,
+	FOM_D2SIG
 };
 
 
@@ -76,6 +78,8 @@ static enum fom get_fom(const char *s)
 	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);
@@ -92,7 +96,8 @@ static void show_help(const char *s)
 "  -p, --pdb=<filename>       Unit cell file to use.\n"
 "      --fom=<FoM>            Calculate this figure of merit  Choose from:\n"
 "                              R1I, R1F, R2, Rsplit, CC, CCstar,\n"
-"			       CCano, CRDano, Rano, Rano/Rsplit\n"
+"			       CCano, CRDano, Rano, Rano/Rsplit, d1sig,\n"
+"                              d2sig\n"
 "      --nshells=<n>          Use <n> resolution shells.\n"
 "  -u                         Force scale factor to 1.\n"
 "      --shell-file=<file>    Write resolution shells to <file>.\n"
@@ -118,6 +123,7 @@ struct fom_context
 {
 	enum fom fom;
 	int nshells;
+	int *cts;
 
 	/* For R-factors */
 	double *num;
@@ -132,6 +138,9 @@ struct fom_context
 	double **vec2;
 	int *n;
 	int nmax;
+
+	/* For "counting" things e.g. d1sig or d2sig */
+	int *n_within;
 };
 
 
@@ -145,6 +154,10 @@ static struct fom_context *init_fom(enum fom fom, int nmax, int nshells)
 
 	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 ) {
 
@@ -190,10 +203,18 @@ static struct fom_context *init_fom(enum fom fom, int nmax, int nshells)
 		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;
@@ -201,11 +222,14 @@ static struct fom_context *init_fom(enum fom fom, int nmax, int nshells)
 
 
 static void add_to_fom(struct fom_context *fctx, double i1, double i2,
-                       double i1bij, double i2bij, int bin)
+                       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);
@@ -260,6 +284,17 @@ static void add_to_fom(struct fom_context *fctx, double i1, double i2,
 		fctx->den[bin] += im + imbij;
 		break;
 
+		case FOM_D1SIG :
+		if ( fabs(i1-i2) < (sig1+sig2)/2.0 ) {
+			fctx->n_within[bin]++;
+		}
+		break;
+
+		case FOM_D2SIG :
+		if ( fabs(i1-i2) < sig1+sig2 ) {  /* = 2 * (sig1+sig2)/2 */
+			fctx->n_within[bin]++;
+		}
+		break;
 
 	}
 }
@@ -357,6 +392,16 @@ static double fom_overall(struct fom_context *fctx)
 		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 ) {
@@ -386,6 +431,10 @@ static double fom_overall(struct fom_context *fctx)
 		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");
@@ -448,6 +497,10 @@ static double fom_shell(struct fom_context *fctx, int i)
 		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");
@@ -614,7 +667,6 @@ static void do_fom(RefList *list1, RefList *list2, UnitCell *cell,
                    int config_intshells, double min_I, double max_I,
                    SymOpList *sym)
 {
-	int *cts;
 	int i;
 	Reflection *refl1;
 	RefListIterator *iter;
@@ -625,18 +677,13 @@ static void do_fom(RefList *list1, RefList *list2, UnitCell *cell,
 	const char *t1, *t2;
 	int n_out;
 
-	cts = malloc(nshells*sizeof(int));
 	fctx = init_fom(fom, num_reflections(list1), nshells);
 
-	if ( (fctx==NULL) || (cts==NULL) ) {
+	if ( fctx==NULL ) {
 		ERROR("Couldn't allocate memory for resolution shells.\n");
 		return;
 	}
 
-	for ( i=0; i<nshells; i++ ) {
-		cts[i] = 0;
-	}
-
 	if ( config_unity ) {
 		scale = 1.0;
 	} else {
@@ -664,6 +711,7 @@ static void do_fom(RefList *list1, RefList *list2, UnitCell *cell,
 		int bin;
 		double i1, i2;
 		double i1bij, i2bij;
+		double sig1, sig2;
 		Reflection *refl2;
 
 		get_indices(refl1, &h, &k, &l);
@@ -679,6 +727,8 @@ static void do_fom(RefList *list1, RefList *list2, UnitCell *cell,
 
 		i1 = get_intensity(refl1);
 		i2 = scale * get_intensity(refl2);
+		sig1 = get_esd_intensity(refl1);
+		sig2 = scale * get_esd_intensity(refl2);
 
 		if ( (fom == FOM_CCANO) || (fom == FOM_CRDANO)
 		  || (fom == FOM_RANO) || (fom == FOM_RANORSPLIT) )
@@ -714,8 +764,7 @@ static void do_fom(RefList *list1, RefList *list2, UnitCell *cell,
 
 		}
 
-		add_to_fom(fctx, i1, i2, i1bij, i2bij, bin);
-		cts[bin]++;
+		add_to_fom(fctx, i1, i2, i1bij, i2bij, sig1, sig2, bin);
 	}
 	if ( n_out)  {
 		ERROR("WARNING: %i reflection pairs outside range.\n", n_out);
@@ -763,6 +812,17 @@ static void do_fom(RefList *list1, RefList *list2, UnitCell *cell,
 		STATUS("Overall Rano/Rsplit =  %.7f\n", fom_overall(fctx));
 		break;
 
+		case FOM_D1SIG :
+		STATUS("Fraction of differences less than 1 sigma = %.7f %%\n",
+		       100.0*fom_overall(fctx));
+		break;
+
+		case FOM_D2SIG :
+		STATUS("Fraction of differences less than 2 sigma = %.7f %%\n",
+		       100.0*fom_overall(fctx));
+		break;
+
+
 	}
 
 	fh = fopen(filename, "w");
@@ -821,6 +881,14 @@ static void do_fom(RefList *list1, RefList *list2, UnitCell *cell,
 		fprintf(fh, "%s Rano/Rsplit       nref%s\n", t1, t2);
 		break;
 
+		case FOM_D1SIG :
+		fprintf(fh, "%s D<1sigma/%%     nref%s\n", t1, t2);
+		break;
+
+		case FOM_D2SIG :
+		fprintf(fh, "%s D<2sigma/%%     nref%s\n", t1, t2);
+		break;
+
 	}
 
 	for ( i=0; i<nshells; i++ ) {
@@ -839,13 +907,14 @@ static void do_fom(RefList *list1, RefList *list2, UnitCell *cell,
 		case FOM_RANO :
 		if ( config_intshells ) {
 			fprintf(fh, "%10.3f %10.2f %10i\n",
-				cen, r*100.0, cts[i]);
+				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, cts[i], (1.0/cen)*1e10,
-				shells->rmins[i]*1.0e-9,
-				shells->rmaxs[i]*1.0e-9);
+			        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;
 
@@ -855,28 +924,43 @@ static void do_fom(RefList *list1, RefList *list2, UnitCell *cell,
 		case FOM_CRDANO :
 		if ( config_intshells ) {
 			fprintf(fh, "%10.3f %10.7f %10i\n",
-				cen, r, cts[i]);
+				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, cts[i], (1.0/cen)*1e10,
-				shells->rmins[i]*1.0e-9,
-				shells->rmaxs[i]*1.0e-9);
+			        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, cts[i]);
+				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, cts[i], (1.0/cen)*1e10,
-				shells->rmins[i]*1.0e-9, shells->rmaxs[i]*1.0e-9);
+			        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;
 
 	}
 
@@ -1085,6 +1169,8 @@ int main(int argc, char *argv[])
 			case FOM_CRDANO :
 			case FOM_RANO :
 			case FOM_RANORSPLIT :
+			case FOM_D1SIG :
+			case FOM_D2SIG :
 			break;
 		}
 	}
@@ -1104,6 +1190,8 @@ int main(int argc, char *argv[])
 			case FOM_RSPLIT :
 			case FOM_CC :
 			case FOM_CCSTAR :
+			case FOM_D1SIG :
+			case FOM_D2SIG :
 			break;
 
 			case FOM_CCANO :