/////////////////////////////////////////////////////////////////////////////
// opt-simann.cc
//
// SIMLIB version: 2.18
// Date: 2004-01-25
//
// Copyright (c) 2000-2004 Petr Peringer 
//
// This library is licensed under GNU Library GPL. See the file COPYING.
//

// EXPERIMENTAL
// simulated annealing optimization method

#include "simlib.h"
#include "optimize.h"		// Param, ParameterVector
#include <cmath>		// exp()

#define debug 1			// 0=NO, 1=OPTn, >1=ALL

//////////////////////////////////////////////////////////////////////////////
// optimization method
//

///////////////////////////////////////////////////////////////////////////////
// decision based on temperature
//
bool accept_bad(double eps)
{
    double p = exp(6.0 * (eps - 1.0));	// go down exponentialy
    double x = Random();	// <0,1)
    return (x <= p);
}

///////////////////////////////////////////////////////////////////////////////
// next random point
void move_to_next_point(ParameterVector & p, double eps)
{
    for (int i = 0; i < p.size(); i++) {
	double range = p[i].Range();
	double delta = range * (Random() - 0.5) * eps;
	p[i] = p[i] + delta;	// limited assignment
    }
}

///////////////////////////////////////////////////////////////////////////////
// simulated annealing method
//
double Optimize_simann(double (*f) (const ParameterVector & p),
		       ParameterVector & p, int MAXT)
{
    int bad_count = 0;
    double opt = 1e30;
    double xopt = opt;
    ParameterVector px(p);
    bad_count = 0;		// statistic
    for (int temp = MAXT; temp > 0; temp--) {	// falling temperature 
	// random step depends on temperature
	double eps = temp / double (MAXT);
	ParameterVector new_p = px;
	move_to_next_point(new_p, eps);
	// evaluate cost function
	double new_x = f(new_p);
#if debug>1			// ALL points
	Print("%g %g %.12g\n", new_p["d"].Value(), new_p["k"].Value(), new_x);
#endif
	bool bad = false;
	if (new_x < xopt || (bad = accept_bad(eps))) {	// move to next point
	    xopt = new_x;	// can be worse 
	    px = new_p;
	    if (bad)
		++bad_count;
#if debug==3
	    Print("# %shil step accepted: ", bad ? "up" : "down");
	    Print("%g %g %.12g\n", new_p["d"].Value(), new_p["k"].Value(),
		  new_x);
#endif
	}
	if (new_x < opt) {	// accept new temporary optimum
	    opt = new_x;
	    p = new_p;
#if debug==1			// optima only
	    Print("%g %g %.12g\n", p["d"].Value(), p["k"].Value(), opt);
#endif
	}
    }
#if debug
    Print("# %d accepted uphill steps\n", bad_count);
#endif
    return opt;			// return optimal value and p
}

// end