#ifndef TMINIMIZER_HXX_SEEN
#define TMINIMIZER_HXX_SEEN
////////////////////////////////////////////////////////////
// $Id: IMinimizer.hxx,v 1.7 2010/06/10 20:45:50 mcgrew Exp $


#include <exception>
#include <string>
#include <vector>

#include <TObject.h>

class TVector3;
class TLorentzVector;

#include <EoaUtility.hxx>

namespace COMET {
    class IMinimizer;
    class IMinimizerFunc;

    OA_EXCEPTION(EMinimizer,EoaUtility);
    /// Exception thrown when there is a parameter count mismatch.
    OA_EXCEPTION(EParameterCount,EMinimizer);
};

/// An abstract base class to simplify access to function minimization
/// routines in T2K software (see \ref minimization).  Notice that several
/// methods are provided as simple stubs and must be over-ridden in the
/// derived classes.  Classes that don't implement the stubs
/// (e.g. FixParameter()) will simply not have the associated capability.
///
/// Function minimizers are in several places by the T2K event reconstruction
/// software, and it is helpful to have a single simple interface that can be
/// used by the entire library.  This A.B.C. (Abstract Base Class) defines a
/// simplified interface that allows minimizers to be interchanged.  See \ref
/// minimization for specific information about how a minimization routine
/// should be used in a routine.
///
/// This class also helps to seperate the T2K software from the specifics of
/// the ROOT minimization routines and provides a generic interface to be used
/// throughout the T2K libraries.  It also corrects some of the misfeatures
/// found in the current generation of ROOT and guarrantees that the
/// minimization interface will remain stable.

class COMET::IMinimizer {
protected:
    IMinimizerFunc& fFunction;
    
    /// The method to minimize IMinimizerFunc that must be implemented by the
    /// concrete class.  See TMinuitMinimizer for an example implementation.
    virtual int Minimize(int npar, double *par) = 0;

public:
    IMinimizer(IMinimizerFunc& func);
    virtual ~IMinimizer();

    /// Get the current estimated minimum value for the function.
    virtual double GetMinimum(void) const = 0;

    /// Get the status of the last run of the minimizer.
    virtual int GetStatus(void) const = 0;

    /// Get the current estimated value of a parameter.  This will only have a
    /// meaningful value after a minimization has been run.
    virtual double GetParameter(int i) const = 0;

    /// Get the i,j element of the covariance matrix
    virtual double GetCovariance(int, int) const {return 0.0;}

    /// Mark a parameter as fixed at its input or current value.
    virtual void FixParameter(int) {}

    /// Get the precision of the minimum that will be found.  See each fitter
    /// for specifics about what this means, but it is approximately
    /// fractional error in the distance to minimum.
    virtual double GetPrecision(void) const {return -1;}

    /// Set the precision for the minimum to be found. See each fitter for
    /// specifics about what this means, but it is approximately the
    /// fractional error in the distance to minimum.
    virtual void SetPrecision(double accuracy) {}

    /// Get the number of iterations done in the last minimization.
    virtual int GetIterations(void) const {return 0;}

    /// Return the IMinimizerFunc that is being minimized.
    IMinimizerFunc& GetFunction(void) const {return fFunction;};

    /// Set the maximum number of iterations for the minimization.
    virtual void SetMaxIterations(int) {};

    /// The minimizer starts at the initial values of par and returns the
    /// values of the parameters at the minimum.
    void Optimize(int npar, double *par);

    /// Convenience function for minimizing with respect to a single
    /// parameter.
    void Optimize(double &par1);

    /// Convenience function for minimizing with respect to two
    /// parameters.
    void Optimize(double &par1,double &par2);

    /// Convenience function for minimizing with respect to three
    /// parameters.
    void Optimize(double &par1,double &par2,double &par3);

    /// Convenience function for minimizing with respect to a TVector3.  In
    /// this case, the TVector3 elements will be passed to the IMinimizerFunc
    /// as the first, through third parameters.
    void Optimize(TVector3& par1);

    /// Convenience function for minimizing with respect to a TLorentzVector.
    /// In this case, the TLorentzVector elements will be passed to the
    /// IMinimizerFunc as the first, through fourth parameters.
    void Optimize(TLorentzVector& par1);

    /// Convenience function for minimizing with respect to two TVector3
    /// objects.  In this case, the first TVector3 object elements will be
    /// passed to the IMinimizerFunc as the first through third parameters and
    /// the second TVector3 object elements will be passed as the fourth
    /// through sixth parameters.
    void Optimize(TVector3& par1, TVector3& par2);

    /// Convenience function for minimizing with respect to a TLorentzVector
    /// and a TVector3.  The TLorentzVector object elements will be passed to
    /// the IMinimizerFunc as the first through fourth parameters and the
    /// TVector3 object elements will be passed as the fifth through seventh
    /// parameters.
    void Optimize(TLorentzVector& par1, TVector3& par2);
};




/// An abstract base class for a function (functor) to be minimized by
/// IMinimizer.
class COMET::IMinimizerFunc: public TObject {
public:
    /// The object must be created with a specific number of parameters.
    explicit IMinimizerFunc(int nPar);
    virtual ~IMinimizerFunc();
    
    /// Return the number of parameters for this function.
    virtual int GetParameterCount(void) const {return fNParameters;};

    /// Initialize the function to it's initial state.  This is called by the
    /// minimizer immediately before the fit is started.  The initial input
    /// parameters are in par and this must return the value of the function.
    virtual double Initialize(const double* par) {return Evaluate(par);};
    
    /// Calculate the value of the function and if grad is non-null, return
    /// the value of the gradient.
    virtual double Evaluate(const double* par, double* grad=NULL) = 0;

    /// Called after IMinimizer has finished the fit.  This can be used to
    /// collect statistics about the fit, summarize the results, etc.  This
    /// must return the value of the function for the parameter values.
    virtual double Finalize(const double* par) {return Evaluate(par);};

    /// Set the scale for a parameter.  The parameter scale can be used to
    /// provide hints to the minimizer about the shape of the function.  The
    /// scale is not precisely defined, but it is used by the minimizer to
    /// decide how big a step should be used when searching for the minimum.
    /// The scale should be large enough so that when the function is
    /// calculated at two points separated by one unit of the scale, the
    /// magnitude will be "substantially" different so that the first
    /// derivative can be estimated.  However, the scale should be small
    /// enough so that the second derivative has not changed appreciably
    /// between the two points.  Usually, the scale should be about the same
    /// size as the "sigma" of the function.  For instance, y = x^2 and y =
    /// exp(x^2) should have a scale of about "1", while y = exp(x^2/50)
    /// should have a scale of about "0.02".
    void SetScale(int i, double scale) {fScale[i] = scale;};

    /// Get the scale for a parameter.
    double GetScale(int i) {return fScale[i];};

    // Set the name of a parameter
    void SetParameterName(int i, const char* name) {fNames[i] = name;};

    const char* GetParameterName(int i) {return fNames[i].c_str();};

private:
    /// The number of parameters for this function.
    int fNParameters;

    /// A pointer to an array of parameter scales.
    std::vector<double> fScale;

    std::vector<std::string> fNames;
};

/*! \page minimization Function Optimization, Minimization and Fitting.

The oaUtility library provides IMinimizer which is an abstract base class
(A.B.C.)  defining a standard interface for function minimization routines
and all "user" routines should access function minimization through this
A.B.C.

The best way to understand how to use IMinimizer is to study a trivial code
snippet like the one below.  This code will minimize a one dimensional
function (TTrivialFunc, we'll return to it in a moment) and report the
best fit parameters and the minimum function value.

\code 
TTrivialFunc trivialFunc(3.14159);
IMinimizer* trivial = new TMinuitMinimizer(trivialFunc);
double a = 0;
trivial->Optimize(a)
if (!trivial->GetStatus()) {
  std::cout << "Optimal parameter is " << a << std::endl;
  std::cout << "Minimum function value is " << trivial->GetMinimum();
}
else {
  std::cout << "Function not optimized" << std::endl;
}
\endcode

A IMinimizer object, \var trivial, is created by constructing new concrete
TMinuitMinimizer, passing a function derived from IMinimizerFunc as the
argument.  Notice that this code doesn't provide the number of parameters
since that is provided by TTrivialFunc when it is constructed.  You can
find the number of parameters by using TTrivialFunc::GetParamterCount(),
but in general you have this information.  You should always use pointers
to the base class IMinimizer so that you can easily change the algorithm
used to optimize the function.  Having created a IMinimizer, the next line
declares a variable for the parameter to be optimized.  This is followed by
a call to the IMinimizer::Optimize() method.  This method is heavily
overloaded and provides several ways to pass the values to be optimized.
After the variable, \var a, is optimized the status is checked with
IMinimizer::GetStatus(), and the optimal value is displayed using
IMinimizer::GetMinimum(). 

Now that you've seen the basic structure of code using the IMinimizer
object, let's consider how to declare the function to be minimized.  The
function to be minimized must be derived from the IMinimizerFunc base
class, and (besides the constructor and destructor), only needs to provide
one method, IMinimizerFunc::Evaluate.  Here is an example of a very trivial
function which implements a parabola with a minimum at x=fA and a minimum
value of zero.  

\code 
class TTrivialFunc: public IMinimizerFunc {
private:
   double fA;
public:
   TTrivialFunc(double a) : fA(a), IMinimizerFunc(1) {};
   ~TTrivialFunc {};
   Evaluate(const double* par, double* grad=NULL) {
      return (par[0]-fA)*(par[0]-fA);
   }
};
\endcode

Notice that the function can keep local data.  In this case, TTrivialFunc
is saving the location of the minimum, but it can be used to store any data
required to calculate the function.  This provides a simple way to pass
information to the function being minimized.  The object also saves the
number of dimensions for this function in a field inherited from
IMinimizerFunc.  The IMinimizerFunc constructor takes a single
parameter which is the number of dimensions in the function being
minimized.  This is the only place that the number of dimensions can be
specified, and you can access it using
IMinimizerFunc::GetParameterCount().  

The IMinimizerFunc::Evaluate method must be overridden by every concrete
function to be minimized.  This method takes two parameters:  The first is
a pointer to an array of parameter values.  The second parameter is a
pointer to an array that be used to return the value of the gradient.
Usually, the grad parameter can be ignored.  The IMinimizerFunc::Evaluate
method must return the value of the function at the parameter values.

There is a useful code idiom that can be used with classes
derived from IMinimizerFunc to pass information into the function.
Consider the case where a function is calculating a goodness for a
COMET::IHitSelection to find an event vertex.  Clearly, the IMinimizer isn't going
to understand the hit selection, and only cares about the goodness at each
potential vertex, but the IMinimizerFunc requires the COMET::IHitSelection.  This
can be handled by a private field.

\code
class TGoodFunc: public IMinimizerFunc {
private:
   COMET::IHitSelection* fHits;
public:
   ...
   void SetHitSelection(COMET::IHitSelection* hits) {fHits = hits};
   ...
};
\endcode

The hit selection can then be passed to the function using the
SetHitSelection method by directly using the object (goodFunc), or by
accessing the function through IMinimizer::GetFunction() and then using
dynamic_cast. 

\code 
...
TVector3 vertex;
TGoodFunc goodFunc();
IMinimizer* goodness = new TMinuitMinimizer(goodFunc);
...
goodFunc.SetHitSelection(hits);
goodness->Optimize(vertex)
...
dynamic_cast<TGoodFunc&>(goodness->GetFunction()).SetHitSelection(hits);
goodness->Optimize(vertex)
...
\endcode
 */

#endif