1 files changed, 138 insertions, 0 deletions
diff --git a/src/refine-brent.c b/src/refine-brent.c
new file mode 100644
index 0000000..06a6f95
--- /dev/null
+++ b/src/refine-brent.c
@@ -0,0 +1,138 @@
+/*
+ * 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);
+
+}
+