diff options
Diffstat (limited to 'src/statistics.c')
| -rw-r--r-- | src/statistics.c | 223 |
1 files changed, 223 insertions, 0 deletions
diff --git a/src/statistics.c b/src/statistics.c new file mode 100644 index 0000000..5fa4634 --- /dev/null +++ b/src/statistics.c @@ -0,0 +1,223 @@ +/* + * statistics.c + * + * Structure Factor Statistics + * + * (c) 2006 Thomas White <taw27@cam.ac.uk> + * Synth2D - two-dimensional Fourier synthesis + * + */ + +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif + +#include <math.h> +#include <stdio.h> +#include <gsl/gsl_errno.h> +#include <gsl/gsl_math.h> +#include <gsl/gsl_min.h> +#include <assert.h> + +#include "reflist.h" +#include "statistics.h" + +/* Return the least-squares estimate of the optimum scaling factor for intensities */ +double stat_scale_intensity(ReflectionList *obs, ReflectionList *calc) { + + unsigned int i; + double top = 0; + double bot = 0; + + if ( obs->n_reflections == 0 ) return 0; + if ( calc->n_reflections == 0 ) return 0; /* No reflections */ + + for ( i=1; i<obs->n_reflections; i++ ) { /* 'hkl' loop */ + + assert(obs->refs[i].amplitude >= 0); + assert(calc->refs[i].amplitude >= 0); + + top += (obs->refs[i].amplitude*obs->refs[i].amplitude) * (calc->refs[i].amplitude*calc->refs[i].amplitude); + bot += (calc->refs[i].amplitude*calc->refs[i].amplitude) * (calc->refs[i].amplitude*calc->refs[i].amplitude); + + } + + return top/bot; + +} + +/* Return the least-squares estimate of the optimum scaling factor */ +double stat_scale(ReflectionList *obs, ReflectionList *calc) { + + unsigned int i; + double top = 0; + double bot = 0; + + if ( obs->n_reflections == 0 ) return 0; + if ( calc->n_reflections == 0 ) return 0; /* No reflections */ + + for ( i=1; i<obs->n_reflections; i++ ) { /* 'hkl' loop */ + + assert(obs->refs[i].amplitude >= 0); + assert(calc->refs[i].amplitude >= 0); + + top += obs->refs[i].amplitude * calc->refs[i].amplitude; + bot += calc->refs[i].amplitude * calc->refs[i].amplitude; + + } + + return top/bot; + +} + +/* R-factor in terms of diffracted intensities */ +double stat_r2(ReflectionList *obs, ReflectionList *calc) { + + unsigned int i; + double scale; + double err; + double den; + + scale = stat_scale_intensity(obs, calc); + err = 0; den = 0; + + for ( i=1; i<obs->n_reflections; i++ ) { /* 'hkl' loop */ + + assert(obs->refs[i].amplitude >= 0); + assert(calc->refs[i].amplitude >= 0); + + err += fabs( (obs->refs[i].amplitude*obs->refs[i].amplitude) - (scale * (calc->refs[i].amplitude*calc->refs[i].amplitude)) ); + den += obs->refs[i].amplitude * obs->refs[i].amplitude; + + } + + return err/den; + +} + +/* R-factor in terms of amplitudes of structure factors */ +double stat_r1(ReflectionList *obs, ReflectionList *calc) { + + unsigned int i; + double scale; + double err; + double den; + + scale = stat_scale(obs, calc); + err = 0; den = 0; + + for ( i=1; i<obs->n_reflections; i++ ) { /* 'hkl' loop */ + + assert(obs->refs[i].amplitude >= 0); + assert(calc->refs[i].amplitude >= 0); + + err += fabs( obs->refs[i].amplitude - (scale * calc->refs[i].amplitude) ); + den += obs->refs[i].amplitude; + + } + + return err/den; + +#if 0 + double scale; + gsl_function F; + gsl_min_fminimizer *s; + int status; + int iter = 0, max_iter = 100; + + printf("Estimate: Scale = %f, R=%.2f%%\n", opt.scale, opt.r*100); + + opt.r = 999999999; + opt.scale = 0; + + F.function = &stat_calc_r; + F.params = &pair; + + s = gsl_min_fminimizer_alloc(gsl_min_fminimizer_brent); + gsl_min_fminimizer_set(s, &F, 10, 1, 1000); + + do { + + double lo, up; + + /* Iterate */ + gsl_min_fminimizer_iterate(s); + iter++; + + /* 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); + if (status == GSL_SUCCESS) { + + OptimumR opt; + + opt.r = stat_calc_r(scale, &pair); + opt.scale = scale; + + printf("Minimum: Scale=%f, R=%.2f%%\n", opt.scale, opt.r*100); + + gsl_min_fminimizer_free(s); + + return opt; + + } + + } while (status == GSL_CONTINUE && iter < max_iter); + + return opt; +#endif + +} + +double stat_sigma_f(ReflectionList *reflections) { + + unsigned int i; + double sigma_f = 0; + + if ( reflections->n_reflections == 0 ) return 0; /* No reflections */ + + for ( i=1; i<reflections->n_reflections; i++ ) { /* 'hkl' loop */ + assert(reflections->refs[i].amplitude >= 0); + sigma_f += reflections->refs[i].amplitude; + } + + return sigma_f; + +} + +double stat_stddev(ReflectionList *a) { + + unsigned int i; + double sigma_f, mean; + double sigma_dev = 0; + + if ( a->n_reflections == 0 ) return 0; /* No reflections */ + + sigma_f = stat_sigma_f(a); + mean = sigma_f / a->n_reflections; + for ( i=1; i<a->n_reflections; i++ ) { /* 'hkl' loop */ + sigma_dev += (a->refs[i].amplitude - mean) * (a->refs[i].amplitude - mean); + } + + return sqrt(sigma_dev/a->n_reflections); + +} + +double stat_maxam(ReflectionList *a) { + + unsigned int i; + double max_am = 0; + + if ( a->n_reflections == 0 ) return 0; /* No reflections */ + + for ( i=1; i<a->n_reflections; i++ ) { + if ( a->refs[i].amplitude > max_am ) max_am = a->refs[i].amplitude; + } + + return max_am; + +} |