Add "-u" (unity scaling) option to compare_hkl

3 files changed, 88 insertions, 69 deletions

diff --git a/src/compare_hkl.c b/src/compare_hkl.c
index 4ab2dc70..02a85d2f 100644
--- a/src/compare_hkl.c
+++ b/src/compare_hkl.c
@@ -277,6 +277,7 @@ int main(int argc, char *argv[])
 	Reflection *refl1;
 	RefListIterator *iter;
 	double *arr2;
+	int config_unity = 0;
 
 	/* Long options */
 	const struct option longopts[] = {
@@ -291,7 +292,9 @@ int main(int argc, char *argv[])
 	};
 
 	/* Short options */
-	while ((c = getopt_long(argc, argv, "ho:y:p:", longopts, NULL)) != -1) {
+	while ((c = getopt_long(argc, argv, "ho:y:p:u",
+	                        longopts, NULL)) != -1)
+	{
 
 		switch (c) {
 		case 'h' :
@@ -310,6 +313,10 @@ int main(int argc, char *argv[])
 			pdb = strdup(optarg);
 			break;
 
+		case 'u' :
+			config_unity = 1;
+			break;
+
 		case 0 :
 			break;
 
@@ -473,33 +480,33 @@ int main(int argc, char *argv[])
 	arr2 = intensities_from_list(list2_transformed);
 	reflist_free(list2_transformed);
 
-	R1 = stat_r1_ignore(list1, arr2, &scale_r1fi);
+	R1 = stat_r1_ignore(list1, arr2, &scale_r1fi, config_unity);
 	STATUS("R1(F) = %5.4f %% (scale=%5.2e)"
 	       " (ignoring negative intensities)\n",
 	       R1*100.0, scale_r1fi);
 
-	R1 = stat_r1_zero(list1, arr2, &scale_r1);
+	R1 = stat_r1_zero(list1, arr2, &scale_r1, config_unity);
 	STATUS("R1(F) = %5.4f %% (scale=%5.2e)"
 	       " (zeroing negative intensities)\n",
 	       R1*100.0, scale_r1);
 
-	R2 = stat_r2(list1, arr2, &scale_r2);
+	R2 = stat_r2(list1, arr2, &scale_r2, config_unity);
 	STATUS("R2(I) = %5.4f %% (scale=%5.2e)\n", R2*100.0, scale_r2);
 
-	R1i = stat_r1_i(list1, arr2, &scale_r1i);
+	R1i = stat_r1_i(list1, arr2, &scale_r1i, config_unity);
 	STATUS("R1(I) = %5.4f %% (scale=%5.2e)\n", R1i*100.0, scale_r1i);
 
-	Rdiff = stat_rdiff_ignore(list1, arr2, &scale_rdig);
+	Rdiff = stat_rdiff_ignore(list1, arr2, &scale_rdig, config_unity);
 	STATUS("Rint(F) = %5.4f %% (scale=%5.2e)"
 	       " (ignoring negative intensities)\n",
 	       Rdiff*100.0, scale_rdig);
 
-	Rdiff = stat_rdiff_zero(list1, arr2, &scale);
+	Rdiff = stat_rdiff_zero(list1, arr2, &scale, config_unity);
 	STATUS("Rint(F) = %5.4f %% (scale=%5.2e)"
 	       " (zeroing negative intensities)\n",
 	       Rdiff*100.0, scale);
 
-	Rdiff = stat_rdiff_intensity(list1, arr2, &scale_rintint);
+	Rdiff = stat_rdiff_intensity(list1, arr2, &scale_rintint, config_unity);
 	STATUS("Rint(I) = %5.4f %% (scale=%5.2e)\n",
 	       Rdiff*100.0, scale_rintint);
 
diff --git a/src/statistics.c b/src/statistics.c
index 2312c4a7..791952b8 100644
--- a/src/statistics.c
+++ b/src/statistics.c
@@ -376,7 +376,8 @@ static double calc_r(double scale, void *params)
 }
 
 
-static double r_minimised(RefList *list1, double *arr2, double *scalep, int fom)
+static double r_minimised(RefList *list1, double *arr2, double *scalep, int fom,
+                          int u)
 {
 	gsl_function F;
 	gsl_min_fminimizer *s;
@@ -389,57 +390,65 @@ static double r_minimised(RefList *list1, double *arr2, double *scalep, int fom)
 	rp.arr2 = arr2;
 	rp.fom = fom;
 
-	F.function = &calc_r;
-	F.params = &rp;
+	if ( u ) {
 
-	s = gsl_min_fminimizer_alloc(gsl_min_fminimizer_brent);
+		scale = 1.0;
 
-	/* Initial guess */
-	switch ( fom ) {
-	case R_1_ZERO :
-	case R_1_IGNORE :
-	case R_DIFF_ZERO :
-	case R_DIFF_IGNORE :
-		scale = stat_scale_sqrti(list1, arr2);
-		break;
-	case R_2 :
-	case R_1_I :
-	case R_DIFF_INTENSITY :
-		scale = stat_scale_intensity(list1, arr2);
-		break;
-	}
-	//STATUS("Initial scale factor estimate: %5.2e\n", scale);
-
-	/* Probably within an order of magnitude either side */
-	gsl_min_fminimizer_set(s, &F, scale, scale/10.0, scale*10.0);
+	} else {
 
-	do {
+		F.function = &calc_r;
+		F.params = &rp;
 
-		double lo, up;
+		s = gsl_min_fminimizer_alloc(gsl_min_fminimizer_brent);
 
-		/* Iterate */
-		if ( gsl_min_fminimizer_iterate(s) ) {
-			ERROR("Failed to find scale factor.\n");
-			return NAN;
+		/* Initial guess */
+		switch ( fom ) {
+		case R_1_ZERO :
+		case R_1_IGNORE :
+		case R_DIFF_ZERO :
+		case R_DIFF_IGNORE :
+			scale = stat_scale_sqrti(list1, arr2);
+			break;
+		case R_2 :
+		case R_1_I :
+		case R_DIFF_INTENSITY :
+			scale = stat_scale_intensity(list1, arr2);
+			break;
 		}
+		//STATUS("Initial scale factor estimate: %5.2e\n", scale);
 
-		/* Get the current estimate */
-		scale = gsl_min_fminimizer_x_minimum(s);
-		lo = gsl_min_fminimizer_x_lower(s);
-		up = gsl_min_fminimizer_x_upper(s);
+		/* Probably within an order of magnitude either side */
+		gsl_min_fminimizer_set(s, &F, scale, scale/10.0, scale*10.0);
 
-		/* Check for convergence */
-		status = gsl_min_test_interval(lo, up, 0.001, 0.0);
+		do {
 
-		iter++;
+			double lo, up;
 
-	} while ( status == GSL_CONTINUE );
+			/* Iterate */
+			if ( gsl_min_fminimizer_iterate(s) ) {
+				ERROR("Failed to find scale factor.\n");
+				return NAN;
+			}
 
-	if (status != GSL_SUCCESS) {
-		ERROR("Scale factor minimisation failed.\n");
-	}
+			/* Get the current estimate */
+			scale = gsl_min_fminimizer_x_minimum(s);
+			lo = gsl_min_fminimizer_x_lower(s);
+			up = gsl_min_fminimizer_x_upper(s);
+
+			/* Check for convergence */
+			status = gsl_min_test_interval(lo, up, 0.001, 0.0);
+
+			iter++;
+
+		} while ( status == GSL_CONTINUE );
 
-	gsl_min_fminimizer_free(s);
+		if ( status != GSL_SUCCESS ) {
+			ERROR("Scale factor minimisation failed.\n");
+		}
+
+		gsl_min_fminimizer_free(s);
+
+	}
 
 	//STATUS("Final scale factor: %5.2e\n", scale);
 	*scalep = scale;
@@ -447,45 +456,45 @@ static double r_minimised(RefList *list1, double *arr2, double *scalep, int fom)
 }
 
 
-double stat_r1_ignore(RefList *list1, double *arr2, double *scalep)
+double stat_r1_ignore(RefList *list1, double *arr2, double *scalep, int u)
 {
-	return r_minimised(list1, arr2, scalep, R_1_IGNORE);
+	return r_minimised(list1, arr2, scalep, R_1_IGNORE, u);
 }
 
 
-double stat_r1_zero(RefList *list1, double *arr2, double *scalep)
+double stat_r1_zero(RefList *list1, double *arr2, double *scalep, int u)
 {
-	return r_minimised(list1, arr2, scalep, R_1_ZERO);
+	return r_minimised(list1, arr2, scalep, R_1_ZERO, u);
 }
 
 
-double stat_r2(RefList *list1, double *arr2, double *scalep)
+double stat_r2(RefList *list1, double *arr2, double *scalep, int u)
 {
-	return r_minimised(list1, arr2, scalep, R_2);
+	return r_minimised(list1, arr2, scalep, R_2, u);
 }
 
 
-double stat_r1_i(RefList *list1, double *arr2, double *scalep)
+double stat_r1_i(RefList *list1, double *arr2, double *scalep, int u)
 {
-	return r_minimised(list1, arr2, scalep, R_1_I);
+	return r_minimised(list1, arr2, scalep, R_1_I, u);
 }
 
 
-double stat_rdiff_zero(RefList *list1, double *arr2, double *scalep)
+double stat_rdiff_zero(RefList *list1, double *arr2, double *scalep, int u)
 {
-	return r_minimised(list1, arr2, scalep, R_DIFF_ZERO);
+	return r_minimised(list1, arr2, scalep, R_DIFF_ZERO, u);
 }
 
 
-double stat_rdiff_ignore(RefList *list1, double *arr2, double *scalep)
+double stat_rdiff_ignore(RefList *list1, double *arr2, double *scalep, int u)
 {
-	return r_minimised(list1, arr2, scalep, R_DIFF_IGNORE);
+	return r_minimised(list1, arr2, scalep, R_DIFF_IGNORE, u);
 }
 
 
-double stat_rdiff_intensity(RefList *list1, double *arr2, double *scalep)
+double stat_rdiff_intensity(RefList *list1, double *arr2, double *scalep, int u)
 {
-	return r_minimised(list1, arr2,  scalep, R_DIFF_INTENSITY);
+	return r_minimised(list1, arr2,  scalep, R_DIFF_INTENSITY, u);
 }
 
 
diff --git a/src/statistics.h b/src/statistics.h
index e1d7244e..ced295de 100644
--- a/src/statistics.h
+++ b/src/statistics.h
@@ -21,17 +21,20 @@
 
 extern double stat_scale_intensity(RefList *list1, double *arr2);
 
-extern double stat_r1_zero(RefList *list1, double *arr2, double *scalep);
-extern double stat_r1_ignore(RefList *list1, double *arr2, double *scalep);
+extern double stat_r1_zero(RefList *list1, double *arr2, double *scalep, int u);
+extern double stat_r1_ignore(RefList *list1, double *arr2,
+                             double *scalep, int u);
 
-extern double stat_r2(RefList *list1, double *arr2, double *scalep);
+extern double stat_r2(RefList *list1, double *arr2, double *scalep, int u);
 
-extern double stat_r1_i(RefList *list1, double *arr2, double *scalep);
+extern double stat_r1_i(RefList *list1, double *arr2, double *scalep, int u);
 
-extern double stat_rdiff_zero(RefList *list1, double *arr2, double *scalep);
-extern double stat_rdiff_ignore(RefList *list1, double *arr2, double *scalep);
+extern double stat_rdiff_zero(RefList *list1, double *arr2,
+                              double *scalep, int u);
+extern double stat_rdiff_ignore(RefList *list1, double *arr2,
+                                double *scalep, int u);
 extern double stat_rdiff_intensity(RefList *list1, double *arr2,
-                                   double *scalep);
+                                   double *scalep, int u);
 
 extern double stat_pearson_i(RefList *list1, double *arr2);
 extern double stat_pearson_f_zero(RefList *list1, double *arr2);