/*
 * refine-rns.c
 *
 * Refinement by Random Neighbour Search
 *
 * (c) 2006-2009 Thomas White <taw27@cam.ac.uk>
 *
 *  synth2d - Two-Dimensional Crystallographic Fourier Synthesis
 *
 */

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

#include <stdio.h>
#include <stdlib.h>

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

void refine_neighboursearch(AtomicModel *model, ReflectionList *reflections,
			    RefinementSpec spec, double shift)
{
	unsigned int iterations_without_change;
	double min_r;
	ReflectionList *model_reflections;

	model_reflections = model_calculate_f(reflections, model, 69);
	if ( spec & REFINE_SPEC_INTENSITIES ) {
		/* stat_r2(Fobs, Fcalc) */
		min_r = stat_r2(reflections, model_reflections);
	} else {
		/* stat_r(Fobs, Fcalc) */
		min_r = stat_r(reflections, model_reflections);
	}

	iterations_without_change = 0;
	while ( (model->refine_window->run_semaphore)
	     && (iterations_without_change < MAX_REFINEMENT_ITERATIONS) ) {

		unsigned int i;
		double x, y, z;
		double new_r;
		AtomicModel *copy;

		/* Copy the structure */
		copy = model_copy(model);

		/* 'Jiggle' the atoms */
		for ( i=1; i<model->n_atoms; i++ ) {

			if ( model->atoms[i].refine ) {

				x = model->atoms[i].x;
				y = model->atoms[i].y;
				z = model->atoms[i].z;

				if ( spec & REFINE_SPEC_X ) {
					if ( random() > (0.9*RAND_MAX) )
						x += shift;
					if ( random() > (0.9*RAND_MAX) )
						x -= shift;
					if ( x<0 ) x = 1+x;
					if ( x>1 ) x = x-1;
					model->atoms[i].x = x;
				}

				if ( spec & REFINE_SPEC_Y ) {
					if ( random() > (0.9*RAND_MAX) )
						y += shift;
					if ( random() > (0.9*RAND_MAX) )
						y -= shift;
					if ( y<0 ) y = 1+y;
					if ( y>1 ) y = y-1;
					model->atoms[i].y = y;
				}

				if ( spec & REFINE_SPEC_Z ) {
					if ( random() > (0.9*RAND_MAX) )
						z += shift;
					if ( random() > (0.9*RAND_MAX) )
						z -= shift;
					if ( z<0 ) z = 1+z;
					if ( z>1 ) z = z-1;
					model->atoms[i].z = z;
				}

			}

		}

		/* Recalculate structure factors and check for improvement */
		model_reflections = model_calculate_f(reflections, model, 69);
		if ( spec & REFINE_SPEC_INTENSITIES ) {
			/* stat_r(Fobs, Fcalc) */
			new_r = stat_r2(reflections, model_reflections);
		} else {
			/* stat_r(Fobs, Fcalc) */
			new_r = stat_r(reflections, model_reflections);
		}
		if ( new_r < min_r ) {

			min_r = new_r;
			refine_schedule_update(model);
			iterations_without_change = 0;

			model_free(copy);

		} else {

			model_move(model, copy);
			model_free(copy);
			iterations_without_change++;

		}
		free(model_reflections);

	}
}