#include "TMinuitMinimizer.h"
#include "Math/IFunction.h"
#include "TMinuit.h"
#include "TGraph.h"
#include "TError.h"
#include <iostream>
#include <cassert>
#include <algorithm>
#include <functional>
#include <cmath>
#define USE_STATIC_TMINUIT
ROOT::Math::IMultiGenFunction * TMinuitMinimizer::fgFunc = 0;
TMinuit * TMinuitMinimizer::fgMinuit = 0;
bool TMinuitMinimizer::fgUsed = false;
ClassImp(TMinuitMinimizer)
TMinuitMinimizer::TMinuitMinimizer(ROOT::Minuit::EMinimizerType type ) :
fUsed(false),
fDim(0),
fStrategy(1),
fType(type),
fMinuit(fgMinuit)
{
}
TMinuitMinimizer::TMinuitMinimizer(const char * type ) :
fUsed(false),
fDim(0),
fStrategy(1),
fMinuit(fgMinuit)
{
std::string algoname(type);
std::transform(algoname.begin(), algoname.end(), algoname.begin(), (int(*)(int)) tolower );
ROOT::Minuit::EMinimizerType algoType = ROOT::Minuit::kMigrad;
if (algoname == "simplex") algoType = ROOT::Minuit::kSimplex;
if (algoname == "minimize" ) algoType = ROOT::Minuit::kCombined;
if (algoname == "migradimproved" ) algoType = ROOT::Minuit::kMigradImproved;
if (algoname == "scan" ) algoType = ROOT::Minuit::kScan;
if (algoname == "seek" ) algoType = ROOT::Minuit::kSeek;
fType = algoType;
}
TMinuitMinimizer::~TMinuitMinimizer()
{
#ifndef USE_STATIC_TMINUIT
if (fMinuit) delete fMinuit;
#endif
}
TMinuitMinimizer::TMinuitMinimizer(const TMinuitMinimizer &) :
Minimizer()
{
}
TMinuitMinimizer & TMinuitMinimizer::operator = (const TMinuitMinimizer &rhs)
{
if (this == &rhs) return *this;
return *this;
}
void TMinuitMinimizer::SetFunction(const ROOT::Math::IMultiGenFunction & func) {
fDim = func.NDim();
#ifdef USE_STATIC_TMINUIT
if (fgMinuit == 0) {
fgUsed = false;
fgMinuit = new TMinuit(fDim);
}
else if (fgMinuit->GetNumPars() != int(fDim) ) {
delete fgMinuit;
fgUsed = false;
fgMinuit = new TMinuit(fDim);
}
fMinuit = fgMinuit;
#else
if (fMinuit) {
delete fMinuit;
}
fMinuit = new TMinuit(fDim);
#endif
fDim = func.NDim();
fgFunc = const_cast<ROOT::Math::IMultiGenFunction *>(&func);
fMinuit->SetFCN(&TMinuitMinimizer::Fcn);
double arglist[1];
arglist[0] = PrintLevel() - 1;
int ierr= 0;
fMinuit->mnexcm("SET PRINT",arglist,1,ierr);
}
void TMinuitMinimizer::SetFunction(const ROOT::Math::IMultiGradFunction & func) {
fDim = func.NDim();
#ifdef USE_STATIC_TMINUIT
if (fgMinuit == 0) {
fgUsed = false;
fgMinuit = new TMinuit(fDim);
}
else if (fgMinuit->GetNumPars() != int(fDim) ) {
delete fgMinuit;
fgUsed = false;
fgMinuit = new TMinuit(fDim);
}
fMinuit = fgMinuit;
#else
if (fMinuit) delete fMinuit;
fMinuit = new TMinuit(fDim);
#endif
fDim = func.NDim();
fgFunc = const_cast<ROOT::Math::IMultiGradFunction *>(&func);
fMinuit->SetFCN(&TMinuitMinimizer::FcnGrad);
double arglist[1];
arglist[0] = PrintLevel() - 1;
int ierr= 0;
fMinuit->mnexcm("SET PRINT",arglist,1,ierr);
fMinuit->mnexcm("SET GRAD",arglist,0,ierr);
}
void TMinuitMinimizer::Fcn( int &, double * , double & f, double * x , int ) {
f = fgFunc->operator()(x);
}
void TMinuitMinimizer::FcnGrad( int &, double * g, double & f, double * x , int iflag ) {
ROOT::Math::IMultiGradFunction * gFunc = dynamic_cast<ROOT::Math::IMultiGradFunction *> ( fgFunc);
assert(gFunc != 0);
f = gFunc->operator()(x);
if (iflag == 2) gFunc->Gradient(x,g);
}
bool TMinuitMinimizer::SetVariable(unsigned int ivar, const std::string & name, double val, double step) {
if (fMinuit == 0) {
Error("SetVariable","invalid TMinuit pointer. Need to call first SetFunction");
return false;
}
#ifdef USE_STATIC_TMINUIT
fUsed = fgUsed;
#endif
if (fUsed) DoClear();
fMinuit->DefineParameter(ivar , name.c_str(), val, step, 0., 0. );
return true;
}
bool TMinuitMinimizer::SetLimitedVariable(unsigned int ivar, const std::string & name, double val, double step, double lower, double upper) {
if (fMinuit == 0) {
Error("SetVariable","invalid TMinuit pointer. Need to call first SetFunction");
return false;
}
#ifdef USE_STATIC_TMINUIT
fUsed = fgUsed;
#endif
if (fUsed) DoClear();
fMinuit->DefineParameter(ivar, name.c_str(), val, step, lower, upper );
return true;
}
#ifdef LATER
bool Minuit2Minimizer::SetLowerLimitedVariable(unsigned int ivar , const std::string & name , double val , double step , double lower ) {
if (fMinuit == 0) {
Error("SetVariable","invalid TMinuit pointer. Need to call first SetFunction");
return false;
}
#ifdef USE_STATIC_TMINUIT
fUsed = fgUsed;
#endif
if (fUsed) DoClear();
double s = val-lower;
double upper = s*1.0E15;
if (s != 0) upper = 1.0E15;
return SetLimitedVariable(ivar, name, val, step, lower,upper);
}
#endif
bool TMinuitMinimizer::SetFixedVariable(unsigned int ivar, const std::string & name, double val) {
if (fMinuit == 0) {
Error("SetVariable","invalid TMinuit pointer. Need to call first SetFunction");
return false;
}
#ifdef USE_STATIC_TMINUIT
fUsed = fgUsed;
#endif
if (fUsed) DoClear();
double step = ( val != 0) ? 0.1 * std::abs(val) : 0.1;
fMinuit->DefineParameter(ivar, name.c_str(), val, step, 0., 0. );
fMinuit->FixParameter(ivar);
return true;
}
bool TMinuitMinimizer::SetVariableValue(unsigned int ivar, double val) {
if (fMinuit == 0) {
Error("SetVariable","invalid TMinuit pointer. Need to call first SetFunction");
return false;
}
double arglist[2];
int ierr = 0;
arglist[0] = ivar+1;
arglist[1] = val;
fMinuit->mnexcm("SET PAR",arglist,2,ierr);
return (ierr==0);
}
bool TMinuitMinimizer::Minimize() {
if (fMinuit == 0) {
Error("Minimize","invalid TMinuit pointer. Need to call first SetFunction and SetVariable");
return false;
}
if (fMinuit->fNu < static_cast<int>(fDim) ) {
Error("Minimize","The total number of defined parameters is different than the function dimension, npar = %d, dim = %d",fMinuit->fNu, fDim);
return false;
}
double arglist[10];
int ierr = 0;
arglist[0] = ErrorDef();
fMinuit->mnexcm("SET Err",arglist,1,ierr);
int printlevel = PrintLevel();
arglist[0] = printlevel - 1;
fMinuit->mnexcm("SET PRINT",arglist,1,ierr);
if (printlevel == 0) fMinuit->mnexcm("SET NOW",arglist,0,ierr);
arglist[0] = MaxFunctionCalls();
arglist[1] = Tolerance();
int nargs = 2;
switch (fType){
case ROOT::Minuit::kMigrad:
fMinuit->mnexcm("MIGRAD",arglist,nargs,ierr);
break;
case ROOT::Minuit::kCombined:
fMinuit->mnexcm("MINIMIZE",arglist,nargs,ierr);
break;
case ROOT::Minuit::kSimplex:
fMinuit->mnexcm("SIMPLEX",arglist,nargs,ierr);
break;
case ROOT::Minuit::kScan:
nargs = 0;
fMinuit->mnexcm("SCAN",arglist,nargs,ierr);
break;
case ROOT::Minuit::kSeek:
nargs = 1;
if (arglist[1] >= 1.) nargs = 2;
fMinuit->mnexcm("SEEK",arglist,nargs,ierr);
break;
default:
fMinuit->mnexcm("MIGRAD",arglist,nargs,ierr);
}
#ifdef USE_STATIC_TMINUIT
fgUsed = true;
#endif
fUsed = true;
fStatus = ierr;
int minErrStatus = ierr;
if (ierr == 0 && fType == ROOT::Minuit::kMigradImproved) {
fMinuit->mnexcm("IMPROVE",arglist,1,ierr);
fStatus += 1000*ierr;
}
if (ierr == 0 && IsValidError() ) {
fMinuit->mnexcm("HESSE",arglist,1,ierr);
fStatus += 100*ierr;
}
RetrieveParams();
if (minErrStatus == 0) {
RetrieveErrorMatrix();
fMinosRun = false;
return true;
}
return false;
}
void TMinuitMinimizer::RetrieveParams() {
assert(fMinuit != 0);
if (fParams.size() != fDim) fParams.resize( fDim);
if (fErrors.size() != fDim) fErrors.resize( fDim);
for (unsigned int i = 0; i < fDim; ++i) {
fMinuit->GetParameter( i, fParams[i], fErrors[i]);
}
}
void TMinuitMinimizer::RetrieveErrorMatrix() {
assert(fMinuit != 0);
unsigned int nfree = NFree();
unsigned int ndim2 = fDim*fDim;
if (fCovar.size() != ndim2 ) fCovar.resize(fDim*fDim);
if (nfree >= fDim) {
fMinuit->mnemat(&fCovar.front(), fDim);
}
else {
std::vector<double> tmpMat(nfree*nfree);
fMinuit->mnemat(&tmpMat.front(), nfree);
unsigned int l = 0;
for (unsigned int i = 0; i < fDim; ++i) {
if ( fMinuit->fNiofex[i] > 0 ) {
unsigned int m = 0;
for (unsigned int j = 0; j <= i; ++j) {
if ( fMinuit->fNiofex[j] > 0 ) {
fCovar[i*fDim + j] = tmpMat[l*nfree + m];
fCovar[j*fDim + i] = fCovar[i*fDim + j];
m++;
}
}
l++;
}
}
}
}
unsigned int TMinuitMinimizer::NCalls() const {
if (fMinuit == 0) return 0;
return fMinuit->fNfcn;
}
double TMinuitMinimizer::MinValue() const {
if (!fMinuit) return 0;
double minval = fMinuit->fAmin;
if (minval == fMinuit->fUndefi) return 0;
return minval;
}
double TMinuitMinimizer::Edm() const {
if (!fMinuit) return -1;
if (fMinuit->fAmin == fMinuit->fUndefi || fMinuit->fEDM == fMinuit->fBigedm) return fMinuit->fUp;
return fMinuit->fEDM;
}
unsigned int TMinuitMinimizer::NFree() const {
if (!fMinuit) return 0;
if (fMinuit->fNpar < 0) return 0;
return fMinuit->fNpar;
}
int TMinuitMinimizer::CovMatrixStatus() const {
if (!fMinuit) return 0;
if (fMinuit->fAmin == fMinuit->fUndefi) return 0;
return fMinuit->fISW[1];
}
double TMinuitMinimizer::GlobalCC(unsigned int i) const {
if (!fMinuit) return 0;
if (!fMinuit->fGlobcc) return 0;
return fMinuit->fGlobcc[i];
}
bool TMinuitMinimizer::GetMinosError(unsigned int i, double & errLow, double & errUp) {
if (fMinuit == 0) {
Error("GetMinosError","invalid TMinuit pointer. Need to call first SetFunction and SetVariable");
return false;
}
double arglist[2];
int ierr = 0;
if (!fMinosRun) {
arglist[0] = ErrorDef();
fMinuit->mnexcm("SET Err",arglist,1,ierr);
arglist[0] = PrintLevel()-1;
fMinuit->mnexcm("SET PRINT",arglist,1,ierr);
if (PrintLevel() == 0) fMinuit->mnexcm("SET NOW",arglist,0,ierr);
}
arglist[0] = MaxFunctionCalls();
arglist[1] = i+1;
int nargs = 2;
fMinuit->mnexcm("MINOS",arglist,nargs,ierr);
fStatus += 10*ierr;
fMinosRun = true;
double errParab = 0;
double gcor = 0;
fMinuit->mnerrs(i,errUp,errLow, errParab, gcor);
if (fStatus%100 != 0 ) return false;
return true;
}
void TMinuitMinimizer::DoClear() {
fMinuit->mncler();
double val = 3;
int inseed = 12345;
fMinuit->mnrn15(val,inseed);
fUsed = false;
#ifdef USE_STATIC_TMINUIT
fgUsed = false;
#endif
}
void TMinuitMinimizer::PrintResults() {
if (fMinuit == 0) return;
if (PrintLevel() > 2)
fMinuit->mnprin(4,fMinuit->fAmin);
else
fMinuit->mnprin(3,fMinuit->fAmin);
}
bool TMinuitMinimizer::Contour(unsigned int ipar, unsigned int jpar, unsigned int &npoints, double * x, double * y) {
if (fMinuit == 0) {
Error("TMinuitMinimizer::Contour"," invalid TMinuit instance");
return false;
}
double arglist[1];
int ierr = 0;
arglist[0] = ErrorDef();
fMinuit->mnexcm("SET Err",arglist,1,ierr);
arglist[0] = PrintLevel()-1;
fMinuit->mnexcm("SET PRINT",arglist,1,ierr);
if (PrintLevel() == 0) fMinuit->mnexcm("SET NOW",arglist,0,ierr);
if (npoints < 4) {
Error("Contour","Cannot make contour with so few points");
return false;
}
int npfound = 0;
npoints -= 1;
fMinuit->mncont( ipar,jpar,npoints, x, y,npfound);
if (npfound<4) {
Error("Contour","Cannot find more than 4 points");
return false;
}
if (npfound!=(int)npoints) {
Warning("Contour","Returning only %d points ",npfound);
npoints = npfound;
}
return true;
}
bool TMinuitMinimizer::Scan(unsigned int ipar, unsigned int & nstep, double * x, double * y, double xmin, double xmax) {
if (fMinuit == 0) {
Error("TMinuitMinimizer::Scan"," invalid TMinuit instance");
return false;
}
if (xmin >= xmax && (int) ipar < fMinuit->GetNumPars() ) {
double val = 0; double err = 0;
TString name;
double xlow = 0; double xup = 0 ;
int iuint = 0;
fMinuit->mnpout( ipar, name, val, err, xlow, xup, iuint);
if (iuint > 0 && err > 0) {
xmin = val - 2.*err;
xmax = val + 2 * err;
}
}
double arglist[4];
int ierr = 0;
arglist[0] = PrintLevel()-1;
fMinuit->mnexcm("SET PRINT",arglist,1,ierr);
if (PrintLevel() == 0) fMinuit->mnexcm("SET NOW",arglist,0,ierr);
if (nstep == 0) return false;
arglist[0] = ipar+1;
arglist[1] = nstep+2;
int nargs = 2;
if (xmax > xmin ) {
arglist[2] = xmin;
arglist[3] = xmax;
nargs = 4;
}
fMinuit->mnexcm("SCAN",arglist,nargs,ierr);
if (ierr) {
Error("TMinuitMinimizer::Scan"," Error executing command SCAN");
return false;
}
TGraph * gr = dynamic_cast<TGraph *>(fMinuit->GetPlot() );
if (!gr) {
Error("TMinuitMinimizer::Scan"," Error in returned graph object");
return false;
}
nstep = std::min(gr->GetN(), (int) nstep);
std::copy(gr->GetX(), gr->GetX()+nstep, x);
std::copy(gr->GetY(), gr->GetY()+nstep, y);
nstep = gr->GetN();
return true;
}
bool TMinuitMinimizer::Hesse() {
if (fMinuit == 0) {
Error("Hesse","invalid TMinuit pointer. Need to call first SetFunction and SetVariable");
return false;
}
double arglist[10];
int ierr = 0;
arglist[0] = ErrorDef();
fMinuit->mnexcm("SET ERR",arglist,1,ierr);
int printlevel = PrintLevel();
arglist[0] = printlevel - 1;
fMinuit->mnexcm("SET PRINT",arglist,1,ierr);
if (printlevel == 0) fMinuit->mnexcm("SET NOW",arglist,0,ierr);
arglist[0] = MaxFunctionCalls();
fMinuit->mnexcm("HESSE",arglist,1,ierr);
fStatus += 100*ierr;
if (ierr != 0) return false;
RetrieveParams();
RetrieveErrorMatrix();
return true;
}