/* ======================================================================== *\
!
! *
! * 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>
!
!   Copyright: MAGIC Software Development, 2000-2003
!
!
\* ======================================================================== */

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

#include <math.h>            // fabs 

#include <TArrayD.h>
#include <TArrayI.h>
#include <TMinuit.h>
#include <TStopwatch.h>

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

ClassImp(MMinuitInterface);

using namespace std;


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

// -----------------------------------------------------------------------
//
// Interface to MINUIT
//
//
Bool_t MMinuitInterface::CallMinuit( 
             void (*fcn)(Int_t &, Double_t *, Double_t &, Double_t *, Int_t),
             const TString* name, const TArrayD &vinit, const TArrayD &step2,
             const TArrayD &limlo, const TArrayD &limup, const TArrayI &fix,
             TObject *objectfit , const TString &method, Bool_t nulloutput)
{
    TArrayD step(step2);

    //
    // Be carefull: This is not thread safe
    //
    if (vinit.GetSize() != step.GetSize()  ||
        vinit.GetSize() != limlo.GetSize() ||
        vinit.GetSize() != limup.GetSize() ||
        vinit.GetSize() != fix.GetSize())
    {
        *fLog << "CallMinuit: Parameter Size Mismatch" << endl;
        return kFALSE;
    }



    //..............................................
    // Set the maximum number of parameters
    //TMinuit *const save = gMinuit;

    TMinuit &minuit = *new TMinuit(vinit.GetSize());
    minuit.SetName("MinuitSupercuts");


    //..............................................
    // 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
    //
    minuit.SetPrintLevel(1);

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

    //..............................................
    // Set the address of the minimization function
    minuit.SetFCN(fcn);

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

    //..............................................
    // Set starting values and step sizes for parameters
    for (Int_t i=0; i<vinit.GetSize(); i++)
    {
        if (minuit.DefineParameter(i, name[i], vinit[i], step[i],
                                   limlo[i], limup[i]))
        {
            *fLog << "CallMinuit: Error in defining parameter "
                << name[i][0] << endl;
            return kFALSE;
        }
    }

    //..............................................
    //Int_t NumPars;
    //NumPars = minuit.GetNumPars();
    //*fLog << "CallMinuit :  number of free parameters = "
    //     << NumPars << endl;


    //..............................................
        // 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
    minuit.SetErrorDef(1.0);

    // Int_t errMigrad;
    // Double_t tmp = 0;
    // minuit.mnexcm("MIGRAD", &tmp, 0, errMigrad);

    //..............................................
    // fix a parameter
    for (Int_t i=0; i<vinit.GetSize(); i++)
    {
        if (fix[i] > 0)
        {
            minuit.FixParameter(i);
            step[i] = 0.0;
        }
    }

    //..............................................
    //NumPars = minuit.GetNumPars();
    //*fLog << "CallMinuit :  number of free parameters = "
    //     << NumPars << endl;

    //..............................................
    // Set maximum number of iterations (default = 500)
    //Int_t maxiter = 100000;
    minuit.SetMaxIterations(10);

    //..............................................
    // minimization by the method of Migrad
    if (method.Contains("Migrad", TString::kIgnoreCase))
    {
        //*fLog << "call MIGRAD" << endl;
        if (nulloutput)
            fLog->SetNullOutput(kTRUE);
        Double_t tmp = 0;
        minuit.mnexcm("MIGRAD", &tmp, 0, fErrMinimize);
        if (nulloutput)
            fLog->SetNullOutput(kFALSE);
        //*fLog << "return from MIGRAD" << endl;
    }

    //..............................................
    // same minimization as by Migrad
    // but switches to the SIMPLEX method if MIGRAD fails to converge
    if (method.Contains("Minimize", TString::kIgnoreCase))
    {
        *fLog << "call MINIMIZE" << endl;
        Double_t tmp = 0;
        minuit.mnexcm("MINIMIZE", &tmp, 0, fErrMinimize);
        *fLog << "return from MINIMIZE" << endl;
    }

    //..............................................
    // minimization by the SIMPLEX method
    if (method.Contains("Simplex", TString::kIgnoreCase))
    {
        *fLog << "call SIMPLEX" << endl;
        if (nulloutput)
            fLog->SetNullOutput(kTRUE);
        Double_t tmp = 0;
        minuit.mnexcm("SIMPLEX", &tmp, 0, fErrMinimize);
        if (nulloutput)
            fLog->SetNullOutput(kFALSE);
        //*fLog << "return from SIMPLEX" << endl;
    }

    //..............................................
    // 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)
    minuit.mnstat(fMin, fEdm, fErrdef, fNpari, fNparx, fIstat);

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

    //*fLog << "CallMinuit : Minimization was successful" << endl;
    //*fLog << "       fMin = " << fMin   << ",   fEdm = "  << fEdm
    //     << ",   fErrdef = "  << fErrdef << ",   fIstat = " << fIstat
    //     << ",   fErrMinimize = " << fErrMinimize << endl;


    //..............................................
    // minimization by the method of Migrad
    if (method.Contains("Hesse", TString::kIgnoreCase))
    {
        //*fLog << "call HESSE" << endl;
        Double_t tmp = 0;
        minuit.mnexcm("HESSE", &tmp, 0, fErrMinimize);
        //*fLog << "return from HESSE" << endl;
    }

    //..............................................
    // Minos error analysis
    if (method.Contains("Minos", TString::kIgnoreCase))
    {
        //*fLog << "call MINOS" << endl;
        Double_t tmp = 0;
        minuit.mnexcm("MINOS", &tmp, 0, fErrMinimize);
        //*fLog << "return from MINOS" << endl;
    }

    //..............................................
    // 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, paraboloc errors, MINOS errors

    //Int_t nkode = 4;
    //minuit.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;
    minuit.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 
    //delete &minuit;
    //gMinuit = save;

    return kTRUE;
}
