#ifndef _utl_Minuit_h_
#define _utl_Minuit_h_

#include <sstream>
#include <TMinuit.h>
#include <utl/AugerException.h>


namespace utl {


  /**
    \class Minuit utl::Minuit "utl/Minuit.h"

    \brief Wrapper class for ROOT TMinuit

    Should be used by passing Function object that derives from
    utl::MinuitFunction helper class and implements all missing
    methods.

    \author Darko Veberic
    \date 25 Aug 2009
    \version $Id$
    \ingroup math
  */

  template<class Function>
  class Minuit {
  public:
    typedef std::pair<double, double> ValueErrorPair;

    Minuit(Function& function)
    {
      if (fFunction)
        throw DoesNotComputeException("Two concurrent instances of the same "
                                      "utl::Minuit<> class not allowed.");
      fFunction = &function;
      const int nPar = fFunction->GetNumParameters();
      fMinuit = new TMinuit(nPar);
      fMinuit->SetFCN(MinuitFunction);
      SetPrintoutLevel(-1);
      fErrorFlag = 0;
      for (int i = 0; i < nPar; ++i) {
        try {
          fMinuit->mnparm(i, fFunction->GetParameterName(i),
                          fFunction->GetInitialValue(i), fFunction->GetInitialStep(i),
                          fFunction->GetLowerBound(i), fFunction->GetUpperBound(i), fErrorFlag);
        } catch (...) {
          Clear();
          throw;
        }
        if (fErrorFlag) {
          std::ostringstream err;
          err << "TMinuit initialization error for parameter " << i << '.';
          Clear();
          throw InvalidConfigurationException(err.str());
        }
        if (fFunction->IsParameterFixed(i))
          fMinuit->FixParameter(i);
      }
      SetErrorDef(fFunction->IsLikelihood() ? 0.5 : 1);
    }

    ~Minuit() { Clear(); }

    bool
    Minimize(const int maxSteps = 500)
    {
      fArgumentList[0] = maxSteps;
      fErrorFlag = 0;
      fMinuit->mnexcm("MINIMIZE", fArgumentList, 1, fErrorFlag);
      return !fErrorFlag;
    }

    bool
    Minos(const int maxSteps = 500)
    {
      fArgumentList[0] = maxSteps;
      fErrorFlag = 0;
      fMinuit->mnexcm("MINOS", fArgumentList, 0, fErrorFlag);
      return !fErrorFlag;
    }

    ValueErrorPair
    GetParameter(const int i)
    {
      if (i < 0 || i >= fFunction->GetNumParameters()) {
        std::ostringstream err;
        err << "Parameter index " << i << " out of bounds.";
        Clear();
        throw OutOfBoundException(err.str());
      }
      ValueErrorPair p;
      fMinuit->GetParameter(i, p.first, p.second);
      return p;
    }

    void
    SetPrintoutLevel(const int level)
    {
      fMinuit->SetPrintLevel(level);
      if (level < 0) {
        fErrorFlag = 0;
        fMinuit->mnexcm("SET NOWARNINGS", fArgumentList, 0, fErrorFlag);
      }
    }

    int GetErrorFlag() const { return fErrorFlag; }

    void SetErrorDef(const double up)
    { fMinuit->SetErrorDef(up); }

    double
    Call(const int flag)
    {
      fArgumentList[0] = flag;
      fErrorFlag = 0;
      fMinuit->mnexcm("CALI", fArgumentList, 1, fErrorFlag);
      return GetMinimum();
    }

    void GetCovariance(double** matrix)
    { fMinuit->mnemat(matrix, fFunction->GetNumParameters()); }

    double GetMinimum()
    { return fMinuit->fAmin; }

  private:
    // prevent copying
    Minuit(const Minuit&);
    Minuit& operator=(const Minuit&);

    void Clear() { fFunction = 0; delete fMinuit; fMinuit = 0; }

    static void MinuitFunction(int& /*nPar*/, double* const /*grad*/,
                               double& value, double* const par,
                               const int /*iFlag*/)
    { value = fFunction->operator()(par); }

    static Function* fFunction;
    TMinuit* fMinuit;
    double fArgumentList[1];
    int fErrorFlag;
  };

  template<class Function>
  Function* Minuit<Function>::fFunction = 0;


}


#endif