/*
 * 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;

}