/***************************************************************************
 *    aon_maximize.cpp - wrapper around gsl minimize functions
 *
 *  Begun      Thu Apr 15 00:48:47 2004
 *  Copyright  2004  Andrew C. O'Neill
 *  Email      oneill@phys.columbia.edu
 ****************************************************************************/

/*
 *  This program is free software; you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation; either version 2 of the License, or
 *  (at your option) any later version.
 *
 *  This program is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU Library General Public License for more details.
 *
 *  You should have received a copy of the GNU General Public License
 *  along with this program; if not, write to the Free Software
 *  Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
 */
 

#include "aon_maximize.h"

#include "gmf_pairproduction.h"

Maximization::Maximization()
{
	dout(MEM) << "Constructing a maximization" << endl;
//	const gsl_min_fminimizer_type *T = gsl_min_fminimizer_brent;
	const gsl_min_fminimizer_type *T = gsl_min_fminimizer_goldensection;
	itsMinimizer = gsl_min_fminimizer_alloc (T);
	
	itsStatus = NOT_READY;
	MaxIterations = 1000;
	itsPrecision = 0.01;
	itsInvertedParams = new InvertedGSLFunctionParams;
	dout(MEM) << "Newed InverteGSLFunctionParams" << endl;
	
}

Maximization::~Maximization()
{
	dout(MEM) << "Destructing a maximization" << endl;
	gsl_min_fminimizer_free(itsMinimizer);	
	delete itsInvertedParams;
	dout(MEM) << "deleted InverteGSLFunctionParams" << endl;
}




int Maximization::setFunction(double (*OriginalFunction)(double, void*), 
                              void* OriginalParams)
{
	// set un-inverted function pointer as given
	if (OriginalFunction != NULL)
	{	
		itsInvertedParams->OriginalFunction = OriginalFunction;
		itsInvertedParams->OriginalParams = OriginalParams;
		
		// setup gsl_function to point to inverted function
		itsGSLFunction.function = &InvertedGSLFunction;
		itsGSLFunction.params = itsInvertedParams;
		
		itsStatus = READY;
		return EXIT_SUCCESS;
	}
	else 
	{
		// do sod all!
		dout(ERR) << "Trying to set a NULL function pointer to be maximized"
		          << endl;		
		itsStatus = NOT_READY;
		return EXIT_FAILURE; 
	}
	
}


double Maximization::FindMaximum(const double& Guess, const double& Lower, 
                                 const double& Upper)
{
	if (itsStatus == NOT_READY)
	{
		dout(ERR) <<  "Trying to use uninitialized maximization" << endl;
		return 0.0;
	}
	else if (itsStatus == READY)
	{
		gsl_min_fminimizer_set(itsMinimizer, &itsGSLFunction, Guess, 
		                                                        Lower, Upper);
		int status;
		double NewGuess = 0.0, NewUpper = 0.0, NewLower = 0.0;
		int Iter = 0;
		do
		{
			Iter++;
			status = gsl_min_fminimizer_iterate (itsMinimizer);

			NewGuess = gsl_min_fminimizer_x_minimum (itsMinimizer);
			NewLower = gsl_min_fminimizer_x_lower (itsMinimizer);
			NewUpper = gsl_min_fminimizer_x_upper (itsMinimizer);
		
			status = gsl_min_test_interval (NewLower, NewUpper, 0.0, itsPrecision);
			
			if (status == GSL_SUCCESS)
      	{
				return -gsl_min_fminimizer_f_minimum(itsMinimizer);
			}

		} while (status == GSL_CONTINUE && Iter < MaxIterations);
		dout(ERR) << "Maximization didn't converge" << endl;
		itsStatus = NOT_READY;
		return 0.0;
	}
	else // status = DONE, no need to find minimum
	{
		return - gsl_min_fminimizer_f_minimum(itsMinimizer);
	}
		
}

double InvertedGSLFunction (double Argument, void* AllParams)
{
	InvertedGSLFunctionParams MyParams = *(InvertedGSLFunctionParams *)AllParams;
	return - (* MyParams.OriginalFunction)(Argument, MyParams.OriginalParams);	
}

double QuickMax(const double& x, const double& y)
{
	if (x > y)
	{
		return x;
	}
	else
	{
		return y;
	}
}