From 189da15810deabd739d7c11c6e95fea55739fe60 Mon Sep 17 00:00:00 2001
From: Thomas White <taw@physics.org>
Date: Sat, 1 Aug 2020 15:13:49 +0200
Subject: Initial import from archive

---
 src/statistics.c | 223 +++++++++++++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 223 insertions(+)
 create mode 100644 src/statistics.c

(limited to 'src/statistics.c')

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;
+
+}
-- 
cgit v1.2.3