aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorThomas White <taw@physics.org>2010-07-07 16:51:26 +0200
committerThomas White <taw@physics.org>2012-02-22 15:26:53 +0100
commitb3854b29262d2f67eae90c1e54b39b09e49ea66b (patch)
tree2e01dad8daeba92789ee0c9a9bbdd1c22edc21f3
parentc30c05aeeadd8caf0cd887fab0024bf894ae7d65 (diff)
compare_hkl: Use minimisation to determine the right scale factor
-rw-r--r--src/compare_hkl.c4
-rw-r--r--src/statistics.c171
2 files changed, 144 insertions, 31 deletions
diff --git a/src/compare_hkl.c b/src/compare_hkl.c
index b2b40da1..f87a7b88 100644
--- a/src/compare_hkl.c
+++ b/src/compare_hkl.c
@@ -142,9 +142,9 @@ int main(int argc, char *argv[])
}
STATUS("%i,%i reflections: %i in common\n", nc1, nc2, ncom);
R2 = stat_r2(ref1, c1, ref2, c2, &scale);
- STATUS("R2 = %5.4f %% (scale=%5.2f)\n", R2*100.0, scale);
+ STATUS("R2 = %5.4f %% (scale=%5.2e)\n", R2*100.0, scale);
Rmerge = stat_rmerge(ref1, c1, ref2, c2, &scale);
- STATUS("Rmerge = %5.4f %% (scale=%5.2f)\n", Rmerge*100.0, scale);
+ STATUS("Rmerge = %5.4f %% (scale=%5.2e)\n", Rmerge*100.0, scale);
if ( outfile != NULL ) {
write_reflections(outfile, NULL, out, NULL, 1, cell, 1);
diff --git a/src/statistics.c b/src/statistics.c
index d78dde9e..da3d95c4 100644
--- a/src/statistics.c
+++ b/src/statistics.c
@@ -16,11 +16,27 @@
#include <math.h>
#include <stdlib.h>
+#include <gsl/gsl_errno.h>
+#include <gsl/gsl_min.h>
#include "statistics.h"
#include "utils.h"
+struct r_params {
+ const double *ref1;
+ const unsigned int *c1;
+ const double *ref2;
+ const unsigned int *c2;
+ int fom;
+};
+
+enum {
+ R_2,
+ R_MERGE,
+};
+
+
double stat_scale_intensity(const double *ref1, const unsigned int *c1,
const double *ref2, const unsigned int *c2)
{
@@ -45,38 +61,39 @@ double stat_scale_intensity(const double *ref1, const unsigned int *c1,
}
-double stat_r2(const double *ref1, const unsigned int *c1,
- const double *ref2, const unsigned int *c2, double *scalep)
+double stat_scale_sqrti(const double *ref1, const unsigned int *c1,
+ const double *ref2, const unsigned int *c2)
{
double top = 0.0;
double bot = 0.0;
- double scale;
int i;
- scale = stat_scale_intensity(ref1, c1, ref2, c2);
- *scalep = scale;
/* Start from i=1 -> skip central beam */
for ( i=1; i<LIST_SIZE; i++ ) {
if ( c1[i] && c2[i] ) {
- double i1, i2;
- i1 = ref1[i] / (scale*(double)c1[i]);
- i2 = ref2[i] / (double)c2[i];
+ double f1, f2;
- top += pow(fabs(i1 - i2), 2.0);
- bot += pow(i1, 2.0);
+ if ( (ref1[i]<0.0) || (ref2[i]<0.0) ) continue;
- } /* else reflection not measured so don't worry about it */
+ f1 = sqrt(ref1[i]) / (double)c1[i];
+ f2 = sqrt(ref2[i]) / (double)c2[i];
+
+ top += f1 * f2;
+ bot += f2 * f2;
+
+ } /* else reflection not common so don't worry about it */
}
- return sqrt(top/bot);
+ return top/bot;
}
-double stat_scale_sqrti(const double *ref1, const unsigned int *c1,
- const double *ref2, const unsigned int *c2)
+static double internal_r2(const double *ref1, const unsigned int *c1,
+ const double *ref2, const unsigned int *c2,
+ double scale)
{
double top = 0.0;
double bot = 0.0;
@@ -87,33 +104,28 @@ double stat_scale_sqrti(const double *ref1, const unsigned int *c1,
if ( c1[i] && c2[i] ) {
- double f1, f2;
-
- if ( (ref1[i]<0.0) || (ref2[i]<0.0) ) continue;
-
- f1 = sqrt(ref1[i]) / (double)c1[i];
- f2 = sqrt(ref2[i]) / (double)c2[i];
+ double i1, i2;
+ i1 = ref1[i] / (scale*(double)c1[i]);
+ i2 = ref2[i] / (double)c2[i];
- top += f1 * f2;
- bot += f2 * f2;
+ top += pow(fabs(i1 - i2), 2.0);
+ bot += pow(i1, 2.0);
- } /* else reflection not common so don't worry about it */
+ } /* else reflection not measured so don't worry about it */
}
- return top/bot;
+ return sqrt(top/bot);
}
-double stat_rmerge(const double *ref1, const unsigned int *c1,
- const double *ref2, const unsigned int *c2, double *scalep)
+static double internal_rmerge(const double *ref1, const unsigned int *c1,
+ const double *ref2, const unsigned int *c2,
+ double scale)
{
double top = 0.0;
double bot = 0.0;
- double scale;
int i;
- scale = stat_scale_sqrti(ref1, c1, ref2, c2);
- *scalep = scale;
/* Start from i=1 -> skip central beam */
for ( i=1; i<LIST_SIZE; i++ ) {
@@ -136,3 +148,104 @@ double stat_rmerge(const double *ref1, const unsigned int *c1,
return 2.0*top/bot;
}
+
+
+static double calc_r(double scale, void *params)
+{
+ struct r_params *rp = params;
+
+ switch ( rp->fom) {
+ case R_MERGE :
+ return internal_rmerge(rp->ref1, rp->c1,
+ rp->ref2, rp->c2, scale);
+ case R_2 :
+ return internal_r2(rp->ref1, rp->c1,
+ rp->ref2, rp->c2, scale);
+ }
+
+ ERROR("No such FoM!\n");
+ abort();
+}
+
+
+static double r_minimised(const double *ref1, const unsigned int *c1,
+ const double *ref2, const unsigned int *c2,
+ double *scalep, int fom)
+{
+ gsl_function F;
+ gsl_min_fminimizer *s;
+ int status;
+ double scale = 1.0;
+ struct r_params rp;
+ int iter = 0;
+
+ rp.ref1 = ref1;
+ rp.ref2 = ref2;
+ rp.c1 = c1;
+ rp.c2 = c2;
+ rp.fom = fom;
+
+ F.function = &calc_r;
+ F.params = &rp;
+
+ s = gsl_min_fminimizer_alloc(gsl_min_fminimizer_brent);
+
+ /* Initial guess */
+ switch ( fom ) {
+ case R_MERGE :
+ scale = stat_scale_sqrti(ref1, c1, ref2, c2);
+ break;
+ case R_2 :
+ scale = stat_scale_intensity(ref1, c1, ref2, c2);
+ 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);
+
+ do {
+
+ double lo, up;
+
+ /* Iterate */
+ gsl_min_fminimizer_iterate(s);
+
+ /* 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 );
+
+ 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;
+ return calc_r(scale, &rp);
+}
+
+
+double stat_rmerge(const double *ref1, const unsigned int *c1,
+ const double *ref2, const unsigned int *c2,
+ double *scalep)
+{
+ return r_minimised(ref1, c1, ref2, c2, scalep, R_MERGE);
+}
+
+
+double stat_r2(const double *ref1, const unsigned int *c1,
+ const double *ref2, const unsigned int *c2,
+ double *scalep)
+{
+ return r_minimised(ref1, c1, ref2, c2, scalep, R_2);
+}