/*
 * refine-brent.c
 *
 * Refinement by Sequential Brent Minimisation
 *
 * (c) 2006-2007 Thomas White <taw27@cam.ac.uk>
 *
 *  synth2d - Two-Dimensional Crystallographic Fourier Synthesis
 *
 */

#ifdef HAVE_CONFIG_H
#include <config.h>
#endif

#include <gsl/gsl_errno.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_min.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_linalg.h>

#include "refine.h"
#include "model.h"
#include "reflist.h"
#include "statistics.h"

/* For the purposes of GSL */
static double refine_calc_r(double loc, void *params) {

	double r;
	ReflectionList *model_reflections;
	RefinementPair *pair = params;
	AtomicModel *model;
	
	model = model_copy(pair->model);
	
	switch ( pair->cur_mod_coor ) {
		case COORDINATE_X : model->atoms[pair->cur_mod_atom].x = loc; break;
		case COORDINATE_Y : model->atoms[pair->cur_mod_atom].y = loc; break;
		case COORDINATE_Z : model->atoms[pair->cur_mod_atom].z = loc; break;
	}
	
	model_reflections = model_calculate_f(pair->reflections, model, 69);
	if ( pair->spec & REFINE_SPEC_INTENSITIES ) {
		r = stat_r2(pair->reflections, model_reflections);	/* stat_r(Fobs, Fcalc) */
	} else {
		r = stat_r(pair->reflections, model_reflections);	/* stat_r(Fobs, Fcalc) */
	}
	free(model_reflections);
	
	model_free(model);
	
	return r;
	
}

void refine_brent(AtomicModel *model, ReflectionList *reflections, RefinementSpec spec) {

	gsl_function F;
	gsl_min_fminimizer *s;
	int status;
	RefinementPair pair;
	
	pair.reflections = reflections;
	pair.model = model;
	pair.cur_mod_coor = COORDINATE_X;
	pair.cur_mod_atom = 1;
	
	F.function = &refine_calc_r;
	F.params = &pair;
	
	s = gsl_min_fminimizer_alloc(gsl_min_fminimizer_brent);

	while ( model->refine_window->run_semaphore ) {
	
		double loc_cur;
		unsigned int iter = 0;
	
		/* Select the next thing to refine */
		pair.cur_mod_coor++;
		if ( pair.cur_mod_coor > COORDINATE_Z ) {
			pair.cur_mod_coor = COORDINATE_X;
			pair.cur_mod_atom++;
			if ( pair.cur_mod_atom >= pair.model->n_atoms ) pair.cur_mod_atom = 1;
		}
	//	if ( random() < (0.3333*RAND_MAX) ) pair.cur_mod_coor = COORDINATE_Y;
	//	if ( random() > (0.6666*RAND_MAX) ) pair.cur_mod_coor = COORDINATE_Z;
	//	atom = 1 + (random()*(double)(pair.model->n_atoms-1))/RAND_MAX;
	//	pair.cur_mod_atom = atom;
		
		loc_cur = 0;
		switch ( pair.cur_mod_coor ) {
			case COORDINATE_X : loc_cur = pair.model->atoms[pair.cur_mod_atom].x; break;
			case COORDINATE_Y : loc_cur = pair.model->atoms[pair.cur_mod_atom].y; break;
			case COORDINATE_Z : loc_cur = pair.model->atoms[pair.cur_mod_atom].z; break;
		}
		
		gsl_min_fminimizer_set(s, &F, loc_cur, loc_cur-0.1, loc_cur+0.1);
		
		do {
		
			double lo, up;
		
			/* Iterate */		
			gsl_min_fminimizer_iterate(s);
			iter++;
			
			/* Get the current estimate */
			loc_cur = 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);
			
		} while ( status == GSL_CONTINUE );
		
		if (status == GSL_SUCCESS) {
			
			if ( loc_cur < 0 ) loc_cur = 1+loc_cur;
			if ( loc_cur >= 1 ) loc_cur = loc_cur-1;
									
			switch ( pair.cur_mod_coor ) {
				case COORDINATE_X : pair.model->atoms[pair.cur_mod_atom].x = loc_cur; break;
				case COORDINATE_Y : pair.model->atoms[pair.cur_mod_atom].y = loc_cur; break;
				case COORDINATE_Z : pair.model->atoms[pair.cur_mod_atom].z = loc_cur; break;
			}
			
			refine_schedule_update(model);
								
		}
	
	}
	
	gsl_min_fminimizer_free(s);

}