/*!\file:  OptimalSearch.cpp
 * \brief optimization algorithm
 */ 

#ifdef HAVE_CONFIG_H
	#include <config.h>
#else
#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
#endif

#include "./numerics.h"
#include "../../classes/objects/objects.h"
#include "../../io/io.h"
#include "../../include/include.h"
#include "../../shared/shared.h"
#include <float.h>

void OptimalSearch(IssmDouble* psearch_scalar,IssmDouble* pJ,OptPars* optpars,IssmDouble (*f)(IssmDouble,OptArgs*), OptArgs* optargs){

	/* This routine is optimizing a given function*/

	/*function values: */
	IssmDouble fx1,fx2,fxbest;
	IssmDouble x1,x2,xmin,xbest;

	/*tolerances: */
	IssmDouble seps;
	IssmDouble tolerance=1.e-4;
	int    maxiter;

	/*counters: */
	int  iter=0;
	bool loop;

	/*Recover parameters:*/
	xmin     =optpars->xmin;
	x1       =optpars->xmin;
	x2       =optpars->xmax;
	maxiter  =optpars->maxiter;

	//get the value of the function at the first boundary
	fx1= (*f)(x1,optargs);
	if (xIsNan<IssmDouble>(fx1)) _error_("Function evaluation returned NaN");
	cout<<setprecision(5);
	if(VerboseControl()) _pprintLine_("");
	if(VerboseControl()) _pprintLine_("       Iteration         x           f(x)       Tolerance");
	if(VerboseControl()) _pprintLine_("");
	if(VerboseControl()) _pprintLine_("           N/A    "<<setw(12)<<x1<<"  "<<setw(12)<<fx1<<"           N/A");

	//update tolerances
	seps=sqrt(DBL_EPSILON);

	loop=true;
	while(loop){

		/*get f(x2)*/
		iter++;
		fx2 = (*f)(x2,optargs);
		if (xIsNan<IssmDouble>(fx2)) _error_("Function evaluation returned NaN");
		if(VerboseControl())
		 _pprintLine_("         "<<setw(5)<<iter<<"    "<<setw(12)<<x2<<"  "<<setw(12)<<fx2<<"  "<<(fabs(x2-x1)>fabs(fx2-fx1)?fabs(fx2-fx1):fabs(x2-x1)));

		//Stop the optimization?
		if ((fabs(x2-x1)+seps)<tolerance || (fabs(fx2-fx1)+seps)<tolerance){
			if(VerboseControl()) _pprintLine_("      " << "optimization terminated: the current x satisfies the termination criteria using 'tolx' of "  << tolerance);
			loop=false;
		}
		else if (iter>=maxiter){
			if(VerboseControl()) _pprintLine_("      " << "exiting: Maximum number of iterations has been exceeded  - increase 'maxiter'");
			loop=false;
		}
		else{
			//continue
			loop=true;
		}

		/*Compute x2*/
		if(fx2<fx1){
			xbest=x2; fxbest=fx2;
			x1=x2;
			x2=xmin+1.1*(x2-xmin); 
			fx1=fx2;
		}
		else{
			xbest=x1; fxbest=fx1;
			x2=xmin+0.5*(x2-xmin);
		}
	}//end while

	/*Assign output pointers: */
	*psearch_scalar=xbest;
	*pJ=fxbest;
}
