/* ======================================================================== *\
!
! *
! * This file is part of MARS, the MAGIC Analysis and Reconstruction
! * Software. It is distributed to you in the hope that it can be a useful
! * and timesaving tool in analysing Data of imaging Cerenkov telescopes.
! * It is distributed WITHOUT ANY WARRANTY.
! *
! * Permission to use, copy, modify and distribute this software and its
! * documentation for any purpose is hereby granted without fee,
! * provided that the above copyright notice appear in all copies and
! * that both that copyright notice and this permission notice appear
! * in supporting documentation. It is provided "as is" without express
! * or implied warranty.
! *
!
!
!   Author(s): Wolfgang Wittek, 7/2003 <mailto:wittek@mppmu.mpg.de>
!   Author(s): Thomas Bretz, 12/2003 <mailto:tbretz@astro.uni-wuerzburg.de>
!
!   Copyright: MAGIC Software Development, 2000-2003
!
!
\* ======================================================================== */

/////////////////////////////////////////////////////////////////////////////
//
// MTMinuit
//
// Class for interfacing with Minuit
//
/////////////////////////////////////////////////////////////////////////////
#include "MTMinuit.h"

#include <TMinuit.h>

#include "MLog.h"
#include "MLogManip.h"

ClassImp(MTMinuit);

using namespace std;

// --------------------------------------------------------------------------
//
// Default constructor.
//
MTMinuit::MTMinuit(const char *name, const char *title)
    : fFcn(NULL), fNames(NULL), fMinuit(NULL)
{
    fName  = name  ? name  : "MTMinuit";
    fTitle = title ? title : "Interface for Minuit";
}

MTMinuit::~MTMinuit()
{
    if (fMinuit)
        delete fMinuit;
}

TMinuit *MTMinuit::GetMinuit() const
{
    return fMinuit;
}

void MTMinuit::SetFcn(void (*fcn)(Int_t &, Double_t *, Double_t &, Double_t *, Int_t))
{
    fFcn = fcn;
}

void MTMinuit::InitParameters(const TArrayD &vinit, const TArrayD *step, const TString *name)
{
    const int n = fVinit.GetSize();

    fVinit = vinit;

    fStep.Set(n);
    memset(fStep.GetArray(), 0, n*sizeof(Double_t));

    if (step && step->GetSize()!=n)
    {
        *fLog << warn << "MTMinuit::SetParameters: number of steps doesn't match number of parameters... ignored." << endl;
        step = NULL;
    }

    for (int i=0; i<n; i++)
        fStep[i] = step ? (*step)[i] : TMath::Abs(fVinit[i]/10.0);

    if (fNames)
        delete fNames;

    fNames = new TString[n];

    // FIXME: Sanity check for TString array size!
    for (int i=0; i<n ; i++)
        fNames[i] = name ? (const char*)name[i] : Form("par%02d", i);

    // ---------------------------------------------------------------

    fLimLo.Set(n);
    fLimUp.Set(n);
    memset(fLimLo.GetArray(), 0, n*sizeof(Double_t));
    memset(fLimUp.GetArray(), 0, n*sizeof(Double_t));

    fFix.Set(n);
    memset(fFix.GetArray(), 0, n*sizeof(Char_t));
}

void MTMinuit::SetLimits(const TArrayD &limlo, const TArrayD &limup)
{
    if (fVinit.GetSize()==limlo.GetSize())
        fLimLo = limlo;
    else
        *fLog << warn << "MTMinuit::SetLimits: size of limlo doesn't match number of parameters... ignored." << endl;

    if (fVinit.GetSize()==limup.GetSize())
        fLimUp = limup;
    else
        *fLog << warn << "MTMinuit::SetLimits: size of limup doesn't match number of parameters... ignored." << endl;
}

void MTMinuit::SetFixedParameters(const TArrayC &fix)
{
    if (fVinit.GetSize()==fix.GetSize())
        fFix = fix;
    else
        *fLog << warn << "MTMinuit::SetFixedParameters: size of fix doesn't match number of parameters... ignored." << endl;
}

// -----------------------------------------------------------------------
//
// Interface to MINUIT
//
Bool_t MTMinuit::CallMinuit(TObject *objectfit , const TString &method, Bool_t nulloutput)
{
    if (!fFcn(NULL))
    {
        *fLog << err << "CallMinuit: ERROR - fFcn not set... abort." << endl;
        return kFALSE;
    }
    if (!fNames)
    {
        *fLog << err << "CallMinuit: ERROR - Parameters not set... abort." << endl;
        return kFALSE;
    }

    // Make a copy of fStep
    TArrayD step(fStep);

    // Save gMinuit
    TMinuit *const save = gMinuit;

    // Set the maximum number of parameters
    if (fMinuit)
        delete fMinuit;
    fMinuit = new TMinuit(fVinit.GetSize());

    //
    // Set the print level
    // -1   no output except SHOW comands
    //  0   minimum output
    //  1   normal output (default)
    //  2   additional ouput giving intermediate results
    //  3   maximum output, showing progress of minimizations
    //
    fMinuit->SetPrintLevel(-1);

    //..............................................
    // Printout for warnings
    //    SET WAR      print warnings
    //    SET NOW      suppress warnings
    Int_t errWarn;
    Double_t tmpwar = 0;
    fMinuit->mnexcm("SET NOW", &tmpwar, 0, errWarn);
    //fMinuit->mnexcm("SET WAR", &tmpwar, 0, errWarn);

    //..............................................
    // Set the address of the minimization function
    fMinuit->SetFCN(fFcn);

    //..............................................
    // Store address of object to be used in fcn
    fMinuit->SetObjectFit(objectfit);

    //..............................................
    // Set starting values and step sizes for parameters
    for (Int_t i=0; i<fVinit.GetSize(); i++)
    {
        if (!fMinuit->DefineParameter(i, fNames[i], fVinit[i], fStep[i], fLimLo[i], fLimUp[i]))
            continue;

        *fLog << err << "CallMinuit: Error in defining parameter " << fNames[i] << endl;
        return kFALSE;
    }

    //
    // Error definition :
    //
    //    for chisquare function :
    //      up = 1.0   means calculate 1-standard deviation error
    //         = 4.0   means calculate 2-standard deviation error
    //
    //    for log(likelihood) function :
    //      up = 0.5   means calculate 1-standard deviation error
    //         = 2.0   means calculate 2-standard deviation error
    //
    fMinuit->SetErrorDef(1.0);

    // Int_t errMigrad;
    // Double_t tmp = 0;
    // fMinuit->mnexcm("MIGRAD", &tmp, 0, errMigrad);

    // fix a parameter
    for (Int_t i=0; i<fVinit.GetSize(); i++)
    {
        if (!fFix[i])
            continue;

        fMinuit->FixParameter(i);
    }

    //..............................................
    // This doesn't seem to have any effect
    // Set maximum number of iterations (default = 500)
    //Int_t maxiter = 100000;
    //fMinuit->SetMaxIterations(maxiter);

    // minimization by the method of Migrad
    Int_t fErrMinimize=0;
    if (method.Contains("Migrad", TString::kIgnoreCase))
    {
        if (nulloutput)
            fLog->SetNullOutput(kTRUE);
        Double_t tmp = 0;
        fMinuit->mnexcm("MIGRAD", &tmp, 0, fErrMinimize);
        if (nulloutput)
            fLog->SetNullOutput(kFALSE);
    }

    //..............................................
    // same minimization as by Migrad
    // but switches to the SIMPLEX method if MIGRAD fails to converge
    if (method.Contains("Minimize", TString::kIgnoreCase))
    {
        Double_t tmp = 0;
        fMinuit->mnexcm("MINIMIZE", &tmp, 0, fErrMinimize);
    }

    //..............................................
    // minimization by the SIMPLEX method
    if (method.Contains("Simplex", TString::kIgnoreCase))
    {
        Double_t tmp[2];
        tmp[0] = 3000; // maxcalls;
        tmp[1] = 0.1;  // tolerance;

        if (nulloutput)
            fLog->SetNullOutput(kTRUE);
        fMinuit->mnexcm("SIMPLEX", &tmp[0], 2, fErrMinimize);
        if (nulloutput)
            fLog->SetNullOutput(kFALSE);
    }

    //..............................................
    // check quality of minimization
    // istat = 0   covariance matrix not calculated
    //         1   diagonal approximation only (not accurate)
    //         2   full matrix, but forced positive-definite
    //         3   full accurate covariance matrix
    //             (indication of normal convergence)
    Double_t fMin,   fEdm,   fErrdef;
    Int_t    fNpari, fNparx, fIstat;
    fMinuit->mnstat(fMin, fEdm, fErrdef, fNpari, fNparx, fIstat);

    if (fErrMinimize != 0)
    {
        *fLog << err << "CallMinuit: Minimization failed with:" << endl;
        *fLog << "  fMin = " << fMin << endl;
        *fLog << "  fEdm = "  << fEdm << endl;
        *fLog << "  fErrdef = "  << fErrdef << endl;
        *fLog << "  fIstat = " << fIstat << endl;
        *fLog << "  fErrMinimize = " << fErrMinimize << endl;
        return kFALSE;
    }

    //..............................................
    // minimization by the method of Migrad
    if (method.Contains("Hesse", TString::kIgnoreCase))
    {
        Double_t tmp = 0;
        fMinuit->mnexcm("HESSE", &tmp, 0, fErrMinimize);
    }

    //..............................................
    // Minos error analysis
    if (method.Contains("Minos", TString::kIgnoreCase))
    {
        Double_t tmp = 0;
        fMinuit->mnexcm("MINOS", &tmp, 0, fErrMinimize);
    }

    //..............................................
    // Print current status of minimization
    // if nkode = 0    only function value
    //            1    parameter values, errors, limits
    //            2    values, errors, step sizes, internal values
    //            3    values, errors, step sizes, 1st derivatives
    //            4    values, parabolic errors, MINOS errors

    //Int_t nkode = 4;
    //fMinuit->mnprin(nkode, fmin);

    //..............................................
    // call fcn with IFLAG = 3 (final calculation : calculate p(chi2))
    // iflag = 1   initial calculations only
    //         2   calculate 1st derivatives and function
    //         3   calculate function only
    //         4   calculate function + final calculations
    Double_t iflag = 3;
    Int_t errfcn3;
    fMinuit->mnexcm("CALL", &iflag, 1, errfcn3);

    // WW : the following statements were commented out because the
    // Minuit object will still be used;
    // this may be changed in the future 
    gMinuit = save;

    return kTRUE;
}