From 241cd73749e2c8812090272bc669f777c9460c05 Mon Sep 17 00:00:00 2001
From: Lorenzo Moneta <Lorenzo.Moneta@cern.ch>
Date: Mon, 3 Jul 2006 22:06:42 +0000
Subject: [PATCH] update doc and identation for coding convention
git-svn-id: http://root.cern.ch/svn/root/trunk@15673 27541ba8-7e3a-0410-8455-c3a389f83636
---
.../inc/Minuit2/MnGlobalCorrelationCoeff.h | 5 +-
minuit2/inc/Minuit2/MnPosDef.h | 6 +-
minuit2/inc/Minuit2/MnUserFcn.h | 6 +-
minuit2/src/MnFumiliMinimize.cxx | 26 +-
minuit2/src/MnFunctionCross.cxx | 584 +++++++++---------
minuit2/src/MnGlobalCorrelationCoeff.cxx | 33 +-
minuit2/src/MnHesse.cxx | 532 ++++++++--------
minuit2/src/MnLineSearch.cxx | 4 +-
minuit2/src/MnMachinePrecision.cxx | 69 ++-
minuit2/src/MnMinos.cxx | 254 ++++----
minuit2/src/MnParabolaFactory.cxx | 4 +-
minuit2/src/MnParameterScan.cxx | 87 +--
minuit2/src/MnPlot.cxx | 64 +-
minuit2/src/MnPosDef.cxx | 121 ++--
minuit2/src/MnPrint.cxx | 466 +++++++-------
minuit2/src/MnScan.cxx | 24 +-
minuit2/src/MnSeedGenerator.cxx | 190 +++---
minuit2/src/MnStrategy.cxx | 69 ++-
minuit2/src/MnTiny.cxx | 9 +-
minuit2/src/MnUserFcn.cxx | 10 +-
minuit2/src/mnlsame.cxx | 88 +--
minuit2/src/mntplot.cxx | 316 +++++-----
22 files changed, 1505 insertions(+), 1462 deletions(-)
diff --git a/minuit2/inc/Minuit2/MnGlobalCorrelationCoeff.h b/minuit2/inc/Minuit2/MnGlobalCorrelationCoeff.h
index 82ac8c0ca6b..757d4dcb260 100644
--- a/minuit2/inc/Minuit2/MnGlobalCorrelationCoeff.h
+++ b/minuit2/inc/Minuit2/MnGlobalCorrelationCoeff.h
@@ -1,4 +1,4 @@
-// @(#)root/minuit2:$Name: $:$Id: MnGlobalCorrelationCoeff.h,v 1.2.2.3 2005/11/29 11:08:34 moneta Exp $
+// @(#)root/minuit2:$Name: $:$Id: MnGlobalCorrelationCoeff.h,v 1.1 2005/11/29 14:42:18 moneta Exp $
// Authors: M. Winkler, F. James, L. Moneta, A. Zsenei 2003-2005
/**********************************************************************
@@ -20,6 +20,9 @@ namespace ROOT {
namespace Minuit2 {
+/**
+ class for global correlation coefficient
+ */
class MnGlobalCorrelationCoeff {
public:
diff --git a/minuit2/inc/Minuit2/MnPosDef.h b/minuit2/inc/Minuit2/MnPosDef.h
index 7da31743233..cb35ca2d4c0 100644
--- a/minuit2/inc/Minuit2/MnPosDef.h
+++ b/minuit2/inc/Minuit2/MnPosDef.h
@@ -1,4 +1,4 @@
-// @(#)root/minuit2:$Name: $:$Id: MnPosDef.h,v 1.5.6.2 2005/11/29 11:08:34 moneta Exp $
+// @(#)root/minuit2:$Name: $:$Id: MnPosDef.h,v 1.1 2005/11/29 14:42:18 moneta Exp $
// Authors: M. Winkler, F. James, L. Moneta, A. Zsenei 2003-2005
/**********************************************************************
@@ -19,6 +19,10 @@ class MinimumState;
class MinimumError;
class MnMachinePrecision;
+/**
+ Force the covariance matrix to be positive defined
+ by adding extra terms in the diagonal
+ */
class MnPosDef {
public:
diff --git a/minuit2/inc/Minuit2/MnUserFcn.h b/minuit2/inc/Minuit2/MnUserFcn.h
index 078bfc38daa..82aa0ac8c68 100644
--- a/minuit2/inc/Minuit2/MnUserFcn.h
+++ b/minuit2/inc/Minuit2/MnUserFcn.h
@@ -1,4 +1,4 @@
-// @(#)root/minuit2:$Name: $:$Id: MnUserFcn.h,v 1.2.6.2 2005/11/29 11:08:34 moneta Exp $
+// @(#)root/minuit2:$Name: $:$Id: MnUserFcn.h,v 1.1 2005/11/29 14:42:18 moneta Exp $
// Authors: M. Winkler, F. James, L. Moneta, A. Zsenei 2003-2005
/**********************************************************************
@@ -19,6 +19,10 @@ namespace ROOT {
class MnUserTransformation;
+ /**
+ Wrapper used by Minuit of FCN interface
+ containing a reference to the transformation object
+ */
class MnUserFcn : public MnFcn {
public:
diff --git a/minuit2/src/MnFumiliMinimize.cxx b/minuit2/src/MnFumiliMinimize.cxx
index cc864306802..88ab9d88a59 100644
--- a/minuit2/src/MnFumiliMinimize.cxx
+++ b/minuit2/src/MnFumiliMinimize.cxx
@@ -1,4 +1,4 @@
-// @(#)root/minuit2:$Name: $:$Id: MnFumiliMinimize.cpp,v 1.1.4.4 2005/11/29 11:08:35 moneta Exp $
+// @(#)root/minuit2:$Name: $:$Id: MnFumiliMinimize.cxx,v 1.1 2005/11/29 14:43:31 moneta Exp $
// Authors: M. Winkler, F. James, L. Moneta, A. Zsenei 2003-2005
/**********************************************************************
@@ -17,20 +17,22 @@ namespace ROOT {
-// need to reimplement otherwise base class method is done
-FunctionMinimum MnFumiliMinimize::operator()(unsigned int maxfcn, double toler) {
- assert(fState.IsValid());
- unsigned int npar = VariableParameters();
-// assert(npar > 0);
- if(maxfcn == 0) maxfcn = 200 + 100*npar + 5*npar*npar;
- FunctionMinimum min = Minimizer().Minimize( Fcnbase(), fState, fStrategy, maxfcn, toler);
- fNumCall += min.NFcn();
- fState = min.UserState();
- return min;
+FunctionMinimum MnFumiliMinimize::operator()(unsigned int maxfcn, double toler) {
+ // minimize using Fumili
+ // need to reimplement otherwise base class method is done
+
+ assert(fState.IsValid());
+ unsigned int npar = VariableParameters();
+ // assert(npar > 0);
+ if(maxfcn == 0) maxfcn = 200 + 100*npar + 5*npar*npar;
+ FunctionMinimum min = Minimizer().Minimize( Fcnbase(), fState, fStrategy, maxfcn, toler);
+ fNumCall += min.NFcn();
+ fState = min.UserState();
+ return min;
}
- } // namespace Minuit2
+ } // namespace Minuit2
} // namespace ROOT
diff --git a/minuit2/src/MnFunctionCross.cxx b/minuit2/src/MnFunctionCross.cxx
index 4607bee25c3..3d1f0cbcdc3 100644
--- a/minuit2/src/MnFunctionCross.cxx
+++ b/minuit2/src/MnFunctionCross.cxx
@@ -1,4 +1,4 @@
-// @(#)root/minuit2:$Name: $:$Id: MnFunctionCross.cpp,v 1.9.2.4 2005/11/29 11:08:35 moneta Exp $
+// @(#)root/minuit2:$Name: $:$Id: MnFunctionCross.cxx,v 1.1 2005/11/29 14:43:31 moneta Exp $
// Authors: M. Winkler, F. James, L. Moneta, A. Zsenei 2003-2005
/**********************************************************************
@@ -25,306 +25,310 @@ namespace ROOT {
#define DEBUG
-MnCross MnFunctionCross::operator()(const std::vector<unsigned int>& par, const std::vector<double>& pmid, const std::vector<double>& pdir, double tlr, unsigned int maxcalls) const {
+MnCross MnFunctionCross::operator()(const std::vector<unsigned int>& par, const std::vector<double>& pmid,
+ const std::vector<double>& pdir, double tlr, unsigned int maxcalls) const {
+ // evaluate crossing point
+ // t.b.documented
+
// double edmmax = 0.5*0.001*toler*fFCN.Up();
- //std::cout<<"MnFunctionCross par "<<par.front()<<std::endl;
-
- unsigned int npar = par.size();
- unsigned int nfcn = 0;
- const MnMachinePrecision& prec = fState.Precision();
-// double tlr = 0.01;
- double tlf = tlr*fFCN.Up();
- double tla = tlr;
- unsigned int maxitr = 15;
- unsigned int ipt = 0;
- double aminsv = fFval;
- double aim = aminsv + fFCN.Up();
- //std::cout<<"aim= "<<aim<<std::endl;
- double aopt = 0.;
- bool limset = false;
- std::vector<double> alsb(3, 0.), flsb(3, 0.);
- double up = fFCN.Up();
-
- double aulim = 100.;
- for(unsigned int i = 0; i < par.size(); i++) {
- unsigned int kex = par[i];
- if(fState.Parameter(kex).HasLimits()) {
- double zmid = pmid[i];
- double zdir = pdir[i];
- if(fabs(zdir) < fState.Precision().Eps()) continue;
-// double zlim = 0.;
- if(zdir > 0. && fState.Parameter(kex).HasUpperLimit()) {
- double zlim = fState.Parameter(kex).UpperLimit();
- aulim = std::min(aulim, (zlim-zmid)/zdir);
- }
- else if(zdir < 0. && fState.Parameter(kex).HasLowerLimit()) {
- double zlim = fState.Parameter(kex).LowerLimit();
- aulim = std::min(aulim, (zlim-zmid)/zdir);
+ //std::cout<<"MnFunctionCross par "<<par.front()<<std::endl;
+
+ unsigned int npar = par.size();
+ unsigned int nfcn = 0;
+ const MnMachinePrecision& prec = fState.Precision();
+ // double tlr = 0.01;
+ double tlf = tlr*fFCN.Up();
+ double tla = tlr;
+ unsigned int maxitr = 15;
+ unsigned int ipt = 0;
+ double aminsv = fFval;
+ double aim = aminsv + fFCN.Up();
+ //std::cout<<"aim= "<<aim<<std::endl;
+ double aopt = 0.;
+ bool limset = false;
+ std::vector<double> alsb(3, 0.), flsb(3, 0.);
+ double up = fFCN.Up();
+
+ double aulim = 100.;
+ for(unsigned int i = 0; i < par.size(); i++) {
+ unsigned int kex = par[i];
+ if(fState.Parameter(kex).HasLimits()) {
+ double zmid = pmid[i];
+ double zdir = pdir[i];
+ if(fabs(zdir) < fState.Precision().Eps()) continue;
+ // double zlim = 0.;
+ if(zdir > 0. && fState.Parameter(kex).HasUpperLimit()) {
+ double zlim = fState.Parameter(kex).UpperLimit();
+ aulim = std::min(aulim, (zlim-zmid)/zdir);
+ }
+ else if(zdir < 0. && fState.Parameter(kex).HasLowerLimit()) {
+ double zlim = fState.Parameter(kex).LowerLimit();
+ aulim = std::min(aulim, (zlim-zmid)/zdir);
+ }
}
- }
- }
- //std::cout<<"100"<<std::endl;
-
- if(aulim < aopt+tla) limset = true;
-
- MnMigrad migrad(fFCN, fState, MnStrategy(std::max(0, int(fStrategy.Strategy()-1))));
-
- for(unsigned int i = 0; i < npar; i++) {
- migrad.SetValue(par[i], pmid[i]);
- }
-
- FunctionMinimum min0 = migrad(maxcalls, tlr);
- nfcn += min0.NFcn();
-
- if(min0.HasReachedCallLimit())
- return MnCross(min0.UserState(), nfcn, MnCross::CrossFcnLimit());
- if(!min0.IsValid()) return MnCross(fState, nfcn);
- if(limset == true && min0.Fval() < aim)
- return MnCross(min0.UserState(), nfcn, MnCross::CrossParLimit());
-
- ipt++;
- alsb[0] = 0.;
- flsb[0] = min0.Fval();
- flsb[0] = std::max(flsb[0], aminsv + 0.1*up);
- aopt = sqrt(up/(flsb[0]-aminsv)) - 1.;
- if(fabs(flsb[0] - aim) < tlf) return MnCross(aopt, min0.UserState(), nfcn);
-
- if(aopt > 1.) aopt = 1.;
- if(aopt < -0.5) aopt = -0.5;
- limset = false;
- if(aopt > aulim) {
- aopt = aulim;
- limset = true;
- }
-
- for(unsigned int i = 0; i < npar; i++) {
- migrad.SetValue(par[i], pmid[i] + (aopt)*pdir[i]);
- }
-
- FunctionMinimum min1 = migrad(maxcalls, tlr);
- nfcn += min1.NFcn();
-
- if(min1.HasReachedCallLimit())
- return MnCross(min1.UserState(), nfcn, MnCross::CrossFcnLimit());
- if(!min1.IsValid()) return MnCross(fState, nfcn);
- if(limset == true && min1.Fval() < aim)
- return MnCross(min1.UserState(), nfcn, MnCross::CrossParLimit());
-
- ipt++;
- alsb[1] = aopt;
- flsb[1] = min1.Fval();
- double dfda = (flsb[1] - flsb[0])/(alsb[1] - alsb[0]);
-// if(dfda > 0.) goto L460;
-
+ }
+ //std::cout<<"100"<<std::endl;
+
+ if(aulim < aopt+tla) limset = true;
+
+ MnMigrad migrad(fFCN, fState, MnStrategy(std::max(0, int(fStrategy.Strategy()-1))));
+
+ for(unsigned int i = 0; i < npar; i++) {
+ migrad.SetValue(par[i], pmid[i]);
+ }
+
+ FunctionMinimum min0 = migrad(maxcalls, tlr);
+ nfcn += min0.NFcn();
+
+ if(min0.HasReachedCallLimit())
+ return MnCross(min0.UserState(), nfcn, MnCross::CrossFcnLimit());
+ if(!min0.IsValid()) return MnCross(fState, nfcn);
+ if(limset == true && min0.Fval() < aim)
+ return MnCross(min0.UserState(), nfcn, MnCross::CrossParLimit());
+
+ ipt++;
+ alsb[0] = 0.;
+ flsb[0] = min0.Fval();
+ flsb[0] = std::max(flsb[0], aminsv + 0.1*up);
+ aopt = sqrt(up/(flsb[0]-aminsv)) - 1.;
+ if(fabs(flsb[0] - aim) < tlf) return MnCross(aopt, min0.UserState(), nfcn);
+
+ if(aopt > 1.) aopt = 1.;
+ if(aopt < -0.5) aopt = -0.5;
+ limset = false;
+ if(aopt > aulim) {
+ aopt = aulim;
+ limset = true;
+ }
+
+ for(unsigned int i = 0; i < npar; i++) {
+ migrad.SetValue(par[i], pmid[i] + (aopt)*pdir[i]);
+ }
+
+ FunctionMinimum min1 = migrad(maxcalls, tlr);
+ nfcn += min1.NFcn();
+
+ if(min1.HasReachedCallLimit())
+ return MnCross(min1.UserState(), nfcn, MnCross::CrossFcnLimit());
+ if(!min1.IsValid()) return MnCross(fState, nfcn);
+ if(limset == true && min1.Fval() < aim)
+ return MnCross(min1.UserState(), nfcn, MnCross::CrossParLimit());
+
+ ipt++;
+ alsb[1] = aopt;
+ flsb[1] = min1.Fval();
+ double dfda = (flsb[1] - flsb[0])/(alsb[1] - alsb[0]);
+ // if(dfda > 0.) goto L460;
+
L300:
- if(dfda < 0.) {
- //std::cout<<"300"<<std::endl;
- unsigned int maxlk = maxitr - ipt;
- for(unsigned int it = 0; it < maxlk; it++) {
- alsb[0] = alsb[1];
- flsb[0] = flsb[1];
- // LM: Add + 1, looking at Fortran code it starts from 1 ( see bug #8396)
- aopt = alsb[0] + 0.2*(it+1);
- limset = false;
- if(aopt > aulim) {
- aopt = aulim;
- limset = true;
- }
- for(unsigned int i = 0; i < npar; i++) {
- migrad.SetValue(par[i], pmid[i] + (aopt)*pdir[i]);
- }
- min1 = migrad(maxcalls, tlr);
- nfcn += min1.NFcn();
-
- if(min1.HasReachedCallLimit())
- return MnCross(min1.UserState(), nfcn, MnCross::CrossFcnLimit());
- if(!min1.IsValid()) return MnCross(fState, nfcn);
- if(limset == true && min1.Fval() < aim)
- return MnCross(min1.UserState(), nfcn, MnCross::CrossParLimit());
- ipt++;
- alsb[1] = aopt;
- flsb[1] = min1.Fval();
- dfda = (flsb[1] - flsb[0])/(alsb[1] - alsb[0]);
-// if(dfda > 0.) goto L460;
- if(dfda > 0.) break;
- }
- if(ipt > maxitr) return MnCross(fState, nfcn);
- } //if(dfda < 0.)
-
+ if(dfda < 0.) {
+ //std::cout<<"300"<<std::endl;
+ unsigned int maxlk = maxitr - ipt;
+ for(unsigned int it = 0; it < maxlk; it++) {
+ alsb[0] = alsb[1];
+ flsb[0] = flsb[1];
+ // LM: Add + 1, looking at Fortran code it starts from 1 ( see bug #8396)
+ aopt = alsb[0] + 0.2*(it+1);
+ limset = false;
+ if(aopt > aulim) {
+ aopt = aulim;
+ limset = true;
+ }
+ for(unsigned int i = 0; i < npar; i++) {
+ migrad.SetValue(par[i], pmid[i] + (aopt)*pdir[i]);
+ }
+ min1 = migrad(maxcalls, tlr);
+ nfcn += min1.NFcn();
+
+ if(min1.HasReachedCallLimit())
+ return MnCross(min1.UserState(), nfcn, MnCross::CrossFcnLimit());
+ if(!min1.IsValid()) return MnCross(fState, nfcn);
+ if(limset == true && min1.Fval() < aim)
+ return MnCross(min1.UserState(), nfcn, MnCross::CrossParLimit());
+ ipt++;
+ alsb[1] = aopt;
+ flsb[1] = min1.Fval();
+ dfda = (flsb[1] - flsb[0])/(alsb[1] - alsb[0]);
+ // if(dfda > 0.) goto L460;
+ if(dfda > 0.) break;
+ }
+ if(ipt > maxitr) return MnCross(fState, nfcn);
+ } //if(dfda < 0.)
+
L460:
- //std::cout<<"460"<<std::endl;
-
- aopt = alsb[1] + (aim-flsb[1])/dfda;
- double fdist = std::min(fabs(aim - flsb[0]), fabs(aim - flsb[1]));
- double adist = std::min(fabs(aopt - alsb[0]), fabs(aopt - alsb[1]));
- tla = tlr;
- if(fabs(aopt) > 1.) tla = tlr*fabs(aopt);
- if(adist < tla && fdist < tlf) return MnCross(aopt, min1.UserState(), nfcn);
- if(ipt > maxitr) return MnCross(fState, nfcn);
- double bmin = std::min(alsb[0], alsb[1]) - 1.;
- if(aopt < bmin) aopt = bmin;
- double bmax = std::max(alsb[0], alsb[1]) + 1.;
- if(aopt > bmax) aopt = bmax;
-
- limset = false;
- if(aopt > aulim) {
- aopt = aulim;
- limset = true;
- }
-
- for(unsigned int i = 0; i < npar; i++) {
- migrad.SetValue(par[i], pmid[i] + (aopt)*pdir[i]);
- }
- FunctionMinimum min2 = migrad(maxcalls, tlr);
- nfcn += min2.NFcn();
-
- if(min2.HasReachedCallLimit())
- return MnCross(min2.UserState(), nfcn, MnCross::CrossFcnLimit());
- if(!min2.IsValid()) return MnCross(fState, nfcn);
- if(limset == true && min2.Fval() < aim)
- return MnCross(min2.UserState(), nfcn, MnCross::CrossParLimit());
-
- ipt++;
- alsb[2] = aopt;
- flsb[2] = min2.Fval();
-
- double ecarmn = fabs(flsb[2] - aim);
- double ecarmx = 0.;
- unsigned int ibest = 2;
- unsigned int iworst = 0;
- unsigned int noless = 0;
-
- for(unsigned int i = 0; i < 3; i++) {
- double ecart = fabs(flsb[i] - aim);
- if(ecart > ecarmx) {
- ecarmx = ecart;
- iworst = i;
- }
- if(ecart < ecarmn) {
- ecarmn = ecart;
- ibest = i;
- }
- if(flsb[i] < aim) noless++;
- }
-
- //std::cout<<"480"<<std::endl;
-
- if(noless == 1 || noless == 2) goto L500;
- if(noless == 0 && ibest != 2) return MnCross(fState, nfcn);
- if(noless == 3 && ibest != 2) {
- alsb[1] = alsb[2];
- flsb[1] = flsb[2];
- goto L300;
- }
-
- flsb[iworst] = flsb[2];
- alsb[iworst] = alsb[2];
- dfda = (flsb[1] - flsb[0])/(alsb[1] - alsb[0]);
- goto L460;
-
-L500:
- //std::cout<<"500"<<std::endl;
- do {
- MnParabola parbol = MnParabolaFactory()(MnParabolaPoint(alsb[0], flsb[0]), MnParabolaPoint(alsb[1], flsb[1]), MnParabolaPoint(alsb[2], flsb[2]));
- // aopt = parbol.X_pos(aim);
- //std::cout<<"alsb1,2,3= "<<alsb[0]<<", "<<alsb[1]<<", "<<alsb[2]<<std::endl;
- //std::cout<<"flsb1,2,3= "<<flsb[0]<<", "<<flsb[1]<<", "<<flsb[2]<<std::endl;
-
- double coeff1 = parbol.c();
- double coeff2 = parbol.b();
- double coeff3 = parbol.a();
- double determ = coeff2*coeff2 - 4.*coeff3*(coeff1 - aim);
- if(determ < prec.Eps()) return MnCross(fState, nfcn);
- double rt = sqrt(determ);
- double x1 = (-coeff2 + rt)/(2.*coeff3);
- double x2 = (-coeff2 - rt)/(2.*coeff3);
- double s1 = coeff2 + 2.*x1*coeff3;
- double s2 = coeff2 + 2.*x2*coeff3;
-
- if(s1*s2 > 0.) std::cout<<"MnFunctionCross problem 1"<<std::endl;
- aopt = x1;
- double slope = s1;
- if(s2 > 0.) {
- aopt = x2;
- slope = s2;
- }
-
- tla = tlr;
- if(fabs(aopt) > 1.) tla = tlr*fabs(aopt);
- if(fabs(aopt - alsb[ibest]) < tla && fabs(flsb[ibest] - aim) < tlf)
- return MnCross(aopt, min2.UserState(), nfcn);
-
-// if(ipt > maxitr) return MnCross();
-
- unsigned int ileft = 3;
- unsigned int iright = 3;
- unsigned int iout = 3;
- ibest = 0;
- ecarmx = 0.;
- ecarmn = fabs(aim-flsb[0]);
- for(unsigned int i = 0; i < 3; i++) {
- double ecart = fabs(flsb[i] - aim);
- if(ecart < ecarmn) {
- ecarmn = ecart;
- ibest = i;
- }
- if(ecart > ecarmx) ecarmx = ecart;
- if(flsb[i] > aim) {
- if(iright == 3) iright = i;
- else if(flsb[i] > flsb[iright]) iout = i;
- else {
- iout = iright;
- iright = i;
- }
- } else if(ileft == 3) ileft = i;
- else if(flsb[i] < flsb[ileft]) iout = i;
- else {
- iout = ileft;
- ileft = i;
- }
- }
- //std::cout<<"550"<<std::endl;
-
- if(ecarmx > 10.*fabs(flsb[iout] - aim))
- aopt = 0.5*(aopt + 0.5*(alsb[iright] + alsb[ileft]));
- double smalla = 0.1*tla;
- if(slope*smalla > tlf) smalla = tlf/slope;
- double aleft = alsb[ileft] + smalla;
- double aright = alsb[iright] - smalla;
- if(aopt < aleft) aopt = aleft;
- if(aopt > aright) aopt = aright;
- if(aleft > aright) aopt = 0.5*(aleft + aright);
-
- limset = false;
- if(aopt > aulim) {
+ //std::cout<<"460"<<std::endl;
+
+ aopt = alsb[1] + (aim-flsb[1])/dfda;
+ double fdist = std::min(fabs(aim - flsb[0]), fabs(aim - flsb[1]));
+ double adist = std::min(fabs(aopt - alsb[0]), fabs(aopt - alsb[1]));
+ tla = tlr;
+ if(fabs(aopt) > 1.) tla = tlr*fabs(aopt);
+ if(adist < tla && fdist < tlf) return MnCross(aopt, min1.UserState(), nfcn);
+ if(ipt > maxitr) return MnCross(fState, nfcn);
+ double bmin = std::min(alsb[0], alsb[1]) - 1.;
+ if(aopt < bmin) aopt = bmin;
+ double bmax = std::max(alsb[0], alsb[1]) + 1.;
+ if(aopt > bmax) aopt = bmax;
+
+ limset = false;
+ if(aopt > aulim) {
aopt = aulim;
limset = true;
- }
-
- for(unsigned int i = 0; i < npar; i++) {
+ }
+
+ for(unsigned int i = 0; i < npar; i++) {
migrad.SetValue(par[i], pmid[i] + (aopt)*pdir[i]);
- }
- min2 = migrad(maxcalls, tlr);
- nfcn += min2.NFcn();
-
- if(min2.HasReachedCallLimit())
+ }
+ FunctionMinimum min2 = migrad(maxcalls, tlr);
+ nfcn += min2.NFcn();
+
+ if(min2.HasReachedCallLimit())
return MnCross(min2.UserState(), nfcn, MnCross::CrossFcnLimit());
- if(!min2.IsValid()) return MnCross(nfcn);
- if(limset == true && min2.Fval() < aim)
+ if(!min2.IsValid()) return MnCross(fState, nfcn);
+ if(limset == true && min2.Fval() < aim)
return MnCross(min2.UserState(), nfcn, MnCross::CrossParLimit());
-
- ipt++;
- alsb[iout] = aopt;
- flsb[iout] = min2.Fval();
- ibest = iout;
- } while(ipt < maxitr);
-
- // goto L500;
-
- return MnCross(fState, nfcn);
+
+ ipt++;
+ alsb[2] = aopt;
+ flsb[2] = min2.Fval();
+
+ double ecarmn = fabs(flsb[2] - aim);
+ double ecarmx = 0.;
+ unsigned int ibest = 2;
+ unsigned int iworst = 0;
+ unsigned int noless = 0;
+
+ for(unsigned int i = 0; i < 3; i++) {
+ double ecart = fabs(flsb[i] - aim);
+ if(ecart > ecarmx) {
+ ecarmx = ecart;
+ iworst = i;
+ }
+ if(ecart < ecarmn) {
+ ecarmn = ecart;
+ ibest = i;
+ }
+ if(flsb[i] < aim) noless++;
+ }
+
+ //std::cout<<"480"<<std::endl;
+
+ if(noless == 1 || noless == 2) goto L500;
+ if(noless == 0 && ibest != 2) return MnCross(fState, nfcn);
+ if(noless == 3 && ibest != 2) {
+ alsb[1] = alsb[2];
+ flsb[1] = flsb[2];
+ goto L300;
+ }
+
+ flsb[iworst] = flsb[2];
+ alsb[iworst] = alsb[2];
+ dfda = (flsb[1] - flsb[0])/(alsb[1] - alsb[0]);
+ goto L460;
+
+L500:
+ //std::cout<<"500"<<std::endl;
+ do {
+ MnParabola parbol = MnParabolaFactory()(MnParabolaPoint(alsb[0], flsb[0]), MnParabolaPoint(alsb[1], flsb[1]), MnParabolaPoint(alsb[2], flsb[2]));
+ // aopt = parbol.X_pos(aim);
+ //std::cout<<"alsb1,2,3= "<<alsb[0]<<", "<<alsb[1]<<", "<<alsb[2]<<std::endl;
+ //std::cout<<"flsb1,2,3= "<<flsb[0]<<", "<<flsb[1]<<", "<<flsb[2]<<std::endl;
+
+ double coeff1 = parbol.c();
+ double coeff2 = parbol.b();
+ double coeff3 = parbol.a();
+ double determ = coeff2*coeff2 - 4.*coeff3*(coeff1 - aim);
+ if(determ < prec.Eps()) return MnCross(fState, nfcn);
+ double rt = sqrt(determ);
+ double x1 = (-coeff2 + rt)/(2.*coeff3);
+ double x2 = (-coeff2 - rt)/(2.*coeff3);
+ double s1 = coeff2 + 2.*x1*coeff3;
+ double s2 = coeff2 + 2.*x2*coeff3;
+
+ if(s1*s2 > 0.) std::cout<<"MnFunctionCross problem 1"<<std::endl;
+ aopt = x1;
+ double slope = s1;
+ if(s2 > 0.) {
+ aopt = x2;
+ slope = s2;
+ }
+
+ tla = tlr;
+ if(fabs(aopt) > 1.) tla = tlr*fabs(aopt);
+ if(fabs(aopt - alsb[ibest]) < tla && fabs(flsb[ibest] - aim) < tlf)
+ return MnCross(aopt, min2.UserState(), nfcn);
+
+ // if(ipt > maxitr) return MnCross();
+
+ unsigned int ileft = 3;
+ unsigned int iright = 3;
+ unsigned int iout = 3;
+ ibest = 0;
+ ecarmx = 0.;
+ ecarmn = fabs(aim-flsb[0]);
+ for(unsigned int i = 0; i < 3; i++) {
+ double ecart = fabs(flsb[i] - aim);
+ if(ecart < ecarmn) {
+ ecarmn = ecart;
+ ibest = i;
+ }
+ if(ecart > ecarmx) ecarmx = ecart;
+ if(flsb[i] > aim) {
+ if(iright == 3) iright = i;
+ else if(flsb[i] > flsb[iright]) iout = i;
+ else {
+ iout = iright;
+ iright = i;
+ }
+ } else if(ileft == 3) ileft = i;
+ else if(flsb[i] < flsb[ileft]) iout = i;
+ else {
+ iout = ileft;
+ ileft = i;
+ }
+ }
+ //std::cout<<"550"<<std::endl;
+
+ if(ecarmx > 10.*fabs(flsb[iout] - aim))
+ aopt = 0.5*(aopt + 0.5*(alsb[iright] + alsb[ileft]));
+ double smalla = 0.1*tla;
+ if(slope*smalla > tlf) smalla = tlf/slope;
+ double aleft = alsb[ileft] + smalla;
+ double aright = alsb[iright] - smalla;
+ if(aopt < aleft) aopt = aleft;
+ if(aopt > aright) aopt = aright;
+ if(aleft > aright) aopt = 0.5*(aleft + aright);
+
+ limset = false;
+ if(aopt > aulim) {
+ aopt = aulim;
+ limset = true;
+ }
+
+ for(unsigned int i = 0; i < npar; i++) {
+ migrad.SetValue(par[i], pmid[i] + (aopt)*pdir[i]);
+ }
+ min2 = migrad(maxcalls, tlr);
+ nfcn += min2.NFcn();
+
+ if(min2.HasReachedCallLimit())
+ return MnCross(min2.UserState(), nfcn, MnCross::CrossFcnLimit());
+ if(!min2.IsValid()) return MnCross(nfcn);
+ if(limset == true && min2.Fval() < aim)
+ return MnCross(min2.UserState(), nfcn, MnCross::CrossParLimit());
+
+ ipt++;
+ alsb[iout] = aopt;
+ flsb[iout] = min2.Fval();
+ ibest = iout;
+ } while(ipt < maxitr);
+
+ // goto L500;
+
+ return MnCross(fState, nfcn);
}
- } // namespace Minuit2
+ } // namespace Minuit2
} // namespace ROOT
diff --git a/minuit2/src/MnGlobalCorrelationCoeff.cxx b/minuit2/src/MnGlobalCorrelationCoeff.cxx
index 3a47298cfee..156b83323d4 100644
--- a/minuit2/src/MnGlobalCorrelationCoeff.cxx
+++ b/minuit2/src/MnGlobalCorrelationCoeff.cxx
@@ -1,4 +1,4 @@
-// @(#)root/minuit2:$Name: $:$Id: MnGlobalCorrelationCoeff.cxx,v 1.1 2005/11/29 14:43:31 moneta Exp $
+// @(#)root/minuit2:$Name: $:$Id: MnGlobalCorrelationCoeff.cxx,v 1.2 2006/04/13 08:39:23 moneta Exp $
// Authors: M. Winkler, F. James, L. Moneta, A. Zsenei 2003-2005
/**********************************************************************
@@ -17,24 +17,25 @@ namespace ROOT {
MnGlobalCorrelationCoeff::MnGlobalCorrelationCoeff(const MnAlgebraicSymMatrix& cov) : fGlobalCC(std::vector<double>()), fValid(true) {
-
- MnAlgebraicSymMatrix inv(cov);
- int ifail = Invert(inv);
- if(ifail != 0) {
+ // constructor: calculate global correlation given a symmetric matrix
+
+ MnAlgebraicSymMatrix inv(cov);
+ int ifail = Invert(inv);
+ if(ifail != 0) {
#ifdef WARNINGMSG
- std::cout<<"MnGlobalCorrelationCoeff: inversion of matrix fails."<<std::endl;
+ std::cout<<"MnGlobalCorrelationCoeff: inversion of matrix fails."<<std::endl;
#endif
- fValid = false;
- } else {
-
- for(unsigned int i = 0; i < cov.Nrow(); i++) {
- double denom = inv(i,i)*cov(i,i);
- if(denom < 1. && denom > 0.) fGlobalCC.push_back(0.);
- else fGlobalCC.push_back(std::sqrt(1. - 1./denom));
- }
- }
+ fValid = false;
+ } else {
+
+ for(unsigned int i = 0; i < cov.Nrow(); i++) {
+ double denom = inv(i,i)*cov(i,i);
+ if(denom < 1. && denom > 0.) fGlobalCC.push_back(0.);
+ else fGlobalCC.push_back(std::sqrt(1. - 1./denom));
+ }
+ }
}
- } // namespace Minuit2
+ } // namespace Minuit2
} // namespace ROOT
diff --git a/minuit2/src/MnHesse.cxx b/minuit2/src/MnHesse.cxx
index 521c7293123..dcb09745b39 100644
--- a/minuit2/src/MnHesse.cxx
+++ b/minuit2/src/MnHesse.cxx
@@ -1,4 +1,4 @@
-// @(#)root/minuit2:$Name: $:$Id: MnHesse.cxx,v 1.2 2006/01/25 12:20:49 moneta Exp $
+// @(#)root/minuit2:$Name: $:$Id: MnHesse.cxx,v 1.3 2006/04/12 16:30:30 moneta Exp $
// Authors: M. Winkler, F. James, L. Moneta, A. Zsenei 2003-2005
/**********************************************************************
@@ -24,301 +24,309 @@ namespace ROOT {
namespace Minuit2 {
-MnUserParameterState MnHesse::operator()(const FCNBase& fcn, const std::vector<double>& par, const std::vector<double>& err, unsigned int maxcalls) const {
- return (*this)(fcn, MnUserParameterState(par, err), maxcalls);
+MnUserParameterState MnHesse::operator()(const FCNBase& fcn, const std::vector<double>& par, const std::vector<double>& err, unsigned int maxcalls) const {
+ // interface from vector of params and errors
+ return (*this)(fcn, MnUserParameterState(par, err), maxcalls);
}
MnUserParameterState MnHesse::operator()(const FCNBase& fcn, const std::vector<double>& par, const std::vector<double>& cov, unsigned int nrow, unsigned int maxcalls) const {
- return (*this)(fcn, MnUserParameterState(par, cov, nrow), maxcalls);
+ // interface from vector of params and covariance
+ return (*this)(fcn, MnUserParameterState(par, cov, nrow), maxcalls);
}
MnUserParameterState MnHesse::operator()(const FCNBase& fcn, const std::vector<double>& par, const MnUserCovariance& cov, unsigned int maxcalls) const {
- return (*this)(fcn, MnUserParameterState(par, cov), maxcalls);
+ // interface from vector of params and covariance
+ return (*this)(fcn, MnUserParameterState(par, cov), maxcalls);
}
MnUserParameterState MnHesse::operator()(const FCNBase& fcn, const MnUserParameters& par, unsigned int maxcalls) const {
- return (*this)(fcn, MnUserParameterState(par), maxcalls);
+ // interface from MnUserParameters
+ return (*this)(fcn, MnUserParameterState(par), maxcalls);
}
MnUserParameterState MnHesse::operator()(const FCNBase& fcn, const MnUserParameters& par, const MnUserCovariance& cov, unsigned int maxcalls) const {
- return (*this)(fcn, MnUserParameterState(par, cov), maxcalls);
+ // interface from MnUserParameters and MnUserCovariance
+ return (*this)(fcn, MnUserParameterState(par, cov), maxcalls);
}
MnUserParameterState MnHesse::operator()(const FCNBase& fcn, const MnUserParameterState& state, unsigned int maxcalls) const {
-
- unsigned int n = state.VariableParameters();
- MnUserFcn mfcn(fcn, state.Trafo());
- MnAlgebraicVector x(n);
- for(unsigned int i = 0; i < n; i++) x(i) = state.IntParameters()[i];
- double amin = mfcn(x);
- Numerical2PGradientCalculator gc(mfcn, state.Trafo(), fStrategy);
- MinimumParameters par(x, amin);
- FunctionGradient gra = gc(par);
- MinimumState tmp = (*this)(mfcn, MinimumState(par, MinimumError(MnAlgebraicSymMatrix(n), 1.), gra, state.Edm(), state.NFcn()), state.Trafo(), maxcalls);
-
- return MnUserParameterState(tmp, fcn.Up(), state.Trafo());
+ // interface from MnUserParameterState
+ // create Minimum state and use that interface
+ unsigned int n = state.VariableParameters();
+ MnUserFcn mfcn(fcn, state.Trafo());
+ MnAlgebraicVector x(n);
+ for(unsigned int i = 0; i < n; i++) x(i) = state.IntParameters()[i];
+ double amin = mfcn(x);
+ Numerical2PGradientCalculator gc(mfcn, state.Trafo(), fStrategy);
+ MinimumParameters par(x, amin);
+ FunctionGradient gra = gc(par);
+ MinimumState tmp = (*this)(mfcn, MinimumState(par, MinimumError(MnAlgebraicSymMatrix(n), 1.), gra, state.Edm(), state.NFcn()), state.Trafo(), maxcalls);
+
+ return MnUserParameterState(tmp, fcn.Up(), state.Trafo());
}
MinimumState MnHesse::operator()(const MnFcn& mfcn, const MinimumState& st, const MnUserTransformation& trafo, unsigned int maxcalls) const {
-
- const MnMachinePrecision& prec = trafo.Precision();
- // make sure starting at the right place
- double amin = mfcn(st.Vec());
- double aimsag = sqrt(prec.Eps2())*(fabs(amin)+mfcn.Up());
-
- // diagonal Elements first
-
- unsigned int n = st.Parameters().Vec().size();
- if(maxcalls == 0) maxcalls = 200 + 100*n + 5*n*n;
-
- MnAlgebraicSymMatrix vhmat(n);
- MnAlgebraicVector g2 = st.Gradient().G2();
- MnAlgebraicVector gst = st.Gradient().Gstep();
- MnAlgebraicVector grd = st.Gradient().Grad();
- MnAlgebraicVector dirin = st.Gradient().Gstep();
- MnAlgebraicVector yy(n);
- if(st.Gradient().IsAnalytical()) {
- InitialGradientCalculator igc(mfcn, trafo, fStrategy);
- FunctionGradient tmp = igc(st.Parameters());
- gst = tmp.Gstep();
- dirin = tmp.Gstep();
- g2 = tmp.G2();
- }
-
- MnAlgebraicVector x = st.Parameters().Vec();
-
- for(unsigned int i = 0; i < n; i++) {
-
- double xtf = x(i);
- double dmin = 8.*prec.Eps2()*(fabs(xtf) + prec.Eps2());
- double d = fabs(gst(i));
- if(d < dmin) d = dmin;
-
- for(unsigned int icyc = 0; icyc < Ncycles(); icyc++) {
- double sag = 0.;
- double fs1 = 0.;
- double fs2 = 0.;
- for(unsigned int multpy = 0; multpy < 5; multpy++) {
- x(i) = xtf + d;
- fs1 = mfcn(x);
- x(i) = xtf - d;
- fs2 = mfcn(x);
- x(i) = xtf;
- sag = 0.5*(fs1+fs2-2.*amin);
- if(fabs(sag) > prec.Eps2()) goto L30; // break;
- if(trafo.Parameter(i).HasLimits()) {
- if(d > 0.5) goto L26;
- d *= 10.;
- if(d > 0.5) d = 0.51;
- continue;
- }
- d *= 10.;
- }
+ // interface from MinimumState and MnUserTransformation
+ // Function who does the real Hessian calculations
+
+ const MnMachinePrecision& prec = trafo.Precision();
+ // make sure starting at the right place
+ double amin = mfcn(st.Vec());
+ double aimsag = sqrt(prec.Eps2())*(fabs(amin)+mfcn.Up());
+
+ // diagonal Elements first
+
+ unsigned int n = st.Parameters().Vec().size();
+ if(maxcalls == 0) maxcalls = 200 + 100*n + 5*n*n;
+
+ MnAlgebraicSymMatrix vhmat(n);
+ MnAlgebraicVector g2 = st.Gradient().G2();
+ MnAlgebraicVector gst = st.Gradient().Gstep();
+ MnAlgebraicVector grd = st.Gradient().Grad();
+ MnAlgebraicVector dirin = st.Gradient().Gstep();
+ MnAlgebraicVector yy(n);
+ if(st.Gradient().IsAnalytical()) {
+ InitialGradientCalculator igc(mfcn, trafo, fStrategy);
+ FunctionGradient tmp = igc(st.Parameters());
+ gst = tmp.Gstep();
+ dirin = tmp.Gstep();
+ g2 = tmp.G2();
+ }
+
+ MnAlgebraicVector x = st.Parameters().Vec();
+
+ for(unsigned int i = 0; i < n; i++) {
+
+ double xtf = x(i);
+ double dmin = 8.*prec.Eps2()*(fabs(xtf) + prec.Eps2());
+ double d = fabs(gst(i));
+ if(d < dmin) d = dmin;
+ for(unsigned int icyc = 0; icyc < Ncycles(); icyc++) {
+ double sag = 0.;
+ double fs1 = 0.;
+ double fs2 = 0.;
+ for(unsigned int multpy = 0; multpy < 5; multpy++) {
+ x(i) = xtf + d;
+ fs1 = mfcn(x);
+ x(i) = xtf - d;
+ fs2 = mfcn(x);
+ x(i) = xtf;
+ sag = 0.5*(fs1+fs2-2.*amin);
+ if(fabs(sag) > prec.Eps2()) goto L30; // break;
+ if(trafo.Parameter(i).HasLimits()) {
+ if(d > 0.5) goto L26;
+ d *= 10.;
+ if(d > 0.5) d = 0.51;
+ continue;
+ }
+ d *= 10.;
+ }
+
L26:
#ifdef WARNINGMSG
- std::cout<<"MnHesse: 2nd derivative zero for Parameter "<<i<<std::endl;
- std::cout<<"MnHesse fails and will return diagonal matrix "<<std::endl;
+ std::cout<<"MnHesse: 2nd derivative zero for Parameter "<<i<<std::endl;
+ std::cout<<"MnHesse fails and will return diagonal matrix "<<std::endl;
#endif
-
- for(unsigned int j = 0; j < n; j++) {
- double tmp = g2(j) < prec.Eps2() ? 1. : 1./g2(j);
- vhmat(j,j) = tmp < prec.Eps2() ? 1. : tmp;
- }
-
- return MinimumState(st.Parameters(), MinimumError(vhmat, MinimumError::MnHesseFailed()), st.Gradient(), st.Edm(), mfcn.NumOfCalls());
-
+
+ for(unsigned int j = 0; j < n; j++) {
+ double tmp = g2(j) < prec.Eps2() ? 1. : 1./g2(j);
+ vhmat(j,j) = tmp < prec.Eps2() ? 1. : tmp;
+ }
+
+ return MinimumState(st.Parameters(), MinimumError(vhmat, MinimumError::MnHesseFailed()), st.Gradient(), st.Edm(), mfcn.NumOfCalls());
+
L30:
- double g2bfor = g2(i);
- g2(i) = 2.*sag/(d*d);
- grd(i) = (fs1-fs2)/(2.*d);
- gst(i) = d;
- dirin(i) = d;
- yy(i) = fs1;
- double dlast = d;
- d = sqrt(2.*aimsag/fabs(g2(i)));
- if(trafo.Parameter(i).HasLimits()) d = std::min(0.5, d);
- if(d < dmin) d = dmin;
-
- // see if converged
- if(fabs((d-dlast)/d) < Tolerstp()) break;
- if(fabs((g2(i)-g2bfor)/g2(i)) < TolerG2()) break;
- d = std::min(d, 10.*dlast);
- d = std::max(d, 0.1*dlast);
- }
- vhmat(i,i) = g2(i);
- if(mfcn.NumOfCalls() > maxcalls) {
- //std::cout<<"maxcalls " << maxcalls << " " << mfcn.NumOfCalls() << " " << st.NFcn() << std::endl;
-
+ double g2bfor = g2(i);
+ g2(i) = 2.*sag/(d*d);
+ grd(i) = (fs1-fs2)/(2.*d);
+ gst(i) = d;
+ dirin(i) = d;
+ yy(i) = fs1;
+ double dlast = d;
+ d = sqrt(2.*aimsag/fabs(g2(i)));
+ if(trafo.Parameter(i).HasLimits()) d = std::min(0.5, d);
+ if(d < dmin) d = dmin;
+
+ // see if converged
+ if(fabs((d-dlast)/d) < Tolerstp()) break;
+ if(fabs((g2(i)-g2bfor)/g2(i)) < TolerG2()) break;
+ d = std::min(d, 10.*dlast);
+ d = std::max(d, 0.1*dlast);
+ }
+ vhmat(i,i) = g2(i);
+ if(mfcn.NumOfCalls() > maxcalls) {
+ //std::cout<<"maxcalls " << maxcalls << " " << mfcn.NumOfCalls() << " " << st.NFcn() << std::endl;
+
#ifdef WARNINGMSG
- std::cout<<"MnHesse: maximum number of allowed function calls exhausted."<<std::endl;
- std::cout<<"MnHesse fails and will return diagonal matrix "<<std::endl;
+ std::cout<<"MnHesse: maximum number of allowed function calls exhausted."<<std::endl;
+ std::cout<<"MnHesse fails and will return diagonal matrix "<<std::endl;
#endif
-
- for(unsigned int j = 0; j < n; j++) {
- double tmp = g2(j) < prec.Eps2() ? 1. : 1./g2(j);
- vhmat(j,j) = tmp < prec.Eps2() ? 1. : tmp;
+
+ for(unsigned int j = 0; j < n; j++) {
+ double tmp = g2(j) < prec.Eps2() ? 1. : 1./g2(j);
+ vhmat(j,j) = tmp < prec.Eps2() ? 1. : tmp;
+ }
+
+ return MinimumState(st.Parameters(), MinimumError(vhmat, MinimumError::MnHesseFailed()), st.Gradient(), st.Edm(), mfcn.NumOfCalls());
}
- return MinimumState(st.Parameters(), MinimumError(vhmat, MinimumError::MnHesseFailed()), st.Gradient(), st.Edm(), mfcn.NumOfCalls());
- }
-
- }
-
- if(fStrategy.Strategy() > 0) {
- // refine first derivative
- HessianGradientCalculator hgc(mfcn, trafo, fStrategy);
- FunctionGradient gr = hgc(st.Parameters(), FunctionGradient(grd, g2, gst));
- grd = gr.Grad();
- }
-
- //off-diagonal Elements
- for(unsigned int i = 0; i < n; i++) {
- x(i) += dirin(i);
- for(unsigned int j = i+1; j < n; j++) {
- x(j) += dirin(j);
- double fs1 = mfcn(x);
- double elem = (fs1 + amin - yy(i) - yy(j))/(dirin(i)*dirin(j));
- vhmat(i,j) = elem;
- x(j) -= dirin(j);
- }
- x(i) -= dirin(i);
- }
-
- //verify if matrix pos-def (still 2nd derivative)
- MinimumError tmp = MnPosDef()(MinimumError(vhmat,1.), prec);
- vhmat = tmp.InvHessian();
- int ifail = Invert(vhmat);
- if(ifail != 0) {
-
+ }
+
+ if(fStrategy.Strategy() > 0) {
+ // refine first derivative
+ HessianGradientCalculator hgc(mfcn, trafo, fStrategy);
+ FunctionGradient gr = hgc(st.Parameters(), FunctionGradient(grd, g2, gst));
+ grd = gr.Grad();
+ }
+
+ //off-diagonal Elements
+ for(unsigned int i = 0; i < n; i++) {
+ x(i) += dirin(i);
+ for(unsigned int j = i+1; j < n; j++) {
+ x(j) += dirin(j);
+ double fs1 = mfcn(x);
+ double elem = (fs1 + amin - yy(i) - yy(j))/(dirin(i)*dirin(j));
+ vhmat(i,j) = elem;
+ x(j) -= dirin(j);
+ }
+ x(i) -= dirin(i);
+ }
+
+ //verify if matrix pos-def (still 2nd derivative)
+ MinimumError tmp = MnPosDef()(MinimumError(vhmat,1.), prec);
+ vhmat = tmp.InvHessian();
+ int ifail = Invert(vhmat);
+ if(ifail != 0) {
+
#ifdef WARNINGMSG
- std::cout<<"MnHesse: matrix inversion fails!"<<std::endl;
- std::cout<<"MnHesse fails and will return diagonal matrix."<<std::endl;
+ std::cout<<"MnHesse: matrix inversion fails!"<<std::endl;
+ std::cout<<"MnHesse fails and will return diagonal matrix."<<std::endl;
#endif
-
- MnAlgebraicSymMatrix tmpsym(vhmat.Nrow());
- for(unsigned int j = 0; j < n; j++) {
- double tmp = g2(j) < prec.Eps2() ? 1. : 1./g2(j);
- tmpsym(j,j) = tmp < prec.Eps2() ? 1. : tmp;
- }
-
- return MinimumState(st.Parameters(), MinimumError(tmpsym, MinimumError::MnHesseFailed()), st.Gradient(), st.Edm(), mfcn.NumOfCalls());
- }
-
- FunctionGradient gr(grd, g2, gst);
-
- // needed this ? (if posdef and inversion ok continue. it is like this in the Fortran version
-// if(tmp.IsMadePosDef()) {
-// std::cout<<"MnHesse: matrix is invalid!"<<std::endl;
-// std::cout<<"MnHesse: matrix is not pos. def.!"<<std::endl;
-// std::cout<<"MnHesse: matrix was forced pos. def."<<std::endl;
-// return MinimumState(st.Parameters(), MinimumError(vhmat, MinimumError::MnMadePosDef()), gr, st.Edm(), mfcn.NumOfCalls());
-// }
-
- //calculate edm
- MinimumError err(vhmat, 0.);
- VariableMetricEDMEstimator estim;
- double edm = estim.Estimate(gr, err);
-
- return MinimumState(st.Parameters(), err, gr, edm, mfcn.NumOfCalls());
+
+ MnAlgebraicSymMatrix tmpsym(vhmat.Nrow());
+ for(unsigned int j = 0; j < n; j++) {
+ double tmp = g2(j) < prec.Eps2() ? 1. : 1./g2(j);
+ tmpsym(j,j) = tmp < prec.Eps2() ? 1. : tmp;
+ }
+
+ return MinimumState(st.Parameters(), MinimumError(tmpsym, MinimumError::MnHesseFailed()), st.Gradient(), st.Edm(), mfcn.NumOfCalls());
+ }
+
+ FunctionGradient gr(grd, g2, gst);
+
+ // needed this ? (if posdef and inversion ok continue. it is like this in the Fortran version
+ // if(tmp.IsMadePosDef()) {
+ // std::cout<<"MnHesse: matrix is invalid!"<<std::endl;
+ // std::cout<<"MnHesse: matrix is not pos. def.!"<<std::endl;
+ // std::cout<<"MnHesse: matrix was forced pos. def."<<std::endl;
+ // return MinimumState(st.Parameters(), MinimumError(vhmat, MinimumError::MnMadePosDef()), gr, st.Edm(), mfcn.NumOfCalls());
+ // }
+
+ //calculate edm
+ MinimumError err(vhmat, 0.);
+ VariableMetricEDMEstimator estim;
+ double edm = estim.Estimate(gr, err);
+
+ return MinimumState(st.Parameters(), err, gr, edm, mfcn.NumOfCalls());
}
/*
-MinimumError MnHesse::Hessian(const MnFcn& mfcn, const MinimumState& st, const MnUserTransformation& trafo) const {
-
- const MnMachinePrecision& prec = trafo.Precision();
- // make sure starting at the right place
- double amin = mfcn(st.Vec());
-// if(fabs(amin - st.Fval()) > prec.Eps2()) std::cout<<"function Value differs from amin by "<<amin - st.Fval()<<std::endl;
-
- double aimsag = sqrt(prec.Eps2())*(fabs(amin)+mfcn.Up());
-
- // diagonal Elements first
-
- unsigned int n = st.Parameters().Vec().size();
- MnAlgebraicSymMatrix vhmat(n);
- MnAlgebraicVector g2 = st.Gradient().G2();
- MnAlgebraicVector gst = st.Gradient().Gstep();
- MnAlgebraicVector grd = st.Gradient().Grad();
- MnAlgebraicVector dirin = st.Gradient().Gstep();
- MnAlgebraicVector yy(n);
- MnAlgebraicVector x = st.Parameters().Vec();
-
- for(unsigned int i = 0; i < n; i++) {
-
- double xtf = x(i);
- double dmin = 8.*prec.Eps2()*fabs(xtf);
- double d = fabs(gst(i));
- if(d < dmin) d = dmin;
- for(int icyc = 0; icyc < Ncycles(); icyc++) {
- double sag = 0.;
- double fs1 = 0.;
- double fs2 = 0.;
- for(int multpy = 0; multpy < 5; multpy++) {
- x(i) = xtf + d;
- fs1 = mfcn(x);
- x(i) = xtf - d;
- fs2 = mfcn(x);
- x(i) = xtf;
- sag = 0.5*(fs1+fs2-2.*amin);
- if(sag > prec.Eps2()) break;
- if(trafo.Parameter(i).HasLimits()) {
- if(d > 0.5) {
- std::cout<<"second derivative zero for Parameter "<<i<<std::endl;
- std::cout<<"return diagonal matrix "<<std::endl;
- for(unsigned int j = 0; j < n; j++) {
- vhmat(j,j) = (g2(j) < prec.Eps2() ? 1. : 1./g2(j));
- return MinimumError(vhmat, 1., false);
- }
- }
- d *= 10.;
- if(d > 0.5) d = 0.51;
- continue;
- }
- d *= 10.;
- }
- if(sag < prec.Eps2()) {
- std::cout<<"MnHesse: internal loop exhausted, return diagonal matrix."<<std::endl;
- for(unsigned int i = 0; i < n; i++)
- vhmat(i,i) = (g2(i) < prec.Eps2() ? 1. : 1./g2(i));
- return MinimumError(vhmat, 1., false);
- }
- double g2bfor = g2(i);
- g2(i) = 2.*sag/(d*d);
- grd(i) = (fs1-fs2)/(2.*d);
- gst(i) = d;
- dirin(i) = d;
- yy(i) = fs1;
- double dlast = d;
- d = sqrt(2.*aimsag/fabs(g2(i)));
- if(trafo.Parameter(i).HasLimits()) d = std::min(0.5, d);
- if(d < dmin) d = dmin;
-
- // see if converged
- if(fabs((d-dlast)/d) < Tolerstp()) break;
- if(fabs((g2(i)-g2bfor)/g2(i)) < TolerG2()) break;
- d = std::min(d, 10.*dlast);
- d = std::max(d, 0.1*dlast);
+ MinimumError MnHesse::Hessian(const MnFcn& mfcn, const MinimumState& st, const MnUserTransformation& trafo) const {
+
+ const MnMachinePrecision& prec = trafo.Precision();
+ // make sure starting at the right place
+ double amin = mfcn(st.Vec());
+ // if(fabs(amin - st.Fval()) > prec.Eps2()) std::cout<<"function Value differs from amin by "<<amin - st.Fval()<<std::endl;
+
+ double aimsag = sqrt(prec.Eps2())*(fabs(amin)+mfcn.Up());
+
+ // diagonal Elements first
+
+ unsigned int n = st.Parameters().Vec().size();
+ MnAlgebraicSymMatrix vhmat(n);
+ MnAlgebraicVector g2 = st.Gradient().G2();
+ MnAlgebraicVector gst = st.Gradient().Gstep();
+ MnAlgebraicVector grd = st.Gradient().Grad();
+ MnAlgebraicVector dirin = st.Gradient().Gstep();
+ MnAlgebraicVector yy(n);
+ MnAlgebraicVector x = st.Parameters().Vec();
+
+ for(unsigned int i = 0; i < n; i++) {
+
+ double xtf = x(i);
+ double dmin = 8.*prec.Eps2()*fabs(xtf);
+ double d = fabs(gst(i));
+ if(d < dmin) d = dmin;
+ for(int icyc = 0; icyc < Ncycles(); icyc++) {
+ double sag = 0.;
+ double fs1 = 0.;
+ double fs2 = 0.;
+ for(int multpy = 0; multpy < 5; multpy++) {
+ x(i) = xtf + d;
+ fs1 = mfcn(x);
+ x(i) = xtf - d;
+ fs2 = mfcn(x);
+ x(i) = xtf;
+ sag = 0.5*(fs1+fs2-2.*amin);
+ if(sag > prec.Eps2()) break;
+ if(trafo.Parameter(i).HasLimits()) {
+ if(d > 0.5) {
+ std::cout<<"second derivative zero for Parameter "<<i<<std::endl;
+ std::cout<<"return diagonal matrix "<<std::endl;
+ for(unsigned int j = 0; j < n; j++) {
+ vhmat(j,j) = (g2(j) < prec.Eps2() ? 1. : 1./g2(j));
+ return MinimumError(vhmat, 1., false);
+ }
+ }
+ d *= 10.;
+ if(d > 0.5) d = 0.51;
+ continue;
+ }
+ d *= 10.;
+ }
+ if(sag < prec.Eps2()) {
+ std::cout<<"MnHesse: internal loop exhausted, return diagonal matrix."<<std::endl;
+ for(unsigned int i = 0; i < n; i++)
+ vhmat(i,i) = (g2(i) < prec.Eps2() ? 1. : 1./g2(i));
+ return MinimumError(vhmat, 1., false);
+ }
+ double g2bfor = g2(i);
+ g2(i) = 2.*sag/(d*d);
+ grd(i) = (fs1-fs2)/(2.*d);
+ gst(i) = d;
+ dirin(i) = d;
+ yy(i) = fs1;
+ double dlast = d;
+ d = sqrt(2.*aimsag/fabs(g2(i)));
+ if(trafo.Parameter(i).HasLimits()) d = std::min(0.5, d);
+ if(d < dmin) d = dmin;
+
+ // see if converged
+ if(fabs((d-dlast)/d) < Tolerstp()) break;
+ if(fabs((g2(i)-g2bfor)/g2(i)) < TolerG2()) break;
+ d = std::min(d, 10.*dlast);
+ d = std::max(d, 0.1*dlast);
+ }
+ vhmat(i,i) = g2(i);
}
- vhmat(i,i) = g2(i);
- }
-
- //off-diagonal Elements
- for(unsigned int i = 0; i < n; i++) {
- x(i) += dirin(i);
- for(unsigned int j = i+1; j < n; j++) {
- x(j) += dirin(j);
- double fs1 = mfcn(x);
- double elem = (fs1 + amin - yy(i) - yy(j))/(dirin(i)*dirin(j));
- vhmat(i,j) = elem;
- x(j) -= dirin(j);
+
+ //off-diagonal Elements
+ for(unsigned int i = 0; i < n; i++) {
+ x(i) += dirin(i);
+ for(unsigned int j = i+1; j < n; j++) {
+ x(j) += dirin(j);
+ double fs1 = mfcn(x);
+ double elem = (fs1 + amin - yy(i) - yy(j))/(dirin(i)*dirin(j));
+ vhmat(i,j) = elem;
+ x(j) -= dirin(j);
+ }
+ x(i) -= dirin(i);
}
- x(i) -= dirin(i);
- }
-
- return MinimumError(vhmat, 0.);
-}
-*/
+
+ return MinimumError(vhmat, 0.);
+ }
+ */
} // namespace Minuit2
diff --git a/minuit2/src/MnLineSearch.cxx b/minuit2/src/MnLineSearch.cxx
index cfea6e6f3bf..fb1fa75fae1 100644
--- a/minuit2/src/MnLineSearch.cxx
+++ b/minuit2/src/MnLineSearch.cxx
@@ -1,4 +1,4 @@
-// @(#)root/minuit2:$Name: $:$Id: MnLineSearch.cxx,v 1.1 2005/11/29 14:43:31 moneta Exp $
+// @(#)root/minuit2:$Name: $:$Id: MnLineSearch.cxx,v 1.2 2006/06/26 11:03:55 moneta Exp $
// Authors: M. Winkler, F. James, L. Moneta, A. Zsenei 2003-2005
/**********************************************************************
@@ -238,6 +238,6 @@ MnParabolaPoint MnLineSearch::operator()(const MnFcn& fcn, const MinimumParamete
return MnParabolaPoint(xvmin, fvmin);
}
- } // namespace Minuit2
+} // namespace Minuit2
} // namespace ROOT
diff --git a/minuit2/src/MnMachinePrecision.cxx b/minuit2/src/MnMachinePrecision.cxx
index 9013d2fe921..7a70a965a49 100644
--- a/minuit2/src/MnMachinePrecision.cxx
+++ b/minuit2/src/MnMachinePrecision.cxx
@@ -1,4 +1,4 @@
-// @(#)root/minuit2:$Name: $:$Id: MnMachinePrecision.cpp,v 1.3.6.3 2005/11/29 11:08:35 moneta Exp $
+// @(#)root/minuit2:$Name: $:$Id: MnMachinePrecision.cxx,v 1.1 2005/11/29 14:43:31 moneta Exp $
// Authors: M. Winkler, F. James, L. Moneta, A. Zsenei 2003-2005
/**********************************************************************
@@ -15,38 +15,39 @@ namespace ROOT {
namespace Minuit2 {
-MnMachinePrecision::MnMachinePrecision() : fEpsMac(4.0E-7),
- fEpsMa2(2.*sqrt(4.0E-7)) {
-
- //determine machine precision
- /*
- char e[] = {"e"};
- fEpsMac = 8.*dlamch_(e);
- fEpsMa2 = 2.*sqrt(fEpsMac);
- */
-
-// std::cout<<"machine precision eps: "<<Eps()<<std::endl;
-
- MnTiny mytiny;
-
- //calculate machine precision
- double epstry = 0.5;
- double epsbak = 0.;
- double epsp1 = 0.;
- double one = 1.0;
- for(int i = 0; i < 100; i++) {
- epstry *= 0.5;
- epsp1 = one + epstry;
- epsbak = mytiny(epsp1);
- if(epsbak < epstry) {
- fEpsMac = 8.*epstry;
- fEpsMa2 = 2.*sqrt(fEpsMac);
- break;
- }
- }
-
+MnMachinePrecision::MnMachinePrecision() :
+ fEpsMac(4.0E-7),
+ fEpsMa2(2.*sqrt(4.0E-7)) {
+
+ //determine machine precision
+ /*
+ char e[] = {"e"};
+ fEpsMac = 8.*dlamch_(e);
+ fEpsMa2 = 2.*sqrt(fEpsMac);
+ */
+
+ // std::cout<<"machine precision eps: "<<Eps()<<std::endl;
+
+ MnTiny mytiny;
+
+ //calculate machine precision
+ double epstry = 0.5;
+ double epsbak = 0.;
+ double epsp1 = 0.;
+ double one = 1.0;
+ for(int i = 0; i < 100; i++) {
+ epstry *= 0.5;
+ epsp1 = one + epstry;
+ epsbak = mytiny(epsp1);
+ if(epsbak < epstry) {
+ fEpsMac = 8.*epstry;
+ fEpsMa2 = 2.*sqrt(fEpsMac);
+ break;
+ }
+ }
+
}
-
- } // namespace Minuit2
-
+
+ } // namespace Minuit2
+
} // namespace ROOT
diff --git a/minuit2/src/MnMinos.cxx b/minuit2/src/MnMinos.cxx
index b4f3ccda0fc..f419a31955f 100644
--- a/minuit2/src/MnMinos.cxx
+++ b/minuit2/src/MnMinos.cxx
@@ -1,4 +1,4 @@
-// @(#)root/minuit2:$Name: $:$Id: MnMinos.cxx,v 1.1 2005/11/29 14:43:31 moneta Exp $
+// @(#)root/minuit2:$Name: $:$Id: MnMinos.cxx,v 1.2 2006/04/12 16:30:31 moneta Exp $
// Authors: M. Winkler, F. James, L. Moneta, A. Zsenei 2003-2005
/**********************************************************************
@@ -21,151 +21,155 @@ namespace ROOT {
std::pair<double,double> MnMinos::operator()(unsigned int par, unsigned int maxcalls) const {
- MinosError mnerr = Minos(par, maxcalls);
- return mnerr();
+ // do Minos analysis given the parameter index returning a pair for (lower,upper) errors
+ MinosError mnerr = Minos(par, maxcalls);
+ return mnerr();
}
double MnMinos::Lower(unsigned int par, unsigned int maxcalls) const {
-
- MnUserParameterState upar = fMinimum.UserState();
- double err = fMinimum.UserState().Error(par);
-
- MnCross aopt = Loval(par, maxcalls);
-
- double Lower = aopt.IsValid() ? -1.*err*(1.+ aopt.Value()) : (aopt.AtLimit() ? upar.Parameter(par).LowerLimit() : upar.Value(par));
-
- return Lower;
+ // get lower error for parameter par
+ MnUserParameterState upar = fMinimum.UserState();
+ double err = fMinimum.UserState().Error(par);
+
+ MnCross aopt = Loval(par, maxcalls);
+
+ double Lower = aopt.IsValid() ? -1.*err*(1.+ aopt.Value()) : (aopt.AtLimit() ? upar.Parameter(par).LowerLimit() : upar.Value(par));
+
+ return Lower;
}
double MnMinos::Upper(unsigned int par, unsigned int maxcalls) const {
- MnCross aopt = Upval(par, maxcalls);
-
- MnUserParameterState upar = fMinimum.UserState();
- double err = fMinimum.UserState().Error(par);
-
- double Upper = aopt.IsValid() ? err*(1.+ aopt.Value()) : (aopt.AtLimit() ? upar.Parameter(par).UpperLimit() : upar.Value(par));
-
- return Upper;
+ // upper error for parameter par
+ MnCross aopt = Upval(par, maxcalls);
+
+ MnUserParameterState upar = fMinimum.UserState();
+ double err = fMinimum.UserState().Error(par);
+
+ double Upper = aopt.IsValid() ? err*(1.+ aopt.Value()) : (aopt.AtLimit() ? upar.Parameter(par).UpperLimit() : upar.Value(par));
+
+ return Upper;
}
MinosError MnMinos::Minos(unsigned int par, unsigned int maxcalls) const {
- assert(fMinimum.IsValid());
- assert(!fMinimum.UserState().Parameter(par).IsFixed());
- assert(!fMinimum.UserState().Parameter(par).IsConst());
-
- MnCross up = Upval(par, maxcalls);
- MnCross lo = Loval(par, maxcalls);
-
- return MinosError(par, fMinimum.UserState().Value(par), lo, up);
+ assert(fMinimum.IsValid());
+ assert(!fMinimum.UserState().Parameter(par).IsFixed());
+ assert(!fMinimum.UserState().Parameter(par).IsConst());
+
+ MnCross up = Upval(par, maxcalls);
+ MnCross lo = Loval(par, maxcalls);
+
+ return MinosError(par, fMinimum.UserState().Value(par), lo, up);
}
MnCross MnMinos::Upval(unsigned int par, unsigned int maxcalls) const {
- assert(fMinimum.IsValid());
- assert(!fMinimum.UserState().Parameter(par).IsFixed());
- assert(!fMinimum.UserState().Parameter(par).IsConst());
- if(maxcalls == 0) {
- unsigned int nvar = fMinimum.UserState().VariableParameters();
- maxcalls = 2*(nvar+1)*(200 + 100*nvar + 5*nvar*nvar);
- }
-
- std::vector<unsigned int> para(1, par);
-
- MnUserParameterState upar = fMinimum.UserState();
- double err = upar.Error(par);
- double val = upar.Value(par) + err;
- std::vector<double> xmid(1, val);
- std::vector<double> xdir(1, err);
-
- double up = fFCN.Up();
- unsigned int ind = upar.IntOfExt(par);
- MnAlgebraicSymMatrix m = fMinimum.Error().Matrix();
- double xunit = sqrt(up/err);
- for(unsigned int i = 0; i < m.Nrow(); i++) {
- if(i == ind) continue;
- double xdev = xunit*m(ind,i);
- unsigned int ext = upar.ExtOfInt(i);
- upar.SetValue(ext, upar.Value(ext) + xdev);
- }
-
- upar.Fix(par);
- upar.SetValue(par, val);
-
-// double edmmax = 0.5*0.1*fFCN.Up()*1.e-3;
- double toler = 0.1;
- MnFunctionCross cross(fFCN, upar, fMinimum.Fval(), fStrategy);
-
- MnCross aopt = cross(para, xmid, xdir, toler, maxcalls);
-
-// std::cout<<"aopt= "<<aopt.Value()<<std::endl;
-
+ // get crossing value in the upper parameter direction
+ assert(fMinimum.IsValid());
+ assert(!fMinimum.UserState().Parameter(par).IsFixed());
+ assert(!fMinimum.UserState().Parameter(par).IsConst());
+ if(maxcalls == 0) {
+ unsigned int nvar = fMinimum.UserState().VariableParameters();
+ maxcalls = 2*(nvar+1)*(200 + 100*nvar + 5*nvar*nvar);
+ }
+
+ std::vector<unsigned int> para(1, par);
+
+ MnUserParameterState upar = fMinimum.UserState();
+ double err = upar.Error(par);
+ double val = upar.Value(par) + err;
+ std::vector<double> xmid(1, val);
+ std::vector<double> xdir(1, err);
+
+ double up = fFCN.Up();
+ unsigned int ind = upar.IntOfExt(par);
+ MnAlgebraicSymMatrix m = fMinimum.Error().Matrix();
+ double xunit = sqrt(up/err);
+ for(unsigned int i = 0; i < m.Nrow(); i++) {
+ if(i == ind) continue;
+ double xdev = xunit*m(ind,i);
+ unsigned int ext = upar.ExtOfInt(i);
+ upar.SetValue(ext, upar.Value(ext) + xdev);
+ }
+
+ upar.Fix(par);
+ upar.SetValue(par, val);
+
+ // double edmmax = 0.5*0.1*fFCN.Up()*1.e-3;
+ double toler = 0.1;
+ MnFunctionCross cross(fFCN, upar, fMinimum.Fval(), fStrategy);
+
+ MnCross aopt = cross(para, xmid, xdir, toler, maxcalls);
+
+ // std::cout<<"aopt= "<<aopt.Value()<<std::endl;
+
#ifdef WARNINGMSG
- if(aopt.AtLimit())
- std::cout<<"MnMinos Parameter "<<par<<" is at Upper limit."<<std::endl;
- if(aopt.AtMaxFcn())
- std::cout<<"MnMinos maximum number of function calls exceeded for Parameter "<<par<<std::endl;
- if(aopt.NewMinimum())
- std::cout<<"MnMinos new Minimum found while looking for Parameter "<<par<<std::endl;
- if(!aopt.IsValid())
- std::cout<<"MnMinos could not find Upper Value for Parameter "<<par<<"."<<std::endl;
+ if(aopt.AtLimit())
+ std::cout<<"MnMinos Parameter "<<par<<" is at Upper limit."<<std::endl;
+ if(aopt.AtMaxFcn())
+ std::cout<<"MnMinos maximum number of function calls exceeded for Parameter "<<par<<std::endl;
+ if(aopt.NewMinimum())
+ std::cout<<"MnMinos new Minimum found while looking for Parameter "<<par<<std::endl;
+ if(!aopt.IsValid())
+ std::cout<<"MnMinos could not find Upper Value for Parameter "<<par<<"."<<std::endl;
#endif
-
- return aopt;
+
+ return aopt;
}
MnCross MnMinos::Loval(unsigned int par, unsigned int maxcalls) const {
- assert(fMinimum.IsValid());
- assert(!fMinimum.UserState().Parameter(par).IsFixed());
- assert(!fMinimum.UserState().Parameter(par).IsConst());
- if(maxcalls == 0) {
- unsigned int nvar = fMinimum.UserState().VariableParameters();
- maxcalls = 2*(nvar+1)*(200 + 100*nvar + 5*nvar*nvar);
- }
- std::vector<unsigned int> para(1, par);
-
- MnUserParameterState upar = fMinimum.UserState();
- double err = upar.Error(par);
- double val = upar.Value(par) - err;
- std::vector<double> xmid(1, val);
- std::vector<double> xdir(1, -err);
-
- double up = fFCN.Up();
- unsigned int ind = upar.IntOfExt(par);
- MnAlgebraicSymMatrix m = fMinimum.Error().Matrix();
- double xunit = sqrt(up/err);
- for(unsigned int i = 0; i < m.Nrow(); i++) {
- if(i == ind) continue;
- double xdev = xunit*m(ind,i);
- unsigned int ext = upar.ExtOfInt(i);
- upar.SetValue(ext, upar.Value(ext) - xdev);
- }
-
- upar.Fix(par);
- upar.SetValue(par, val);
-
-// double edmmax = 0.5*0.1*fFCN.Up()*1.e-3;
- double toler = 0.1;
- MnFunctionCross cross(fFCN, upar, fMinimum.Fval(), fStrategy);
-
- MnCross aopt = cross(para, xmid, xdir, toler, maxcalls);
-
-// std::cout<<"aopt= "<<aopt.Value()<<std::endl;
-
+ // return crossing in the lower parameter direction
+ assert(fMinimum.IsValid());
+ assert(!fMinimum.UserState().Parameter(par).IsFixed());
+ assert(!fMinimum.UserState().Parameter(par).IsConst());
+ if(maxcalls == 0) {
+ unsigned int nvar = fMinimum.UserState().VariableParameters();
+ maxcalls = 2*(nvar+1)*(200 + 100*nvar + 5*nvar*nvar);
+ }
+ std::vector<unsigned int> para(1, par);
+
+ MnUserParameterState upar = fMinimum.UserState();
+ double err = upar.Error(par);
+ double val = upar.Value(par) - err;
+ std::vector<double> xmid(1, val);
+ std::vector<double> xdir(1, -err);
+
+ double up = fFCN.Up();
+ unsigned int ind = upar.IntOfExt(par);
+ MnAlgebraicSymMatrix m = fMinimum.Error().Matrix();
+ double xunit = sqrt(up/err);
+ for(unsigned int i = 0; i < m.Nrow(); i++) {
+ if(i == ind) continue;
+ double xdev = xunit*m(ind,i);
+ unsigned int ext = upar.ExtOfInt(i);
+ upar.SetValue(ext, upar.Value(ext) - xdev);
+ }
+
+ upar.Fix(par);
+ upar.SetValue(par, val);
+
+ // double edmmax = 0.5*0.1*fFCN.Up()*1.e-3;
+ double toler = 0.1;
+ MnFunctionCross cross(fFCN, upar, fMinimum.Fval(), fStrategy);
+
+ MnCross aopt = cross(para, xmid, xdir, toler, maxcalls);
+
+ // std::cout<<"aopt= "<<aopt.Value()<<std::endl;
+
#ifdef WARNINGMSG
- if(aopt.AtLimit())
- std::cout<<"MnMinos Parameter "<<par<<" is at Lower limit."<<std::endl;
- if(aopt.AtMaxFcn())
- std::cout<<"MnMinos maximum number of function calls exceeded for Parameter "<<par<<std::endl;
- if(aopt.NewMinimum())
- std::cout<<"MnMinos new Minimum found while looking for Parameter "<<par<<std::endl;
- if(!aopt.IsValid())
- std::cout<<"MnMinos could not find Lower Value for Parameter "<<par<<"."<<std::endl;
+ if(aopt.AtLimit())
+ std::cout<<"MnMinos Parameter "<<par<<" is at Lower limit."<<std::endl;
+ if(aopt.AtMaxFcn())
+ std::cout<<"MnMinos maximum number of function calls exceeded for Parameter "<<par<<std::endl;
+ if(aopt.NewMinimum())
+ std::cout<<"MnMinos new Minimum found while looking for Parameter "<<par<<std::endl;
+ if(!aopt.IsValid())
+ std::cout<<"MnMinos could not find Lower Value for Parameter "<<par<<"."<<std::endl;
#endif
-
- return aopt;
-
+
+ return aopt;
+
}
- } // namespace Minuit2
+ } // namespace Minuit2
} // namespace ROOT
diff --git a/minuit2/src/MnParabolaFactory.cxx b/minuit2/src/MnParabolaFactory.cxx
index 55e357ebf0e..08d3cc0e1f7 100644
--- a/minuit2/src/MnParabolaFactory.cxx
+++ b/minuit2/src/MnParabolaFactory.cxx
@@ -1,4 +1,4 @@
-// @(#)root/minuit2:$Name: $:$Id: MnParabolaFactory.cxx,v 1.1 2005/11/29 14:43:31 moneta Exp $
+// @(#)root/minuit2:$Name: $:$Id: MnParabolaFactory.cxx,v 1.2 2006/06/26 11:03:55 moneta Exp $
// Authors: M. Winkler, F. James, L. Moneta, A. Zsenei 2003-2005
/**********************************************************************
@@ -69,7 +69,7 @@ MnParabola MnParabolaFactory::operator()(const MnParabolaPoint& p1,
double c = y1 - a*xx1 - b*x1;
return MnParabola(a, b, c);
- }
+}
} // namespace Minuit2
diff --git a/minuit2/src/MnParameterScan.cxx b/minuit2/src/MnParameterScan.cxx
index 563be09bd8e..f6bf51bcfd6 100644
--- a/minuit2/src/MnParameterScan.cxx
+++ b/minuit2/src/MnParameterScan.cxx
@@ -1,4 +1,4 @@
-// @(#)root/minuit2:$Name: $:$Id: MnParameterScan.cpp,v 1.1.6.4 2005/11/29 11:08:35 moneta Exp $
+// @(#)root/minuit2:$Name: $:$Id: MnParameterScan.cxx,v 1.1 2005/11/29 14:43:31 moneta Exp $
// Authors: M. Winkler, F. James, L. Moneta, A. Zsenei 2003-2005
/**********************************************************************
@@ -20,49 +20,50 @@ MnParameterScan::MnParameterScan(const FCNBase& fcn, const MnUserParameters& par
MnParameterScan::MnParameterScan(const FCNBase& fcn, const MnUserParameters& par, double fval) : fFCN(fcn), fParameters(par), fAmin(fval) {}
std::vector<std::pair<double, double> > MnParameterScan::operator()(unsigned int par, unsigned int maxsteps, double low, double high) {
-
- if(maxsteps > 101) maxsteps = 101;
- std::vector<std::pair<double, double> > result; result.reserve(maxsteps+1);
- std::vector<double> params = fParameters.Params();
- result.push_back(std::pair<double, double>(params[par], fAmin));
-
- if(low > high) return result;
- if(maxsteps < 2) return result;
-
- if(low == 0. && high == 0.) {
- low = params[par] - 2.*fParameters.Error(par);
- high = params[par] + 2.*fParameters.Error(par);
- }
-
- if(low == 0. && high == 0. && fParameters.Parameter(par).HasLimits()) {
- if(fParameters.Parameter(par).HasLowerLimit())
- low = fParameters.Parameter(par).LowerLimit();
- if(fParameters.Parameter(par).HasUpperLimit())
- high = fParameters.Parameter(par).UpperLimit();
- }
-
- if(fParameters.Parameter(par).HasLimits()) {
- if(fParameters.Parameter(par).HasLowerLimit())
- low = std::max(low, fParameters.Parameter(par).LowerLimit());
- if(fParameters.Parameter(par).HasUpperLimit())
- high = std::min(high, fParameters.Parameter(par).UpperLimit());
- }
-
- double x0 = low;
- double stp = (high - low)/double(maxsteps - 1);
- for(unsigned int i = 0; i < maxsteps; i++) {
- params[par] = x0 + double(i)*stp;
- double fval = fFCN(params);
- if(fval < fAmin) {
- fParameters.SetValue(par, params[par]);
- fAmin = fval;
- }
- result.push_back(std::pair<double, double>(params[par], fval));
- }
-
- return result;
+ // do the scan for parameter par between low and high values
+
+ if(maxsteps > 101) maxsteps = 101;
+ std::vector<std::pair<double, double> > result; result.reserve(maxsteps+1);
+ std::vector<double> params = fParameters.Params();
+ result.push_back(std::pair<double, double>(params[par], fAmin));
+
+ if(low > high) return result;
+ if(maxsteps < 2) return result;
+
+ if(low == 0. && high == 0.) {
+ low = params[par] - 2.*fParameters.Error(par);
+ high = params[par] + 2.*fParameters.Error(par);
+ }
+
+ if(low == 0. && high == 0. && fParameters.Parameter(par).HasLimits()) {
+ if(fParameters.Parameter(par).HasLowerLimit())
+ low = fParameters.Parameter(par).LowerLimit();
+ if(fParameters.Parameter(par).HasUpperLimit())
+ high = fParameters.Parameter(par).UpperLimit();
+ }
+
+ if(fParameters.Parameter(par).HasLimits()) {
+ if(fParameters.Parameter(par).HasLowerLimit())
+ low = std::max(low, fParameters.Parameter(par).LowerLimit());
+ if(fParameters.Parameter(par).HasUpperLimit())
+ high = std::min(high, fParameters.Parameter(par).UpperLimit());
+ }
+
+ double x0 = low;
+ double stp = (high - low)/double(maxsteps - 1);
+ for(unsigned int i = 0; i < maxsteps; i++) {
+ params[par] = x0 + double(i)*stp;
+ double fval = fFCN(params);
+ if(fval < fAmin) {
+ fParameters.SetValue(par, params[par]);
+ fAmin = fval;
+ }
+ result.push_back(std::pair<double, double>(params[par], fval));
+ }
+
+ return result;
}
- } // namespace Minuit2
+ } // namespace Minuit2
} // namespace ROOT
diff --git a/minuit2/src/MnPlot.cxx b/minuit2/src/MnPlot.cxx
index ad49fc2c0be..88c03446f79 100644
--- a/minuit2/src/MnPlot.cxx
+++ b/minuit2/src/MnPlot.cxx
@@ -1,4 +1,4 @@
-// @(#)root/minuit2:$Name: $:$Id: MnPlot.cpp,v 1.2.6.3 2005/11/29 11:08:35 moneta Exp $
+// @(#)root/minuit2:$Name: $:$Id: MnPlot.cxx,v 1.1 2005/11/29 14:43:31 moneta Exp $
// Authors: M. Winkler, F. James, L. Moneta, A. Zsenei 2003-2005
/**********************************************************************
@@ -17,40 +17,42 @@ namespace ROOT {
void mnplot(double* xpt, double* ypt, char* chpt, int nxypt, int npagwd, int npagln);
void MnPlot::operator()(const std::vector<std::pair<double,double> >& points) const {
- std::vector<double> x; x.reserve(points.size());
- std::vector<double> y; y.reserve(points.size());
- std::vector<char> chpt; chpt.reserve(points.size());
-
- for(std::vector<std::pair<double,double> >::const_iterator ipoint = points.begin(); ipoint != points.end(); ipoint++) {
- x.push_back((*ipoint).first);
- y.push_back((*ipoint).second);
- chpt.push_back('*');
- }
-
- mnplot(&(x.front()), &(y.front()), &(chpt.front()), points.size(), Width(), Length());
-
+ // call routine from Fortran minuit (mnplot) to plot the vector of (x,y) points
+ std::vector<double> x; x.reserve(points.size());
+ std::vector<double> y; y.reserve(points.size());
+ std::vector<char> chpt; chpt.reserve(points.size());
+
+ for(std::vector<std::pair<double,double> >::const_iterator ipoint = points.begin(); ipoint != points.end(); ipoint++) {
+ x.push_back((*ipoint).first);
+ y.push_back((*ipoint).second);
+ chpt.push_back('*');
+ }
+
+ mnplot(&(x.front()), &(y.front()), &(chpt.front()), points.size(), Width(), Length());
+
}
void MnPlot::operator()(double xmin, double ymin, const std::vector<std::pair<double,double> >& points) const {
- std::vector<double> x; x.reserve(points.size()+2);
- x.push_back(xmin);
- x.push_back(xmin);
- std::vector<double> y; y.reserve(points.size()+2);
- y.push_back(ymin);
- y.push_back(ymin);
- std::vector<char> chpt; chpt.reserve(points.size()+2);
- chpt.push_back(' ');
- chpt.push_back('X');
-
- for(std::vector<std::pair<double,double> >::const_iterator ipoint = points.begin(); ipoint != points.end(); ipoint++) {
- x.push_back((*ipoint).first);
- y.push_back((*ipoint).second);
- chpt.push_back('*');
- }
-
- mnplot(&(x.front()), &(y.front()), &(chpt.front()), points.size()+2, Width(), Length());
+ // call routine from Fortran minuit (mnplot) to plot the vector of (x,y) points + minimum values
+ std::vector<double> x; x.reserve(points.size()+2);
+ x.push_back(xmin);
+ x.push_back(xmin);
+ std::vector<double> y; y.reserve(points.size()+2);
+ y.push_back(ymin);
+ y.push_back(ymin);
+ std::vector<char> chpt; chpt.reserve(points.size()+2);
+ chpt.push_back(' ');
+ chpt.push_back('X');
+
+ for(std::vector<std::pair<double,double> >::const_iterator ipoint = points.begin(); ipoint != points.end(); ipoint++) {
+ x.push_back((*ipoint).first);
+ y.push_back((*ipoint).second);
+ chpt.push_back('*');
+ }
+
+ mnplot(&(x.front()), &(y.front()), &(chpt.front()), points.size()+2, Width(), Length());
}
- } // namespace Minuit2
+ } // namespace Minuit2
} // namespace ROOT
diff --git a/minuit2/src/MnPosDef.cxx b/minuit2/src/MnPosDef.cxx
index 23e6f6142da..8fb2f4958ae 100644
--- a/minuit2/src/MnPosDef.cxx
+++ b/minuit2/src/MnPosDef.cxx
@@ -1,4 +1,4 @@
-// @(#)root/minuit2:$Name: $:$Id: MnPosDef.cxx,v 1.1 2005/11/29 14:43:31 moneta Exp $
+// @(#)root/minuit2:$Name: $:$Id: MnPosDef.cxx,v 1.2 2006/04/12 16:30:31 moneta Exp $
// Authors: M. Winkler, F. James, L. Moneta, A. Zsenei 2003-2005
/**********************************************************************
@@ -26,79 +26,80 @@ LAVector eigenvalues(const LASymMatrix&);
MinimumState MnPosDef::operator()(const MinimumState& st, const MnMachinePrecision& prec) const {
-
- MinimumError err = (*this)(st.Error(), prec);
- return MinimumState(st.Parameters(), err, st.Gradient(), st.Edm(), st.NFcn());
+ // interface from minimum state
+ MinimumError err = (*this)(st.Error(), prec);
+ return MinimumState(st.Parameters(), err, st.Gradient(), st.Edm(), st.NFcn());
}
MinimumError MnPosDef::operator()(const MinimumError& e, const MnMachinePrecision& prec) const {
-
- MnAlgebraicSymMatrix err(e.InvHessian());
- if(err.size() == 1 && err(0,0) < prec.Eps()) {
- err(0,0) = 1.;
- return MinimumError(err, MinimumError::MnMadePosDef());
- }
- if(err.size() == 1 && err(0,0) > prec.Eps()) {
- return e;
- }
-// std::cout<<"MnPosDef init matrix= "<<err<<std::endl;
-
- double epspdf = std::max(1.e-6, prec.Eps2());
- double dgmin = err(0,0);
-
- for(unsigned int i = 0; i < err.Nrow(); i++) {
+ // make error matrix positive defined returning a new corrected minimum error state
+
+ MnAlgebraicSymMatrix err(e.InvHessian());
+ if(err.size() == 1 && err(0,0) < prec.Eps()) {
+ err(0,0) = 1.;
+ return MinimumError(err, MinimumError::MnMadePosDef());
+ }
+ if(err.size() == 1 && err(0,0) > prec.Eps()) {
+ return e;
+ }
+ // std::cout<<"MnPosDef init matrix= "<<err<<std::endl;
+
+ double epspdf = std::max(1.e-6, prec.Eps2());
+ double dgmin = err(0,0);
+
+ for(unsigned int i = 0; i < err.Nrow(); i++) {
#ifdef WARNINGMSG
- if(err(i,i) < prec.Eps2()) std::cout<<"negative or zero diagonal element "<<i<<" in covariance matrix"<<std::endl;
+ if(err(i,i) < prec.Eps2()) std::cout<<"negative or zero diagonal element "<<i<<" in covariance matrix"<<std::endl;
#endif
- if(err(i,i) < dgmin) dgmin = err(i,i);
- }
- double dg = 0.;
- if(dgmin < prec.Eps2()) {
- //dg = 1. + epspdf - dgmin;
- dg = 0.5 + epspdf - dgmin;
-// dg = 0.5*(1. + epspdf - dgmin);
+ if(err(i,i) < dgmin) dgmin = err(i,i);
+ }
+ double dg = 0.;
+ if(dgmin < prec.Eps2()) {
+ //dg = 1. + epspdf - dgmin;
+ dg = 0.5 + epspdf - dgmin;
+ // dg = 0.5*(1. + epspdf - dgmin);
#ifdef WARNINGMSG
- std::cout<<"added "<<dg<<" to diagonal of Error matrix"<<std::endl;
+ std::cout<<"added "<<dg<<" to diagonal of Error matrix"<<std::endl;
#endif
- //std::cout << "Error matrix " << err << std::endl;
- }
-
- MnAlgebraicVector s(err.Nrow());
- MnAlgebraicSymMatrix p(err.Nrow());
- for(unsigned int i = 0; i < err.Nrow(); i++) {
- err(i,i) += dg;
- if(err(i,i) < 0.) err(i,i) = 1.;
- s(i) = 1./sqrt(err(i,i));
- for(unsigned int j = 0; j <= i; j++) {
- p(i,j) = err(i,j)*s(i)*s(j);
- }
- }
-
- //std::cout<<"MnPosDef p: "<<p<<std::endl;
- MnAlgebraicVector eval = eigenvalues(p);
- double pmin = eval(0);
- double pmax = eval(eval.size() - 1);
- //std::cout<<"pmin= "<<pmin<<" pmax= "<<pmax<<std::endl;
- pmax = std::max(fabs(pmax), 1.);
- if(pmin > epspdf*pmax) return MinimumError(err, e.Dcovar());
-
- double padd = 0.001*pmax - pmin;
+ //std::cout << "Error matrix " << err << std::endl;
+ }
+
+ MnAlgebraicVector s(err.Nrow());
+ MnAlgebraicSymMatrix p(err.Nrow());
+ for(unsigned int i = 0; i < err.Nrow(); i++) {
+ err(i,i) += dg;
+ if(err(i,i) < 0.) err(i,i) = 1.;
+ s(i) = 1./sqrt(err(i,i));
+ for(unsigned int j = 0; j <= i; j++) {
+ p(i,j) = err(i,j)*s(i)*s(j);
+ }
+ }
+
+ //std::cout<<"MnPosDef p: "<<p<<std::endl;
+ MnAlgebraicVector eval = eigenvalues(p);
+ double pmin = eval(0);
+ double pmax = eval(eval.size() - 1);
+ //std::cout<<"pmin= "<<pmin<<" pmax= "<<pmax<<std::endl;
+ pmax = std::max(fabs(pmax), 1.);
+ if(pmin > epspdf*pmax) return MinimumError(err, e.Dcovar());
+
+ double padd = 0.001*pmax - pmin;
#ifdef DEBUG
- std::cout<<"eigenvalues: "<<std::endl;
+ std::cout<<"eigenvalues: "<<std::endl;
#endif
- for(unsigned int i = 0; i < err.Nrow(); i++) {
- err(i,i) *= (1. + padd);
+ for(unsigned int i = 0; i < err.Nrow(); i++) {
+ err(i,i) *= (1. + padd);
#ifdef DEBUG
- std::cout<<eval(i)<<std::endl;
+ std::cout<<eval(i)<<std::endl;
#endif
- }
-// std::cout<<"MnPosDef final matrix: "<<err<<std::endl;
+ }
+ // std::cout<<"MnPosDef final matrix: "<<err<<std::endl;
#ifdef WARNINGMSG
- std::cout<<"matrix forced pos-def by adding "<<padd<<" to diagonal"<<std::endl;
+ std::cout<<"matrix forced pos-def by adding "<<padd<<" to diagonal"<<std::endl;
#endif
- return MinimumError(err, MinimumError::MnMadePosDef());
+ return MinimumError(err, MinimumError::MnMadePosDef());
}
- } // namespace Minuit2
+ } // namespace Minuit2
} // namespace ROOT
diff --git a/minuit2/src/MnPrint.cxx b/minuit2/src/MnPrint.cxx
index 6cc8231d67d..ef24e8e24fa 100644
--- a/minuit2/src/MnPrint.cxx
+++ b/minuit2/src/MnPrint.cxx
@@ -1,4 +1,4 @@
-// @(#)root/minuit2:$Name: $:$Id: MnPrint.cpp,v 1.2.4.4 2005/11/29 11:08:35 moneta Exp $
+// @(#)root/minuit2:$Name: $:$Id: MnPrint.cxx,v 1.1 2005/11/29 14:43:31 moneta Exp $
// Authors: M. Winkler, F. James, L. Moneta, A. Zsenei 2003-2005
/**********************************************************************
@@ -29,277 +29,279 @@ namespace ROOT {
std::ostream& operator<<(std::ostream& os, const LAVector& vec) {
- os << "LAVector parameters:" << std::endl;
- {
- os << std::endl;
- int nrow = vec.size();
- for (int i = 0; i < nrow; i++) {
- os.precision(6); os.width(13);
- os << vec(i) << std::endl;
- }
- os << std::endl;
- }
- return os;
+ // print a vector
+ os << "LAVector parameters:" << std::endl;
+ {
+ os << std::endl;
+ int nrow = vec.size();
+ for (int i = 0; i < nrow; i++) {
+ os.precision(6); os.width(13);
+ os << vec(i) << std::endl;
+ }
+ os << std::endl;
+ }
+ return os;
}
std::ostream& operator<<(std::ostream& os, const LASymMatrix& matrix) {
- os << "LASymMatrix parameters:" << std::endl;
+ // print a matrix
+ os << "LASymMatrix parameters:" << std::endl;
{
- os << std::endl;
- int n = matrix.Nrow();
- for (int i = 0; i < n; i++) {
- for (int j = 0; j < n; j++) {
- os.precision(6); os.width(13); os << matrix(i,j);
- }
os << std::endl;
- }
- }
- return os;
+ int n = matrix.Nrow();
+ for (int i = 0; i < n; i++) {
+ for (int j = 0; j < n; j++) {
+ os.precision(6); os.width(13); os << matrix(i,j);
+ }
+ os << std::endl;
+ }
+ }
+ return os;
}
std::ostream& operator<<(std::ostream& os, const MnUserParameters& par) {
-
- os << std::endl;
-
- os << "# ext. |" << "| Name |" << "| type |" << "| Value |" << "| Error +/- " << std::endl;
-
- os << std::endl;
-
- bool atLoLim = false;
- bool atHiLim = false;
- for(std::vector<MinuitParameter>::const_iterator ipar = par.Parameters().begin(); ipar != par.Parameters().end(); ipar++) {
- os << std::setw(4) << (*ipar).Number() << std::setw(5) << "||";
- os << std::setw(10) << (*ipar).Name() << std::setw(3) << "||";
- if((*ipar).IsConst()) {
- os << " const ||" << std::setprecision(4) << std::setw(10) << (*ipar).Value() << " ||" << std::endl;
- } else if((*ipar).IsFixed()) {
- os << " fixed ||" << std::setprecision(4) << std::setw(10) << (*ipar).Value() << " ||" << std::endl;
- } else if((*ipar).HasLimits()) {
- if((*ipar).Error() > 0.) {
- os << " limited ||" << std::setprecision(4) << std::setw(10) << (*ipar).Value();
- if(fabs((*ipar).Value() - (*ipar).LowerLimit()) < par.Precision().Eps2()) {
- os <<"*";
- atLoLim = true;
- }
- if(fabs((*ipar).Value() - (*ipar).UpperLimit()) < par.Precision().Eps2()) {
- os <<"**";
- atHiLim = true;
- }
- os << " ||" << std::setw(8) << (*ipar).Error() << std::endl;
- } else
- os << " free ||" << std::setprecision(4) << std::setw(10) << (*ipar).Value() << " ||" << std::setw(8) << "no" << std::endl;
- } else {
- if((*ipar).Error() > 0.)
- os << " free ||" << std::setprecision(4) << std::setw(10) << (*ipar).Value() << " ||" << std::setw(8) << (*ipar).Error() << std::endl;
- else
- os << " free ||" << std::setprecision(4) << std::setw(10) << (*ipar).Value() << " ||" << std::setw(8) << "no" << std::endl;
-
- }
- }
- os << std::endl;
- if(atLoLim) os << "* Parameter is at Lower limit" << std::endl;
- if(atHiLim) os << "** Parameter is at Upper limit" << std::endl;
- os << std::endl;
-
- return os;
+ // print the MnUserParameter object
+ os << std::endl;
+
+ os << "# ext. |" << "| Name |" << "| type |" << "| Value |" << "| Error +/- " << std::endl;
+
+ os << std::endl;
+
+ bool atLoLim = false;
+ bool atHiLim = false;
+ for(std::vector<MinuitParameter>::const_iterator ipar = par.Parameters().begin(); ipar != par.Parameters().end(); ipar++) {
+ os << std::setw(4) << (*ipar).Number() << std::setw(5) << "||";
+ os << std::setw(10) << (*ipar).Name() << std::setw(3) << "||";
+ if((*ipar).IsConst()) {
+ os << " const ||" << std::setprecision(4) << std::setw(10) << (*ipar).Value() << " ||" << std::endl;
+ } else if((*ipar).IsFixed()) {
+ os << " fixed ||" << std::setprecision(4) << std::setw(10) << (*ipar).Value() << " ||" << std::endl;
+ } else if((*ipar).HasLimits()) {
+ if((*ipar).Error() > 0.) {
+ os << " limited ||" << std::setprecision(4) << std::setw(10) << (*ipar).Value();
+ if(fabs((*ipar).Value() - (*ipar).LowerLimit()) < par.Precision().Eps2()) {
+ os <<"*";
+ atLoLim = true;
+ }
+ if(fabs((*ipar).Value() - (*ipar).UpperLimit()) < par.Precision().Eps2()) {
+ os <<"**";
+ atHiLim = true;
+ }
+ os << " ||" << std::setw(8) << (*ipar).Error() << std::endl;
+ } else
+ os << " free ||" << std::setprecision(4) << std::setw(10) << (*ipar).Value() << " ||" << std::setw(8) << "no" << std::endl;
+ } else {
+ if((*ipar).Error() > 0.)
+ os << " free ||" << std::setprecision(4) << std::setw(10) << (*ipar).Value() << " ||" << std::setw(8) << (*ipar).Error() << std::endl;
+ else
+ os << " free ||" << std::setprecision(4) << std::setw(10) << (*ipar).Value() << " ||" << std::setw(8) << "no" << std::endl;
+
+ }
+ }
+ os << std::endl;
+ if(atLoLim) os << "* Parameter is at Lower limit" << std::endl;
+ if(atHiLim) os << "** Parameter is at Upper limit" << std::endl;
+ os << std::endl;
+
+ return os;
}
std::ostream& operator<<(std::ostream& os, const MnUserCovariance& matrix) {
-
- os << std::endl;
-
- os << "MnUserCovariance: " << std::endl;
-
- {
- os << std::endl;
- unsigned int n = matrix.Nrow();
- for (unsigned int i = 0; i < n; i++) {
- for (unsigned int j = 0; j < n; j++) {
- os.precision(6); os.width(13); os << matrix(i,j);
- }
+ // print the MnUserCovariance
+ os << std::endl;
+
+ os << "MnUserCovariance: " << std::endl;
+
+ {
os << std::endl;
- }
- }
-
+ unsigned int n = matrix.Nrow();
+ for (unsigned int i = 0; i < n; i++) {
+ for (unsigned int j = 0; j < n; j++) {
+ os.precision(6); os.width(13); os << matrix(i,j);
+ }
+ os << std::endl;
+ }
+ }
+
os << std::endl;
os << "MnUserCovariance Parameter correlations: " << std::endl;
-
+
{
- os << std::endl;
- unsigned int n = matrix.Nrow();
- for (unsigned int i = 0; i < n; i++) {
- double di = matrix(i,i);
- for (unsigned int j = 0; j < n; j++) {
- double dj = matrix(j,j);
- os.precision(6); os.width(13); os << matrix(i,j)/sqrt(fabs(di*dj));
- }
os << std::endl;
- }
- }
-
- return os;
+ unsigned int n = matrix.Nrow();
+ for (unsigned int i = 0; i < n; i++) {
+ double di = matrix(i,i);
+ for (unsigned int j = 0; j < n; j++) {
+ double dj = matrix(j,j);
+ os.precision(6); os.width(13); os << matrix(i,j)/sqrt(fabs(di*dj));
+ }
+ os << std::endl;
+ }
+ }
+
+ return os;
}
std::ostream& operator<<(std::ostream& os, const MnGlobalCorrelationCoeff& coeff) {
-
- os << std::endl;
-
- os << "MnGlobalCorrelationCoeff: " << std::endl;
-
+ // print the global correlation coefficient
+ os << std::endl;
+
+ os << "MnGlobalCorrelationCoeff: " << std::endl;
+
{
- os << std::endl;
- for (unsigned int i = 0; i < coeff.GlobalCC().size(); i++) {
- os.precision(6); os.width(13); os << coeff.GlobalCC()[i];
os << std::endl;
- }
- }
-
- return os;
+ for (unsigned int i = 0; i < coeff.GlobalCC().size(); i++) {
+ os.precision(6); os.width(13); os << coeff.GlobalCC()[i];
+ os << std::endl;
+ }
+ }
+
+ return os;
}
std::ostream& operator<<(std::ostream& os, const MnUserParameterState& state) {
-
- os << std::endl;
-
- if(!state.IsValid()) {
- os << std::endl;
- os <<"WARNING: MnUserParameterState is not valid."<<std::endl;
- os << std::endl;
- }
-
- os <<"# of function calls: "<<state.NFcn()<<std::endl;
- os <<"function Value: "<<state.Fval()<<std::endl;
- os <<"expected distance to the Minimum (edm): "<<state.Edm()<<std::endl;
- os <<"external parameters: "<<state.Parameters()<<std::endl;
- if(state.HasCovariance())
- os <<"covariance matrix: "<<state.Covariance()<<std::endl;
- if(state.HasGlobalCC())
- os <<"global correlation coefficients : "<<state.GlobalCC()<<std::endl;
-
- if(!state.IsValid())
- os <<"WARNING: MnUserParameterState is not valid."<<std::endl;
-
- os << std::endl;
-
- return os;
+ // print the MnUserParameterState
+ os << std::endl;
+
+ if(!state.IsValid()) {
+ os << std::endl;
+ os <<"WARNING: MnUserParameterState is not valid."<<std::endl;
+ os << std::endl;
+ }
+
+ os <<"# of function calls: "<<state.NFcn()<<std::endl;
+ os <<"function Value: "<<state.Fval()<<std::endl;
+ os <<"expected distance to the Minimum (edm): "<<state.Edm()<<std::endl;
+ os <<"external parameters: "<<state.Parameters()<<std::endl;
+ if(state.HasCovariance())
+ os <<"covariance matrix: "<<state.Covariance()<<std::endl;
+ if(state.HasGlobalCC())
+ os <<"global correlation coefficients : "<<state.GlobalCC()<<std::endl;
+
+ if(!state.IsValid())
+ os <<"WARNING: MnUserParameterState is not valid."<<std::endl;
+
+ os << std::endl;
+
+ return os;
}
std::ostream& operator<<(std::ostream& os, const FunctionMinimum& min) {
-
- os << std::endl;
- if(!min.IsValid()) {
- os << std::endl;
- os <<"WARNING: Minuit did not converge."<<std::endl;
- os << std::endl;
- } else {
- os << std::endl;
- os <<"Minuit did successfully converge."<<std::endl;
- os << std::endl;
- }
-
- os <<"# of function calls: "<<min.NFcn()<<std::endl;
- os <<"minimum function Value: "<<min.Fval()<<std::endl;
- os <<"minimum edm: "<<min.Edm()<<std::endl;
- os <<"minimum internal state vector: "<<min.Parameters().Vec()<<std::endl;
- if(min.HasValidCovariance())
- os <<"minimum internal covariance matrix: "<<min.Error().Matrix()<<std::endl;
-
- os << min.UserParameters() << std::endl;
- os << min.UserCovariance() << std::endl;
- os << min.UserState().GlobalCC() << std::endl;
-
- if(!min.IsValid())
- os <<"WARNING: FunctionMinimum is invalid."<<std::endl;
-
- os << std::endl;
-
- return os;
+ // print the FunctionMinimum
+ os << std::endl;
+ if(!min.IsValid()) {
+ os << std::endl;
+ os <<"WARNING: Minuit did not converge."<<std::endl;
+ os << std::endl;
+ } else {
+ os << std::endl;
+ os <<"Minuit did successfully converge."<<std::endl;
+ os << std::endl;
+ }
+
+ os <<"# of function calls: "<<min.NFcn()<<std::endl;
+ os <<"minimum function Value: "<<min.Fval()<<std::endl;
+ os <<"minimum edm: "<<min.Edm()<<std::endl;
+ os <<"minimum internal state vector: "<<min.Parameters().Vec()<<std::endl;
+ if(min.HasValidCovariance())
+ os <<"minimum internal covariance matrix: "<<min.Error().Matrix()<<std::endl;
+
+ os << min.UserParameters() << std::endl;
+ os << min.UserCovariance() << std::endl;
+ os << min.UserState().GlobalCC() << std::endl;
+
+ if(!min.IsValid())
+ os <<"WARNING: FunctionMinimum is invalid."<<std::endl;
+
+ os << std::endl;
+
+ return os;
}
std::ostream& operator<<(std::ostream& os, const MinimumState& min) {
-
- os << std::endl;
-
- os <<"minimum function Value: "<<min.Fval()<<std::endl;
- os <<"minimum edm: "<<min.Edm()<<std::endl;
- os <<"minimum internal state vector: "<<min.Vec()<<std::endl;
- os <<"minimum internal Gradient vector: "<<min.Gradient().Vec()<<std::endl;
- if(min.HasCovariance())
- os <<"minimum internal covariance matrix: "<<min.Error().Matrix()<<std::endl;
-
- os << std::endl;
-
- return os;
+
+ os << std::endl;
+
+ os <<"minimum function Value: "<<min.Fval()<<std::endl;
+ os <<"minimum edm: "<<min.Edm()<<std::endl;
+ os <<"minimum internal state vector: "<<min.Vec()<<std::endl;
+ os <<"minimum internal Gradient vector: "<<min.Gradient().Vec()<<std::endl;
+ if(min.HasCovariance())
+ os <<"minimum internal covariance matrix: "<<min.Error().Matrix()<<std::endl;
+
+ os << std::endl;
+
+ return os;
}
std::ostream& operator<<(std::ostream& os, const MnMachinePrecision& prec) {
-
- os << std::endl;
-
- os <<"current machine precision is set to "<<prec.Eps()<<std::endl;
-
- os << std::endl;
-
- return os;
+ // print the Precision
+ os << std::endl;
+
+ os <<"current machine precision is set to "<<prec.Eps()<<std::endl;
+
+ os << std::endl;
+
+ return os;
}
std::ostream& operator<<(std::ostream& os, const MinosError& me) {
-
- os << std::endl;
-
- os <<"Minos # of function calls: "<<me.NFcn()<<std::endl;
-
- if(!me.IsValid())
- os << "Minos Error is not valid." <<std::endl;
- if(!me.LowerValid())
- os << "lower Minos Error is not valid." <<std::endl;
- if(!me.UpperValid())
- os << "upper Minos Error is not valid." <<std::endl;
- if(me.AtLowerLimit())
- os << "Minos Error is Lower limit of Parameter "<<me.Parameter()<<"." <<std::endl;
- if(me.AtUpperLimit())
- os << "Minos Error is Upper limit of Parameter "<<me.Parameter()<<"." <<std::endl;
- if(me.AtLowerMaxFcn())
- os << "Minos number of function calls for Lower Error exhausted."<<std::endl;
- if(me.AtUpperMaxFcn())
- os << "Minos number of function calls for Upper Error exhausted."<<std::endl;
- if(me.LowerNewMin()) {
- os << "Minos found a new Minimum in negative direction."<<std::endl;
- os << me.LowerState() <<std::endl;
- }
- if(me.UpperNewMin()) {
- os << "Minos found a new Minimum in positive direction."<<std::endl;
- os << me.UpperState() <<std::endl;
- }
-
- os << "# ext. |" << "| Name |" << "| Value@min |" << "| negative |" << "| positive " << std::endl;
- os << std::setw(4) << me.Parameter() << std::setw(5) << "||";
- os << std::setw(10) << me.LowerState().Name(me.Parameter()) << std::setw(3) << "||";
- os << std::setprecision(4) << std::setw(10) << me.Min() << " ||" << std::setprecision(4) << std::setw(10) << me.Lower() << " ||" << std::setw(8) << me.Upper() << std::endl;
-
- os << std::endl;
-
- return os;
+ // print the Minos Error
+ os << std::endl;
+
+ os <<"Minos # of function calls: "<<me.NFcn()<<std::endl;
+
+ if(!me.IsValid())
+ os << "Minos Error is not valid." <<std::endl;
+ if(!me.LowerValid())
+ os << "lower Minos Error is not valid." <<std::endl;
+ if(!me.UpperValid())
+ os << "upper Minos Error is not valid." <<std::endl;
+ if(me.AtLowerLimit())
+ os << "Minos Error is Lower limit of Parameter "<<me.Parameter()<<"." <<std::endl;
+ if(me.AtUpperLimit())
+ os << "Minos Error is Upper limit of Parameter "<<me.Parameter()<<"." <<std::endl;
+ if(me.AtLowerMaxFcn())
+ os << "Minos number of function calls for Lower Error exhausted."<<std::endl;
+ if(me.AtUpperMaxFcn())
+ os << "Minos number of function calls for Upper Error exhausted."<<std::endl;
+ if(me.LowerNewMin()) {
+ os << "Minos found a new Minimum in negative direction."<<std::endl;
+ os << me.LowerState() <<std::endl;
+ }
+ if(me.UpperNewMin()) {
+ os << "Minos found a new Minimum in positive direction."<<std::endl;
+ os << me.UpperState() <<std::endl;
+ }
+
+ os << "# ext. |" << "| Name |" << "| Value@min |" << "| negative |" << "| positive " << std::endl;
+ os << std::setw(4) << me.Parameter() << std::setw(5) << "||";
+ os << std::setw(10) << me.LowerState().Name(me.Parameter()) << std::setw(3) << "||";
+ os << std::setprecision(4) << std::setw(10) << me.Min() << " ||" << std::setprecision(4) << std::setw(10) << me.Lower() << " ||" << std::setw(8) << me.Upper() << std::endl;
+
+ os << std::endl;
+
+ return os;
}
std::ostream& operator<<(std::ostream& os, const ContoursError& ce) {
-
- os << std::endl;
- os <<"Contours # of function calls: "<<ce.NFcn()<<std::endl;
- os << "MinosError in x: "<<std::endl;
- os << ce.XMinosError() << std::endl;
- os << "MinosError in y: "<<std::endl;
- os << ce.YMinosError() << std::endl;
- MnPlot plot;
- plot(ce.XMin(), ce.YMin(), ce());
- for(std::vector<std::pair<double,double> >::const_iterator ipar = ce().begin(); ipar != ce().end(); ipar++) {
- os << ipar - ce().begin() <<" "<< (*ipar).first <<" "<< (*ipar).second <<std::endl;
- }
- os << std::endl;
-
- return os;
+ // print the ContoursError
+ os << std::endl;
+ os <<"Contours # of function calls: "<<ce.NFcn()<<std::endl;
+ os << "MinosError in x: "<<std::endl;
+ os << ce.XMinosError() << std::endl;
+ os << "MinosError in y: "<<std::endl;
+ os << ce.YMinosError() << std::endl;
+ MnPlot plot;
+ plot(ce.XMin(), ce.YMin(), ce());
+ for(std::vector<std::pair<double,double> >::const_iterator ipar = ce().begin(); ipar != ce().end(); ipar++) {
+ os << ipar - ce().begin() <<" "<< (*ipar).first <<" "<< (*ipar).second <<std::endl;
+ }
+ os << std::endl;
+
+ return os;
}
- } // namespace Minuit2
+ } // namespace Minuit2
} // namespace ROOT
diff --git a/minuit2/src/MnScan.cxx b/minuit2/src/MnScan.cxx
index abbf6b2f60b..183bd013e00 100644
--- a/minuit2/src/MnScan.cxx
+++ b/minuit2/src/MnScan.cxx
@@ -1,4 +1,4 @@
-// @(#)root/minuit2:$Name: $:$Id: MnScan.cpp,v 1.3.6.4 2005/11/29 11:08:35 moneta Exp $
+// @(#)root/minuit2:$Name: $:$Id: MnScan.cxx,v 1.1 2005/11/29 14:43:31 moneta Exp $
// Authors: M. Winkler, F. James, L. Moneta, A. Zsenei 2003-2005
/**********************************************************************
@@ -16,19 +16,19 @@ namespace ROOT {
std::vector<std::pair<double, double> > MnScan::Scan(unsigned int par, unsigned int maxsteps, double low, double high) {
-
- MnParameterScan Scan(fFCN, fState.Parameters());
- double amin = Scan.Fval();
-
- std::vector<std::pair<double, double> > result = Scan(par, maxsteps, low, high);
- if(Scan.Fval() < amin) {
- fState.SetValue(par, Scan.Parameters().Value(par));
- amin = Scan.Fval();
- }
+ // perform a scan of the function in the parameter par
+ MnParameterScan Scan(fFCN, fState.Parameters());
+ double amin = Scan.Fval();
- return result;
+ std::vector<std::pair<double, double> > result = Scan(par, maxsteps, low, high);
+ if(Scan.Fval() < amin) {
+ fState.SetValue(par, Scan.Parameters().Value(par));
+ amin = Scan.Fval();
+ }
+
+ return result;
}
- } // namespace Minuit2
+ } // namespace Minuit2
} // namespace ROOT
diff --git a/minuit2/src/MnSeedGenerator.cxx b/minuit2/src/MnSeedGenerator.cxx
index cde57502ffd..8a22ea3f17f 100644
--- a/minuit2/src/MnSeedGenerator.cxx
+++ b/minuit2/src/MnSeedGenerator.cxx
@@ -1,4 +1,4 @@
-// @(#)root/minuit2:$Name: $:$Id: MnSeedGenerator.cpp,v 1.12.2.4 2005/11/29 11:08:35 moneta Exp $
+// @(#)root/minuit2:$Name: $:$Id: MnSeedGenerator.cxx,v 1.1 2005/11/29 14:43:31 moneta Exp $
// Authors: M. Winkler, F. James, L. Moneta, A. Zsenei 2003-2005
/**********************************************************************
@@ -40,106 +40,106 @@ namespace ROOT {
MinimumSeed MnSeedGenerator::operator()(const MnFcn& fcn, const GradientCalculator& gc, const MnUserParameterState& st, const MnStrategy& stra) const {
-
- unsigned int n = st.VariableParameters();
- const MnMachinePrecision& prec = st.Precision();
-
- // initial starting values
- MnAlgebraicVector x(n);
- for(unsigned int i = 0; i < n; i++) x(i) = st.IntParameters()[i];
- double fcnmin = fcn(x);
- MinimumParameters pa(x, fcnmin);
- FunctionGradient dgrad = gc(pa);
- MnAlgebraicSymMatrix mat(n);
- double dcovar = 1.;
- if(st.HasCovariance()) {
- for(unsigned int i = 0; i < n; i++)
- for(unsigned int j = i; j < n; j++) mat(i,j) = st.IntCovariance()(i,j);
- dcovar = 0.;
- } else {
- for(unsigned int i = 0; i < n; i++)
- mat(i,i) = (fabs(dgrad.G2()(i)) > prec.Eps2() ? 1./dgrad.G2()(i) : 1.);
- }
- MinimumError err(mat, dcovar);
- double edm = VariableMetricEDMEstimator().Estimate(dgrad, err);
- MinimumState state(pa, err, dgrad, edm, fcn.NumOfCalls());
-
- NegativeG2LineSearch ng2ls;
- if(ng2ls.HasNegativeG2(dgrad, prec)) {
- state = ng2ls(fcn, state, gc, prec);
- }
-
- if(stra.Strategy() == 2 && !st.HasCovariance()) {
- //calculate full 2nd derivative
- MinimumState tmp = MnHesse(stra)(fcn, state, st.Trafo());
- return MinimumSeed(tmp, st.Trafo());
- }
-
- return MinimumSeed(state, st.Trafo());
+ // find seed (initial minimization point) using the calculated gradient
+ unsigned int n = st.VariableParameters();
+ const MnMachinePrecision& prec = st.Precision();
+
+ // initial starting values
+ MnAlgebraicVector x(n);
+ for(unsigned int i = 0; i < n; i++) x(i) = st.IntParameters()[i];
+ double fcnmin = fcn(x);
+ MinimumParameters pa(x, fcnmin);
+ FunctionGradient dgrad = gc(pa);
+ MnAlgebraicSymMatrix mat(n);
+ double dcovar = 1.;
+ if(st.HasCovariance()) {
+ for(unsigned int i = 0; i < n; i++)
+ for(unsigned int j = i; j < n; j++) mat(i,j) = st.IntCovariance()(i,j);
+ dcovar = 0.;
+ } else {
+ for(unsigned int i = 0; i < n; i++)
+ mat(i,i) = (fabs(dgrad.G2()(i)) > prec.Eps2() ? 1./dgrad.G2()(i) : 1.);
+ }
+ MinimumError err(mat, dcovar);
+ double edm = VariableMetricEDMEstimator().Estimate(dgrad, err);
+ MinimumState state(pa, err, dgrad, edm, fcn.NumOfCalls());
+
+ NegativeG2LineSearch ng2ls;
+ if(ng2ls.HasNegativeG2(dgrad, prec)) {
+ state = ng2ls(fcn, state, gc, prec);
+ }
+
+ if(stra.Strategy() == 2 && !st.HasCovariance()) {
+ //calculate full 2nd derivative
+ MinimumState tmp = MnHesse(stra)(fcn, state, st.Trafo());
+ return MinimumSeed(tmp, st.Trafo());
+ }
+
+ return MinimumSeed(state, st.Trafo());
}
MinimumSeed MnSeedGenerator::operator()(const MnFcn& fcn, const AnalyticalGradientCalculator& gc, const MnUserParameterState& st, const MnStrategy& stra) const {
-
- unsigned int n = st.VariableParameters();
- const MnMachinePrecision& prec = st.Precision();
-
- // initial starting values
- MnAlgebraicVector x(n);
- for(unsigned int i = 0; i < n; i++) x(i) = st.IntParameters()[i];
- double fcnmin = fcn(x);
- MinimumParameters pa(x, fcnmin);
-
- InitialGradientCalculator igc(fcn, st.Trafo(), stra);
- FunctionGradient tmp = igc(pa);
- FunctionGradient grd = gc(pa);
- FunctionGradient dgrad(grd.Grad(), tmp.G2(), tmp.Gstep());
-
- if(gc.CheckGradient()) {
- bool good = true;
- HessianGradientCalculator hgc(fcn, st.Trafo(), MnStrategy(2));
- std::pair<FunctionGradient, MnAlgebraicVector> hgrd = hgc.DeltaGradient(pa, dgrad);
- for(unsigned int i = 0; i < n; i++) {
- if(fabs(hgrd.first.Grad()(i) - grd.Grad()(i)) > hgrd.second(i)) {
- std::cout<<"gradient discrepancy of external Parameter "<<st.Trafo().ExtOfInt(i)<<" (internal Parameter "<<i<<") too large."<<std::endl;
- good = false;
+ // find seed (initial point for minimization) using analytical gradient
+ unsigned int n = st.VariableParameters();
+ const MnMachinePrecision& prec = st.Precision();
+
+ // initial starting values
+ MnAlgebraicVector x(n);
+ for(unsigned int i = 0; i < n; i++) x(i) = st.IntParameters()[i];
+ double fcnmin = fcn(x);
+ MinimumParameters pa(x, fcnmin);
+
+ InitialGradientCalculator igc(fcn, st.Trafo(), stra);
+ FunctionGradient tmp = igc(pa);
+ FunctionGradient grd = gc(pa);
+ FunctionGradient dgrad(grd.Grad(), tmp.G2(), tmp.Gstep());
+
+ if(gc.CheckGradient()) {
+ bool good = true;
+ HessianGradientCalculator hgc(fcn, st.Trafo(), MnStrategy(2));
+ std::pair<FunctionGradient, MnAlgebraicVector> hgrd = hgc.DeltaGradient(pa, dgrad);
+ for(unsigned int i = 0; i < n; i++) {
+ if(fabs(hgrd.first.Grad()(i) - grd.Grad()(i)) > hgrd.second(i)) {
+ std::cout<<"gradient discrepancy of external Parameter "<<st.Trafo().ExtOfInt(i)<<" (internal Parameter "<<i<<") too large."<<std::endl;
+ good = false;
+ }
}
- }
- if(!good) {
- std::cout<<"Minuit does not accept user specified Gradient. To force acceptance, override 'virtual bool CheckGradient() const' of FCNGradientBase.h in the derived class."<<std::endl;
- assert(good);
- }
- }
-
- MnAlgebraicSymMatrix mat(n);
- double dcovar = 1.;
- if(st.HasCovariance()) {
- for(unsigned int i = 0; i < n; i++)
- for(unsigned int j = i; j < n; j++) mat(i,j) = st.IntCovariance()(i,j);
- dcovar = 0.;
- } else {
- for(unsigned int i = 0; i < n; i++)
- mat(i,i) = (fabs(dgrad.G2()(i)) > prec.Eps2() ? 1./dgrad.G2()(i) : 1.);
- }
- MinimumError err(mat, dcovar);
- double edm = VariableMetricEDMEstimator().Estimate(dgrad, err);
- MinimumState state(pa, err, dgrad, edm, fcn.NumOfCalls());
-
- NegativeG2LineSearch ng2ls;
- if(ng2ls.HasNegativeG2(dgrad, prec)) {
- Numerical2PGradientCalculator ngc(fcn, st.Trafo(), stra);
- state = ng2ls(fcn, state, ngc, prec);
- }
-
- if(stra.Strategy() == 2 && !st.HasCovariance()) {
- //calculate full 2nd derivative
- MinimumState tmp = MnHesse(stra)(fcn, state, st.Trafo());
- return MinimumSeed(tmp, st.Trafo());
- }
-
- return MinimumSeed(state, st.Trafo());
+ if(!good) {
+ std::cout<<"Minuit does not accept user specified Gradient. To force acceptance, override 'virtual bool CheckGradient() const' of FCNGradientBase.h in the derived class."<<std::endl;
+ assert(good);
+ }
+ }
+
+ MnAlgebraicSymMatrix mat(n);
+ double dcovar = 1.;
+ if(st.HasCovariance()) {
+ for(unsigned int i = 0; i < n; i++)
+ for(unsigned int j = i; j < n; j++) mat(i,j) = st.IntCovariance()(i,j);
+ dcovar = 0.;
+ } else {
+ for(unsigned int i = 0; i < n; i++)
+ mat(i,i) = (fabs(dgrad.G2()(i)) > prec.Eps2() ? 1./dgrad.G2()(i) : 1.);
+ }
+ MinimumError err(mat, dcovar);
+ double edm = VariableMetricEDMEstimator().Estimate(dgrad, err);
+ MinimumState state(pa, err, dgrad, edm, fcn.NumOfCalls());
+
+ NegativeG2LineSearch ng2ls;
+ if(ng2ls.HasNegativeG2(dgrad, prec)) {
+ Numerical2PGradientCalculator ngc(fcn, st.Trafo(), stra);
+ state = ng2ls(fcn, state, ngc, prec);
+ }
+
+ if(stra.Strategy() == 2 && !st.HasCovariance()) {
+ //calculate full 2nd derivative
+ MinimumState tmp = MnHesse(stra)(fcn, state, st.Trafo());
+ return MinimumSeed(tmp, st.Trafo());
+ }
+
+ return MinimumSeed(state, st.Trafo());
}
- } // namespace Minuit2
+ } // namespace Minuit2
} // namespace ROOT
diff --git a/minuit2/src/MnStrategy.cxx b/minuit2/src/MnStrategy.cxx
index 83c1ea6440e..c5790e29cc6 100644
--- a/minuit2/src/MnStrategy.cxx
+++ b/minuit2/src/MnStrategy.cxx
@@ -1,4 +1,4 @@
-// @(#)root/minuit2:$Name: $:$Id: MnStrategy.cpp,v 1.3.6.3 2005/11/29 11:08:35 moneta Exp $
+// @(#)root/minuit2:$Name: $:$Id: MnStrategy.cxx,v 1.1 2005/11/29 14:43:31 moneta Exp $
// Authors: M. Winkler, F. James, L. Moneta, A. Zsenei 2003-2005
/**********************************************************************
@@ -14,51 +14,56 @@ namespace ROOT {
namespace Minuit2 {
-//default strategy
+
MnStrategy::MnStrategy() {
- SetMediumStrategy();
+ //default strategy
+ SetMediumStrategy();
}
-//user defined strategy (0, 1, >=2)
+
MnStrategy::MnStrategy(unsigned int stra) {
- if(stra == 0) SetLowStrategy();
- else if(stra == 1) SetMediumStrategy();
- else SetHighStrategy();
+ //user defined strategy (0, 1, >=2)
+ if(stra == 0) SetLowStrategy();
+ else if(stra == 1) SetMediumStrategy();
+ else SetHighStrategy();
}
void MnStrategy::SetLowStrategy() {
- fStrategy = 0;
- SetGradientNCycles(2);
- SetGradientStepTolerance(0.5);
- SetGradientTolerance(0.1);
- SetHessianNCycles(3);
- SetHessianStepTolerance(0.5);
- SetHessianG2Tolerance(0.1);
- SetHessianGradientNCycles(1);
+ // set low strategy (0) values
+ fStrategy = 0;
+ SetGradientNCycles(2);
+ SetGradientStepTolerance(0.5);
+ SetGradientTolerance(0.1);
+ SetHessianNCycles(3);
+ SetHessianStepTolerance(0.5);
+ SetHessianG2Tolerance(0.1);
+ SetHessianGradientNCycles(1);
}
void MnStrategy::SetMediumStrategy() {
- fStrategy = 1;
- SetGradientNCycles(3);
- SetGradientStepTolerance(0.3);
- SetGradientTolerance(0.05);
- SetHessianNCycles(5);
- SetHessianStepTolerance(0.3);
- SetHessianG2Tolerance(0.05);
- SetHessianGradientNCycles(2);
+ // set minimum strategy (1) the default
+ fStrategy = 1;
+ SetGradientNCycles(3);
+ SetGradientStepTolerance(0.3);
+ SetGradientTolerance(0.05);
+ SetHessianNCycles(5);
+ SetHessianStepTolerance(0.3);
+ SetHessianG2Tolerance(0.05);
+ SetHessianGradientNCycles(2);
}
void MnStrategy::SetHighStrategy() {
- fStrategy = 2;
- SetGradientNCycles(5);
- SetGradientStepTolerance(0.1);
- SetGradientTolerance(0.02);
- SetHessianNCycles(7);
- SetHessianStepTolerance(0.1);
- SetHessianG2Tolerance(0.02);
- SetHessianGradientNCycles(6);
+ // set high strategy (2)
+ fStrategy = 2;
+ SetGradientNCycles(5);
+ SetGradientStepTolerance(0.1);
+ SetGradientTolerance(0.02);
+ SetHessianNCycles(7);
+ SetHessianStepTolerance(0.1);
+ SetHessianG2Tolerance(0.02);
+ SetHessianGradientNCycles(6);
}
- } // namespace Minuit2
+ } // namespace Minuit2
} // namespace ROOT
diff --git a/minuit2/src/MnTiny.cxx b/minuit2/src/MnTiny.cxx
index be0a32ce3b9..819cac37157 100644
--- a/minuit2/src/MnTiny.cxx
+++ b/minuit2/src/MnTiny.cxx
@@ -1,4 +1,4 @@
-// @(#)root/minuit2:$Name: $:$Id: MnTiny.cpp,v 1.1.6.3 2005/11/29 11:08:35 moneta Exp $
+// @(#)root/minuit2:$Name: $:$Id: MnTiny.cxx,v 1.1 2005/11/29 14:43:31 moneta Exp $
// Authors: M. Winkler, F. James, L. Moneta, A. Zsenei 2003-2005
/**********************************************************************
@@ -17,10 +17,11 @@ namespace ROOT {
double MnTiny::One() const {return fOne;}
double MnTiny::operator()(double epsp1) const {
- double result = epsp1 - One();
- return result;
+ // evaluate minimal diference between two floating points
+ double result = epsp1 - One();
+ return result;
}
- } // namespace Minuit2
+ } // namespace Minuit2
} // namespace ROOT
diff --git a/minuit2/src/MnUserFcn.cxx b/minuit2/src/MnUserFcn.cxx
index 744c07dc87b..e676c3ae2ed 100644
--- a/minuit2/src/MnUserFcn.cxx
+++ b/minuit2/src/MnUserFcn.cxx
@@ -1,4 +1,4 @@
-// @(#)root/minuit2:$Name: $:$Id: MnUserFcn.cpp,v 1.2.6.3 2005/11/29 11:08:35 moneta Exp $
+// @(#)root/minuit2:$Name: $:$Id: MnUserFcn.cxx,v 1.1 2005/11/29 14:43:31 moneta Exp $
// Authors: M. Winkler, F. James, L. Moneta, A. Zsenei 2003-2005
/**********************************************************************
@@ -17,11 +17,11 @@ namespace ROOT {
double MnUserFcn::operator()(const MnAlgebraicVector& v) const {
-
- fNumCall++;
- return Fcn()( fTransform(v) );
+ // call Fcn function transforming from internal parameter vector to external std::vector
+ fNumCall++;
+ return Fcn()( fTransform(v) );
}
- } // namespace Minuit2
+ } // namespace Minuit2
} // namespace ROOT
diff --git a/minuit2/src/mnlsame.cxx b/minuit2/src/mnlsame.cxx
index 83f848f7a47..d4e7dd117a8 100644
--- a/minuit2/src/mnlsame.cxx
+++ b/minuit2/src/mnlsame.cxx
@@ -1,4 +1,4 @@
-// @(#)root/minuit2:$Name: $:$Id: mnlsame.cpp,v 1.1.6.2 2005/11/29 11:08:35 moneta Exp $
+// @(#)root/minuit2:$Name: $:$Id: mnlsame.cxx,v 1.1 2005/11/29 14:43:31 moneta Exp $
// Authors: M. Winkler, F. James, L. Moneta, A. Zsenei 2003-2005
/**********************************************************************
@@ -20,51 +20,51 @@ namespace ROOT {
bool mnlsame(const char* ca, const char* cb) {
- /* System generated locals */
- bool ret_val = false;
-
- /* Local variables */
-// integer inta, intb, zcode;
-
-
-/* -- LAPACK auxiliary routine (version 2.0) -- */
-/* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., */
-/* Courant Institute, Argonne National Lab, and Rice University */
-/* January 31, 1994 */
-
-/* .. Scalar Arguments .. */
-/* .. */
-
-/* Purpose */
-/* ======= */
-
-/* LSAME returns .TRUE. if CA is the same letter as CB regardless of */
-/* case. */
-
-/* Arguments */
-/* ========= */
-
-/* CA (input) CHARACTER*1 */
-/* CB (input) CHARACTER*1 */
-/* CA and CB specify the single characters to be compared. */
-
-/* ===================================================================== */
-
-/* .. Intrinsic Functions .. */
-/* .. */
-/* .. Local Scalars .. */
-/* .. */
-/* .. Executable Statements .. */
-
-/* Test if the characters are equal */
-
- int comp = strcmp(ca, cb);
- if(comp == 0) ret_val = true;
-
- return ret_val;
+ /* System generated locals */
+ bool ret_val = false;
+
+ /* Local variables */
+ // integer inta, intb, zcode;
+
+
+ /* -- LAPACK auxiliary routine (version 2.0) -- */
+ /* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., */
+ /* Courant Institute, Argonne National Lab, and Rice University */
+ /* January 31, 1994 */
+
+ /* .. Scalar Arguments .. */
+ /* .. */
+
+ /* Purpose */
+ /* ======= */
+
+ /* LSAME returns .TRUE. if CA is the same letter as CB regardless of */
+ /* case. */
+
+ /* Arguments */
+ /* ========= */
+
+ /* CA (input) CHARACTER*1 */
+ /* CB (input) CHARACTER*1 */
+ /* CA and CB specify the single characters to be compared. */
+
+ /* ===================================================================== */
+
+ /* .. Intrinsic Functions .. */
+ /* .. */
+ /* .. Local Scalars .. */
+ /* .. */
+ /* .. Executable Statements .. */
+
+ /* Test if the characters are equal */
+
+ int comp = strcmp(ca, cb);
+ if(comp == 0) ret_val = true;
+
+ return ret_val;
} /* lsame_ */
- } // namespace Minuit2
+ } // namespace Minuit2
} // namespace ROOT
diff --git a/minuit2/src/mntplot.cxx b/minuit2/src/mntplot.cxx
index 37f0b3767ff..ddc596d7ee9 100644
--- a/minuit2/src/mntplot.cxx
+++ b/minuit2/src/mntplot.cxx
@@ -1,4 +1,4 @@
-// @(#)root/minuit2:$Name: $:$Id: mntplot.cpp,v 1.1.6.3 2005/11/29 11:08:35 moneta Exp $
+// @(#)root/minuit2:$Name: $:$Id: mntplot.cxx,v 1.1 2005/11/29 14:43:31 moneta Exp $
// Authors: M. Winkler, F. James, L. Moneta, A. Zsenei 2003-2005
/**********************************************************************
@@ -37,164 +37,164 @@ void mnplot(double* xpt, double* ypt, char* chpt, int nxypt, int npagwd, int npa
// char cdot[] = ".";
// char cslash[] = "/";
- /* Local variables */
- double xmin, ymin, xmax, ymax, savx, savy, yprt;
- double bwidx, bwidy, xbest, ybest, ax, ay, bx, by;
- double xvalus[12], any, dxx, dyy;
- int iten, i, j, k, maxnx, maxny, iquit, ni, linodd;
- int nxbest, nybest, km1, ibk, isp1, nx, ny, ks, ix;
- char ctemp[120];
- bool overpr;
- char cline[120];
- char chsav, chbest;
-
- /* Function Body */
- //*-* Computing MIN
- maxnx = npagwd-20 < 100 ? npagwd-20 : 100;
- if (maxnx < 10) maxnx = 10;
- maxny = npagln;
- if (maxny < 10) maxny = 10;
- if (nxypt <= 1) return;
- xbest = xpt[0];
- ybest = ypt[0];
- chbest = chpt[0];
- //*-*- order the points by decreasing y
- km1 = nxypt - 1;
- for (i = 1; i <= km1; ++i) {
- iquit = 0;
- ni = nxypt - i;
- for (j = 1; j <= ni; ++j) {
- if (ypt[j-1] > ypt[j]) continue;
- savx = xpt[j-1];
- xpt[j-1] = xpt[j];
- xpt[j] = savx;
- savy = ypt[j-1];
- ypt[j-1] = ypt[j];
- ypt[j] = savy;
- chsav = chpt[j-1];
- chpt[j-1]= chpt[j];
- chpt[j] = chsav;
- iquit = 1;
- }
- if (iquit == 0) break;
- }
- //*-*- find extreme values
- xmax = xpt[0];
- xmin = xmax;
- for (i = 1; i <= nxypt; ++i) {
- if (xpt[i-1] > xmax) xmax = xpt[i-1];
- if (xpt[i-1] < xmin) xmin = xpt[i-1];
- }
- dxx = (xmax - xmin)*.001;
- xmax += dxx;
- xmin -= dxx;
- mnbins(xmin, xmax, maxnx, xmin, xmax, nx, bwidx);
- ymax = ypt[0];
- ymin = ypt[nxypt-1];
- if (ymax == ymin) ymax = ymin + 1;
- dyy = (ymax - ymin)*.001;
- ymax += dyy;
- ymin -= dyy;
- mnbins(ymin, ymax, maxny, ymin, ymax, ny, bwidy);
- any = (double) ny;
- //*-*- if first point is blank, it is an 'origin'
- if (chbest == ' ') goto L50;
- xbest = (xmax + xmin)*.5;
- ybest = (ymax + ymin)*.5;
- L50:
- //*-*- find Scale constants
- ax = 1 / bwidx;
- ay = 1 / bwidy;
- bx = -ax*xmin + 2;
- by = -ay*ymin - 2;
- //*-*- convert points to grid positions
- for (i = 1; i <= nxypt; ++i) {
- xpt[i-1] = ax*xpt[i-1] + bx;
- ypt[i-1] = any - ay*ypt[i-1] - by;
- }
- nxbest = int((ax*xbest + bx));
- nybest = int((any - ay*ybest - by));
- //*-*- print the points
- ny += 2;
- nx += 2;
- isp1 = 1;
- linodd = 1;
- overpr = false;
- for (i = 1; i <= ny; ++i) {
- for (ibk = 1; ibk <= nx; ++ibk) { cline[ibk-1] = ' '; }
- cline[nx] = '\0';
- cline[nx+1] = '\0';
- cline[0] = '.';
- cline[nx-1] = '.';
- cline[nxbest-1] = '.';
- if (i != 1 && i != nybest && i != ny) goto L320;
- for (j = 1; j <= nx; ++j) { cline[j-1] = '.'; }
- L320:
- yprt = ymax - double(i-1)*bwidy;
- if (isp1 > nxypt) goto L350;
- //*-*- find the points to be plotted on this line
- for (k = isp1; k <= nxypt; ++k) {
- ks = int(ypt[k-1]);
- if (ks > i) goto L345;
- ix = int(xpt[k-1]);
- if (cline[ix-1] == '.') goto L340;
- if (cline[ix-1] == ' ') goto L340;
- if (cline[ix-1] == chpt[k-1]) continue;
- overpr = true;
- //*-*- OVERPR is true if one or more positions contains more than
- //*-*- one point
- cline[ix-1] = '&';
- continue;
- L340:
- cline[ix-1] = chpt[k-1];
- }
- isp1 = nxypt + 1;
- goto L350;
- L345:
- isp1 = k;
- L350:
- if (linodd == 1 || i == ny) goto L380;
- linodd = 1;
- memcpy(ctemp, cline, 120);
- printf(" %s",(const char*)ctemp);
- goto L400;
- L380:
-// ctemp = cline;
- memcpy(ctemp, cline, 120);
- printf(" %14.7g ..%s",yprt,(const char*)ctemp);
- linodd = 0;
- L400:
- printf("\n");
- }
- //*-*- print labels on x-axis every ten columns
- for (ibk = 1; ibk <= nx; ++ibk) {
- cline[ibk-1] = ' ';
- if (ibk % 10 == 1) cline[ibk-1] = '/';
- }
- printf(" %s",cline);
- printf("\n");
-
- for (ibk = 1; ibk <= 12; ++ibk) {
- xvalus[ibk-1] = xmin + double(ibk-1)*10*bwidx;
- }
- printf(" ");
- iten = (nx + 9) / 10;
- for (ibk = 1; ibk <= iten; ++ibk) {
- printf(" %9.4g", xvalus[ibk-1]);
- }
- printf("\n");
-
- if (overpr) {
- char chmess[] = " Overprint character is &";
- printf(" ONE COLUMN=%13.7g%s",bwidx,(const char*)chmess);
- } else {
- char chmess[] = " ";
- printf(" ONE COLUMN=%13.7g%s",bwidx,(const char*)chmess);
- }
- printf("\n");
-
+ /* Local variables */
+ double xmin, ymin, xmax, ymax, savx, savy, yprt;
+ double bwidx, bwidy, xbest, ybest, ax, ay, bx, by;
+ double xvalus[12], any, dxx, dyy;
+ int iten, i, j, k, maxnx, maxny, iquit, ni, linodd;
+ int nxbest, nybest, km1, ibk, isp1, nx, ny, ks, ix;
+ char ctemp[120];
+ bool overpr;
+ char cline[120];
+ char chsav, chbest;
+
+ /* Function Body */
+ //*-* Computing MIN
+ maxnx = npagwd-20 < 100 ? npagwd-20 : 100;
+ if (maxnx < 10) maxnx = 10;
+ maxny = npagln;
+ if (maxny < 10) maxny = 10;
+ if (nxypt <= 1) return;
+ xbest = xpt[0];
+ ybest = ypt[0];
+ chbest = chpt[0];
+ //*-*- order the points by decreasing y
+ km1 = nxypt - 1;
+ for (i = 1; i <= km1; ++i) {
+ iquit = 0;
+ ni = nxypt - i;
+ for (j = 1; j <= ni; ++j) {
+ if (ypt[j-1] > ypt[j]) continue;
+ savx = xpt[j-1];
+ xpt[j-1] = xpt[j];
+ xpt[j] = savx;
+ savy = ypt[j-1];
+ ypt[j-1] = ypt[j];
+ ypt[j] = savy;
+ chsav = chpt[j-1];
+ chpt[j-1]= chpt[j];
+ chpt[j] = chsav;
+ iquit = 1;
+ }
+ if (iquit == 0) break;
+ }
+ //*-*- find extreme values
+ xmax = xpt[0];
+ xmin = xmax;
+ for (i = 1; i <= nxypt; ++i) {
+ if (xpt[i-1] > xmax) xmax = xpt[i-1];
+ if (xpt[i-1] < xmin) xmin = xpt[i-1];
+ }
+ dxx = (xmax - xmin)*.001;
+ xmax += dxx;
+ xmin -= dxx;
+ mnbins(xmin, xmax, maxnx, xmin, xmax, nx, bwidx);
+ ymax = ypt[0];
+ ymin = ypt[nxypt-1];
+ if (ymax == ymin) ymax = ymin + 1;
+ dyy = (ymax - ymin)*.001;
+ ymax += dyy;
+ ymin -= dyy;
+ mnbins(ymin, ymax, maxny, ymin, ymax, ny, bwidy);
+ any = (double) ny;
+ //*-*- if first point is blank, it is an 'origin'
+ if (chbest == ' ') goto L50;
+ xbest = (xmax + xmin)*.5;
+ ybest = (ymax + ymin)*.5;
+L50:
+ //*-*- find Scale constants
+ ax = 1 / bwidx;
+ ay = 1 / bwidy;
+ bx = -ax*xmin + 2;
+ by = -ay*ymin - 2;
+ //*-*- convert points to grid positions
+ for (i = 1; i <= nxypt; ++i) {
+ xpt[i-1] = ax*xpt[i-1] + bx;
+ ypt[i-1] = any - ay*ypt[i-1] - by;
+ }
+ nxbest = int((ax*xbest + bx));
+ nybest = int((any - ay*ybest - by));
+ //*-*- print the points
+ ny += 2;
+ nx += 2;
+ isp1 = 1;
+ linodd = 1;
+ overpr = false;
+ for (i = 1; i <= ny; ++i) {
+ for (ibk = 1; ibk <= nx; ++ibk) { cline[ibk-1] = ' '; }
+ cline[nx] = '\0';
+ cline[nx+1] = '\0';
+ cline[0] = '.';
+ cline[nx-1] = '.';
+ cline[nxbest-1] = '.';
+ if (i != 1 && i != nybest && i != ny) goto L320;
+ for (j = 1; j <= nx; ++j) { cline[j-1] = '.'; }
+L320:
+ yprt = ymax - double(i-1)*bwidy;
+ if (isp1 > nxypt) goto L350;
+ //*-*- find the points to be plotted on this line
+ for (k = isp1; k <= nxypt; ++k) {
+ ks = int(ypt[k-1]);
+ if (ks > i) goto L345;
+ ix = int(xpt[k-1]);
+ if (cline[ix-1] == '.') goto L340;
+ if (cline[ix-1] == ' ') goto L340;
+ if (cline[ix-1] == chpt[k-1]) continue;
+ overpr = true;
+ //*-*- OVERPR is true if one or more positions contains more than
+ //*-*- one point
+ cline[ix-1] = '&';
+ continue;
+L340:
+ cline[ix-1] = chpt[k-1];
+ }
+ isp1 = nxypt + 1;
+ goto L350;
+L345:
+ isp1 = k;
+L350:
+ if (linodd == 1 || i == ny) goto L380;
+ linodd = 1;
+ memcpy(ctemp, cline, 120);
+ printf(" %s",(const char*)ctemp);
+ goto L400;
+L380:
+ // ctemp = cline;
+ memcpy(ctemp, cline, 120);
+ printf(" %14.7g ..%s",yprt,(const char*)ctemp);
+ linodd = 0;
+L400:
+ printf("\n");
+ }
+ //*-*- print labels on x-axis every ten columns
+ for (ibk = 1; ibk <= nx; ++ibk) {
+ cline[ibk-1] = ' ';
+ if (ibk % 10 == 1) cline[ibk-1] = '/';
+ }
+ printf(" %s",cline);
+ printf("\n");
+
+ for (ibk = 1; ibk <= 12; ++ibk) {
+ xvalus[ibk-1] = xmin + double(ibk-1)*10*bwidx;
+ }
+ printf(" ");
+ iten = (nx + 9) / 10;
+ for (ibk = 1; ibk <= iten; ++ibk) {
+ printf(" %9.4g", xvalus[ibk-1]);
+ }
+ printf("\n");
+
+ if (overpr) {
+ char chmess[] = " Overprint character is &";
+ printf(" ONE COLUMN=%13.7g%s",bwidx,(const char*)chmess);
+ } else {
+ char chmess[] = " ";
+ printf(" ONE COLUMN=%13.7g%s",bwidx,(const char*)chmess);
+ }
+ printf("\n");
+
} /* mnplot_ */
- } // namespace Minuit2
+ } // namespace Minuit2
} // namespace ROOT
--
GitLab