/////////////////////////////////////////////////////////////////////////////
// opt-hooke.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
// Hooke-Jeves optimization method

#include "simlib.h"
#include "optimize.h"

#define debug 1			// print values

//## repair
#include <math.h>

//////////////////////////////////////////////////////////////////////////////
// optimization method
//
// look for a better minimum, change one coord at a time 
static double hooke_step(double *delta, opt_function_t f, ParameterVector & p,
			 double min0)
{
    int n = p.size();
    double fmin = min0;
    for (int i = 0; i < n; i++) {
	if (delta[i] == 0.0)
	    continue;		// ignore zero delta
	double old = p[i];
	p[i] = old + delta[i];
	double ftmp = (p[i] == old) ? fmin : f(p);
	if (ftmp < fmin)
	    fmin = ftmp;
	else {
	    delta[i] = -delta[i];	// opposite direction
	    p[i] = old + delta[i];
	    ftmp = (p[i] == old) ? fmin : f(p);
	    if (ftmp < fmin)
		fmin = ftmp;
	    else
		p[i] = old;	// no change
	}
    }
    return fmin;
}

//////////////////////////////////////////////////////////////////////////////

double Optimize_hooke(opt_function_t f, ParameterVector & parameter,
		      double rho, double epsilon, int itermax)
{
// assert(rho>0.01 && rho <1.0);
// assert(epsilon>1e-12 && epsilon < rho);  // 1e-12 > 100*DBL_EPS
// assert(itermax>0);
    int n = parameter.size();	// number of parameters
// assert(n>0);
    double *delta = new double[n];
    ParameterVector oldx(parameter);
    ParameterVector newx(parameter);
    for (int i = 0; i < n; i++)
	delta[i] = fabs(parameter[i].Range() / 10);	// initial step 10==MPARAMETER-???
    int iteration = 0;
    double steplength = rho;	// 1.0 ???
    double newf = f(newx);
#if debug
    newx.PrintValues();
    Print("%g\n", newf);
#endif
    double oldf = newf;
    while (iteration < itermax && steplength > epsilon) {
	iteration++;
	newx = oldx;
	newf = hooke_step(delta, f, newx, oldf);
	// if we made some improvements, continue that direction 
	while (newf < oldf) {
#if debug
	    newx.PrintValues();
	    Print("%g\n", newf);
#endif
	    for (int i = 0; i < n; i++) {
		double dxi = newx[i] - oldx[i];	// actual difference
		// arrange the sign of delta[]
		delta[i] = (dxi <= 0.0) ? -fabs(delta[i]) : fabs(delta[i]);
		// move further in this direction 
		oldx[i] = newx[i];
		newx[i] = newx[i] + dxi;
	    }
	    oldf = newf;
	    newf = hooke_step(delta, f, newx, oldf);
	    /* if the further (optimistic) move was bad.... */
	    if (newf >= oldf)	// worse
		break;		///////////////////// break
	    /* !!!! make sure that the differences between the new and the old
	     * points are due to actual displacements; beware of roundoff
	     * errors that might cause newf < oldf 
	     */
	    int i;
	    for (i = 0; i < n; i++)
		if (fabs(newx[i] - oldx[i]) > (0.5 * fabs(delta[i])))
		    break;
	    // if there is NOT actual change in at least one axis
	    if (i == n)
		break;		///////////////////// break
	}
	if (steplength >= epsilon && newf >= oldf) {
	    steplength *= rho;
	    for (int i = 0; i < n; i++)
		delta[i] *= rho;	// decrease delta
	}
    }
    delete[]delta;		// 
    parameter = oldx;		// copy result
//    iterations = iteration;   // number of iterations
    return oldf;		// return last minimum value
}

// end