#include "TROOT.h"
#include "TFitterMinuit.h"
#include "TF1.h"
#include "TH1.h"
#include "TGraph.h"
#include "TChi2FCN.h"
#include "TChi2ExtendedFCN.h"
#include "TBinLikelihoodFCN.h"
#include "Minuit2/MnMigrad.h"
#include "Minuit2/MnMinos.h"
#include "Minuit2/MnHesse.h"
#include "Minuit2/MinuitParameter.h"
#include "Minuit2/MnPrint.h"
#include "Minuit2/FunctionMinimum.h"
#include "Minuit2/VariableMetricMinimizer.h"
#include "Minuit2/SimplexMinimizer.h"
#include "Minuit2/CombinedMinimizer.h"
#include "Minuit2/ScanMinimizer.h"
#include <iomanip>
using namespace ROOT::Minuit2;
#ifndef ROOT_TMethodCall
#include "TMethodCall.h"
#endif
#include "Api.h"
#ifdef OLDFCN_INTERFACE
extern void H1FitChisquare(Int_t &npar, Double_t *gin, Double_t &f, Double_t *u, Int_t flag);
extern void H1FitLikelihood(Int_t &npar, Double_t *gin, Double_t &f, Double_t *u, Int_t flag);
extern void Graph2DFitChisquare(Int_t &npar, Double_t *gin, Double_t &f, Double_t *u, Int_t flag);
extern void MultiGraphFitChisquare(Int_t &npar, Double_t *gin, Double_t &f, Double_t *u, Int_t flag);
#endif
ClassImp(TFitterMinuit);
TFitterMinuit* gMinuit2 = 0;
TFitterMinuit::TFitterMinuit() : fErrorDef(1.) , fEDMVal(0.), fGradient(false),
fState(MnUserParameterState()), fMinosErrors(std::vector<MinosError>()), fMinimizer(0), fMinuitFCN(0), fDebug(1), fStrategy(1), fMinTolerance(0) {
Initialize();
}
TFitterMinuit::TFitterMinuit(Int_t ) : fErrorDef(1.) , fEDMVal(0.), fGradient(false), fState(MnUserParameterState()), fMinosErrors(std::vector<MinosError>()), fMinimizer(0), fMinuitFCN(0), fDebug(1), fStrategy(1), fMinTolerance(0) {
Initialize();
}
void TFitterMinuit::Initialize() {
SetName("Minuit2");
gMinuit2 = this;
gROOT->GetListOfSpecials()->Add(gMinuit2);
}
void TFitterMinuit::CreateMinimizer(EMinimizerType type) {
#ifdef DEBUG
std::cout<<"TFitterMinuit:create minimizer of type "<< type << std::endl;
#endif
if (fMinimizer) delete fMinimizer;
switch (type) {
case kMigrad:
SetMinimizer( new VariableMetricMinimizer() );
return;
case kSimplex:
SetMinimizer( new SimplexMinimizer() );
return;
case kCombined:
SetMinimizer( new CombinedMinimizer() );
return;
case kScan:
SetMinimizer( new ScanMinimizer() );
return;
case kFumili:
std::cout << "TFitterMinuit::Error - Fumili Minimizer must be created from TFitterFumili " << std::endl;
SetMinimizer(0);
return;
default:
SetMinimizer( new VariableMetricMinimizer() );
}
}
TFitterMinuit::~TFitterMinuit() {
#ifdef DEBUG
std::cout << "delete minimizer and FCN" << std::endl;
#endif
if (fMinuitFCN) delete fMinuitFCN;
if (fMinimizer) delete fMinimizer;
}
Double_t TFitterMinuit::Chisquare(Int_t npar, Double_t *params) const {
const TBinLikelihoodFCN * fcn = dynamic_cast<const TBinLikelihoodFCN *> (GetMinuitFCN() );
if (fcn == 0) return 0;
std::vector<double> p(npar);
for (int i = 0; i < npar; ++i)
p[i] = params[i];
return fcn->Chi2(p);
}
void TFitterMinuit::Clear(Option_t*) {
fErrorDef = 1.;
fEDMVal = 0;
fGradient = false;
State() = MnUserParameterState();
fMinosErrors.clear();
fStrategy = 1;
fMinTolerance = 0;
fCovar.clear();
if (fMinimizer) {
delete fMinimizer;
fMinimizer = 0;
}
}
FunctionMinimum TFitterMinuit::DoMinimization( int nfcn, double edmval) {
assert(GetMinuitFCN() != 0 );
assert(GetMinimizer() != 0 );
fMinuitFCN->SetErrorDef(fErrorDef);
if (fDebug >=1) {
std::cout << "TFitterMinuit - Minimize with max iterations = " << nfcn << " edmval = " << edmval << " errorDef = " << fMinuitFCN->Up() << std::endl;
}
return GetMinimizer()->Minimize(*GetMinuitFCN(), State(), MnStrategy(fStrategy), nfcn, edmval);
}
int TFitterMinuit::Minimize( int nfcn, double edmval) {
edmval = std::max(fMinTolerance, edmval);
int prevLevel = gErrorIgnoreLevel;
if (fDebug < 0)
gErrorIgnoreLevel = 1001;
FunctionMinimum min = DoMinimization(nfcn,edmval);
if (fDebug < 0) gErrorIgnoreLevel = prevLevel;
fState = min.UserState();
return ExamineMinimum(min);
}
Int_t TFitterMinuit::ExecuteCommand(const char *command, Double_t *args, Int_t nargs){
#ifdef DEBUG
std::cout<<"Execute command= "<<command<<std::endl;
#endif
if (strncmp(command,"MIG",3) == 0 || strncmp(command,"mig",3) == 0) {
int nfcn = 0;
double edmval = 0.1;
if (nargs > 0) nfcn = int ( args[0] );
if (nargs > 1) edmval = args[1];
CreateMinimizer(kMigrad);
return Minimize(nfcn, edmval);
}
if (strncmp(command,"MINI",4) == 0 || strncmp(command,"mini",4) == 0) {
int nfcn = 0;
double edmval = 0.1;
if (nargs > 0) nfcn = int ( args[0] );
if (nargs > 1) edmval = args[1];
CreateMinimizer(kCombined);
return Minimize(nfcn, edmval);
}
if (strncmp(command,"SIM",3) == 0 || strncmp(command,"sim",3) == 0) {
int nfcn = 0;
double edmval = 0.1;
if (nargs > 0) nfcn = int ( args[0] );
if (nargs > 1) edmval = args[1];
CreateMinimizer(kSimplex);
return Minimize(nfcn, edmval);
}
if (strncmp(command,"SCA",3) == 0 || strncmp(command,"sca",3) == 0) {
int nfcn = 0;
double edmval = 0.1;
if (nargs > 0) nfcn = int ( args[0] );
if (nargs > 1) edmval = args[1];
CreateMinimizer(kScan);
return Minimize(nfcn, edmval);
}
else if (strncmp(command,"MINO",4) == 0 || strncmp(command,"mino",4) == 0) {
int prevLevel = gErrorIgnoreLevel;
if (fDebug < 0)
gErrorIgnoreLevel = 1001;
FunctionMinimum min = DoMinimization();
if (!min.IsValid() ) {
std::cout << "TFitterMinuit::MINOS failed due to invalid function minimum" << std::endl;
if (fDebug < 0) gErrorIgnoreLevel = prevLevel;
return -10;
}
MnMinos minos( *fMinuitFCN, min);
fMinosErrors.clear();
for(unsigned int i = 0; i < State().VariableParameters(); i++) {
if (fDebug>=3) std::cout << "Running Minos for parameter (ext#) " << State().ExtOfInt(i) << std::endl;
MinosError me = minos.Minos(State().ExtOfInt(i));
if (fDebug >= 0) {
if ( !me.IsValid() )
std::cout << "Error running Minos for parameter " << State().ExtOfInt(i) << std::endl;
}
if (fDebug >= 1) {
if (!me.LowerValid() )
std::cout << "Minos: Invalid lower error for parameter " << State().ExtOfInt(i) << std::endl;
if(me.AtLowerLimit())
std::cout << "Minos: Parameter is at Lower limit."<<std::endl;
if(me.AtLowerMaxFcn())
std::cout << "Minos: Maximum number of function calls exceeded when running for lower error" <<std::endl;
if(me.LowerNewMin() )
std::cout << "Minos: New Minimum found while running Minos for lower error" <<std::endl;
if (!me.UpperValid() )
std::cout << "Minos: Invalid upper error for parameter " << State().ExtOfInt(i) << std::endl;
if(me.AtUpperLimit())
std::cout << "Minos: Parameter is at Upper limit."<<std::endl;
if(me.AtUpperMaxFcn())
std::cout << "Minos: Maximum number of function calls exceeded when running for upper error" <<std::endl;
if(me.UpperNewMin() )
std::cout << "Minos: New Minimum found while running Minos for upper error" <<std::endl;
}
fMinosErrors.push_back(me);
}
if (fDebug >= 3) {
for(std::vector<MinosError>::const_iterator ime = fMinosErrors.begin();
ime != fMinosErrors.end(); ime++)
std::cout<<*ime<<std::endl;
}
if (fDebug < 0) gErrorIgnoreLevel = prevLevel;
return 0;
}
else if (strncmp(command,"HES",3) == 0 || strncmp(command,"hes",3) == 0) {
int prevLevel = gErrorIgnoreLevel;
if (fDebug < 0)
gErrorIgnoreLevel = 1001;
MnHesse hesse( GetStrategy() );
const FCNBase * fcn = GetMinuitFCN();
assert(fcn != 0);
fState = hesse(*fcn, State(),100000 );
if (fDebug < 0) gErrorIgnoreLevel = prevLevel;
if (!fState.IsValid() ) {
std::cout << "TFitterMinuit::Hesse calculation failed " << std::endl;
return -10;
}
return 0;
}
else if (strncmp(command,"FIX",3) == 0 || strncmp(command,"fix",3) == 0) {
for(int i = 0; i < nargs; i++) {
FixParameter(int(args[i])-1);
}
return 0;
}
else if (strncmp(command,"SET LIM",7) == 0 || strncmp(command,"set lim",7) == 0) {
assert(nargs >= 3);
int ipar = int(args[0]);
double low = args[1];
double up = args[2];
State().SetLimits(ipar, low, up);
return 0;
}
else if (strncmp(command,"SET PRINT",9) == 0 || strncmp(command,"set print",9) == 0) {
fDebug = int(args[0]);
#ifdef DEBUG
fDebug = 3;
#endif
return 0;
}
else if (strncmp(command,"SET Err",7) == 0 || strncmp(command,"set err",7) == 0) {
fErrorDef = args[0];
return 0;
}
else if (strncmp(command,"SET STR",7) == 0 || strncmp(command,"set str",7) == 0) {
fStrategy = int(args[0]);
return 0;
}
else if (strncmp(command,"SET GRA",7) == 0 || strncmp(command,"set gra",7) == 0) {
return -1;
}
else if (strncmp(command,"CALL FCN",8) == 0 || strncmp(command,"call fcn",8) == 0) {
if (nargs < 1 || fFCN == 0) return -1;
const std::vector<double> & params = State().Params();
std::cout << State() << std::endl;
int npar = params.size();
double fval = 0;
(*fFCN)(npar, 0, fval, const_cast<double *>(¶ms.front()),int(args[0]) ) ;
return 0;
}
else {
return 0;
}
}
int TFitterMinuit::ExamineMinimum(const FunctionMinimum & min) {
if (fDebug >= 3) {
#ifdef LATER
const std::vector<MinimumState>& iterationStates = min.States();
std::cout << "Number of iterations " << iterationStates.size() << std::endl;
for (unsigned int i = 0; i < iterationStates.size(); ++i) {
const MinimumState & st = iterationStates[i];
std::cout << "----------> Iteration " << i << std::endl;
int pr = std::cout.precision(18);
std::cout << " FVAL = " << st.Fval() << " Edm = " << st.Edm() << " Nfcn = " << st.NFcn() << std::endl;
std::cout.precision(pr);
std::cout << " Error matrix change = " << st.Error().Dcovar() << std::endl;
std::cout << " Internal parameters : ";
for (int j = 0; j < st.size() ; ++j) std::cout << " p" << j << " = " << st.Vec()(j);
std::cout << std::endl;
}
#endif
}
if (min.IsValid() ) {
if (fDebug >=1 ) {
std::cout << "Minimum Found" << std::endl;
int pr = std::cout.precision(18);
std::cout << "FVAL = " << State().Fval() << std::endl;
std::cout << "Edm = " << State().Edm() << std::endl;
std::cout.precision(pr);
std::cout << "Nfcn = " << State().NFcn() << std::endl;
std::vector<double> par = State().Params();
std::vector<double> err = State().Errors();
for (unsigned int i = 0; i < State().Params().size(); ++i)
std::cout << State().Parameter(i).Name() << "\t = " << par[i] << "\t +/- " << err[i] << std::endl;
}
return 0;
}
else {
if (fDebug >= 1) {
std::cout << "TFitterMinuit::Minimization DID not converge !" << std::endl;
std::cout << "FVAL = " << State().Fval() << std::endl;
std::cout << "Edm = " << State().Edm() << std::endl;
std::cout << "Nfcn = " << State().NFcn() << std::endl;
}
if (min.HasMadePosDefCovar() ) {
if (fDebug >= 1) std::cout << " Covar was made pos def" << std::endl;
return -11;
}
if (min.HesseFailed() ) {
if (fDebug >= 1) std::cout << " Hesse is not valid" << std::endl;
return -12;
}
if (min.IsAboveMaxEdm() ) {
if (fDebug >= 1) std::cout << " Edm is above max" << std::endl;
return -13;
}
if (min.HasReachedCallLimit() ) {
if (fDebug >= 1) std::cout << " Reached call limit" << std::endl;
return -14;
}
return -10;
}
return 0;
}
void TFitterMinuit::FixParameter(Int_t ipar) {
State().Fix(ipar);
}
Double_t* TFitterMinuit::GetCovarianceMatrix() const {
unsigned int npar = State().Covariance().Nrow();
if ( int(npar) != GetNumberFreeParameters() ) {
std::cout << "TFitterMinuit::GetCovarianceMatrix Error - return null pointer" << std::endl;
return 0;
}
if (fCovar.size() != npar )
fCovar.resize(npar*npar);
for (unsigned int i = 0; i < npar; ++i) {
for (unsigned int j = 0; j < npar; ++j) {
fCovar[j + npar*i] = State().Covariance()(i,j);
}
}
return &(fCovar.front());
}
Double_t TFitterMinuit::GetCovarianceMatrixElement(Int_t i, Int_t j) const {
return State().Covariance()(i,j);
}
Int_t TFitterMinuit::GetErrors(Int_t ipar,Double_t &eplus, Double_t &eminus, Double_t &eparab, Double_t &globcc) const {
eplus = 0.;
eminus = 0.;
MinuitParameter mpar = State().Parameters().Parameter(ipar);
if (mpar.IsFixed() || mpar.IsConst() ) return 0;
if(fMinosErrors.empty()) return 0;
unsigned int nintern = State().IntOfExt(ipar);
eplus = fMinosErrors[nintern].Upper();
eminus = fMinosErrors[nintern].Lower();
eparab = State().Error(ipar);
globcc = State().GlobalCC().GlobalCC()[ipar];
return 0;
}
Int_t TFitterMinuit::GetNumberTotalParameters() const {
return State().Parameters().Parameters().size();
}
Int_t TFitterMinuit::GetNumberFreeParameters() const {
return State().Parameters().VariableParameters();
}
Double_t TFitterMinuit::GetParError(Int_t ipar) const {
return State().Error(ipar);
}
Double_t TFitterMinuit::GetParameter(Int_t ipar) const {
return State().Value(ipar);
}
Int_t TFitterMinuit::GetParameter(Int_t ipar,char *name,Double_t &value,Double_t &verr,Double_t &vlow, Double_t &vhigh) const {
const MinuitParameter& mp = State().Parameter(ipar);
name = const_cast<char*>(mp.Name());
value = mp.Value();
verr = mp.Error();
vlow = mp.LowerLimit();
vhigh = mp.UpperLimit();
return 0;
}
const char *TFitterMinuit::GetParName(Int_t ipar) const {
const MinuitParameter& mp = State().Parameter(ipar);
return mp.Name();
}
Int_t TFitterMinuit::GetStats(Double_t &amin, Double_t &edm, Double_t &errdef, Int_t &nvpar, Int_t &nparx) const {
amin = State().Fval();
edm = State().Edm();
errdef = fErrorDef;
nvpar = State().VariableParameters();
nparx = State().Parameters().Parameters().size();
return 0;
}
Double_t TFitterMinuit::GetSumLog(Int_t) {
return 0.;
}
Bool_t TFitterMinuit::IsFixed(Int_t ipar) const {
return State().Parameter(ipar).IsFixed();
}
void TFitterMinuit::PrintResults(Int_t level, Double_t) const {
if (fDebug >= 0 || level > 3) {
std::cout<<State()<<std::endl;
}
else {
if(!State().IsValid()) {
std::cout << std::endl;
std::cout << "WARNING: Minimum is not valid."<<std::endl;
std::cout << std::endl;
}
std::cout << "# of function calls: "<<State().NFcn()<<std::endl;
std::cout << "function Value: "<< std::setprecision(12) << State().Fval()<<std::endl;
std::cout << "expected distance to the Minimum (edm): "<< std::setprecision(8) << State().Edm()<<std::endl;
std::cout << "fitted parameters: "<<State().Parameters()<<std::endl;
}
if (level > 3) {
for(std::vector<MinosError>::const_iterator ime = fMinosErrors.begin();
ime != fMinosErrors.end(); ime++) {
std::cout<<*ime<<std::endl;
}
}
}
void TFitterMinuit::ReleaseParameter(Int_t ipar) {
State().Release(ipar);
}
void TFitterMinuit::SetFitMethod(const char *name) {
#ifdef DEBUG
std::cout<<"SetFitMethod to "<< name << std::endl;
#endif
if (!strcmp(name,"H1FitChisquare")) {
#ifdef OLDFCN_INTERFACE
SetFCN(H1FitChisquare);
#else
CreateChi2FCN();
#endif
return;
}
if (!strcmp(name,"GraphFitChisquare")) {
#ifdef OLDFCN_INTERFACE
SetFCN(GraphFitChisquare);
#else
if (!GetFitOption().W1)
CreateChi2ExtendedFCN( );
else
CreateChi2FCN( );
#endif
return;
}
if (!strcmp(name, "Graph2DFitChisquare")) {
#ifdef OLDFCN_INTERFACE
SetFCN(Graph2DFitChisquare);
#else
CreateChi2FCN( );
#endif
return;
}
if (!strcmp(name, "MultiGraphFitChisquare")) {
fErrorDef = 1.;
#ifdef OLDFCN_INTERFACE
SetFCN(MultiGraphFitChisquare);
#else
CreateChi2FCN();
#endif
return;
}
if (!strcmp(name, "H1FitLikelihood")) {
#ifdef OLDFCN_INTERFACE
SetFCN(H1FitLikelihood);
#else
CreateBinLikelihoodFCN();
#endif
return;
}
std::cout << "TFitterMinuit::fit method " << name << " is not supported !" << std::endl;
assert(fMinuitFCN != 0);
}
Int_t TFitterMinuit::SetParameter(Int_t,const char *parname,Double_t value,Double_t verr,Double_t vlow, Double_t vhigh) {
#ifdef DEBUG
std::cout<<"SetParameter";
std::cout << parname<<" value = " << value << " verr= "<<verr<<std::endl;
#endif
if(vlow < vhigh) {
State().Add(parname, value, verr, vlow, vhigh);
}
else
State().Add(parname, value, verr);
if (verr == 0) State().Fix(parname);
return 0;
}
void TFitterMinuit::SetFCN(void (*fcn)(Int_t &, Double_t *, Double_t &f, Double_t *, Int_t))
{
fFCN = fcn;
if (fMinuitFCN) delete fMinuitFCN;
fMinuitFCN = new TFcnAdapter(fFCN);
}
void InteractiveFCNm2(Int_t &npar, Double_t *gin, Double_t &f, Double_t *u, Int_t flag)
{
TMethodCall *m = gMinuit2->GetMethodCall();
if (!m) return;
Long_t args[5];
args[0] = (Long_t)∦
args[1] = (Long_t)gin;
args[2] = (Long_t)&f;
args[3] = (Long_t)u;
args[4] = (Long_t)flag;
m->SetParamPtrs(args);
Double_t result;
m->Execute(result);
}
void TFitterMinuit::SetFCN(void *fcn)
{
if (!fcn) return;
char *funcname = G__p2f2funcname(fcn);
if (funcname) {
fMethodCall = new TMethodCall();
fMethodCall->InitWithPrototype(funcname,"Int_t&,Double_t*,Double_t&,Double_t*,Int_t");
}
fFCN = InteractiveFCNm2;
gMinuit2 = this;
if (fMinuitFCN) delete fMinuitFCN;
fMinuitFCN = new TFcnAdapter(fFCN);
}
void TFitterMinuit::SetMinuitFCN( FCNBase * f) {
if (fMinuitFCN) delete fMinuitFCN;
fMinuitFCN = f;
}
void TFitterMinuit::CreateChi2FCN() {
SetMinuitFCN(new TChi2FCN( *this ) );
}
void TFitterMinuit::CreateChi2ExtendedFCN() {
SetMinuitFCN(new TChi2ExtendedFCN( *this ) );
}
void TFitterMinuit::CreateBinLikelihoodFCN() {
SetMinuitFCN(new TBinLikelihoodFCN( *this ) );
}
Last update: Thu Jan 17 08:50:16 2008
This page has been automatically generated. If you have any comments or suggestions about the page layout send a mail to ROOT support, or contact the developers with any questions or problems regarding ROOT.