Skip to content

Commit

Permalink
NumericalDerivatorMinuit2 now also uses exact same internal gradients
Browse files Browse the repository at this point in the history 
Another important step in #11.

Also removed many debug prints.
  • Loading branch information
@egpbos
egpbos committed Nov 21, 2017
1 parent 52c5fcc commit 4fccb6c
Show file tree
Hide file tree
Showing 3 changed files with 58 additions and 41 deletions.
36 changes: 18 additions & 18 deletions math/minuit2/src/MnSeedGenerator.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -73,19 +73,19 @@ MinimumSeed MnSeedGenerator::operator()(const MnFcn& fcn, const GradientCalculat
}

MinimumParameters pa(x, fcnmin);
std::cout << "... doing gc(pa) ..." << std::endl;
// std::cout << "... doing gc(pa) ..." << std::endl;
FunctionGradient dgrad = gc(pa);
std::cout << "dgrad.Vec: " << dgrad.Vec() << std::endl;
std::cout << "dgrad.G2: " << dgrad.G2() << std::endl;
// std::cout << "dgrad.Vec: " << dgrad.Vec() << std::endl;
// std::cout << "dgrad.G2: " << dgrad.G2() << std::endl;
MnAlgebraicSymMatrix mat(n);
double dcovar = 1.;
if(st.HasCovariance()) {
std::cout << "has covariance" << std::endl;
// std::cout << "has covariance" << std::endl;
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 {
std::cout << "has no covariance" << std::endl;
// std::cout << "has no covariance" << std::endl;
for(unsigned int i = 0; i < n; i++)
mat(i,i) = (fabs(dgrad.G2()(i)) > prec.Eps2() ? 1./dgrad.G2()(i) : 1.);
}
Expand All @@ -99,10 +99,10 @@ MinimumSeed MnSeedGenerator::operator()(const MnFcn& fcn, const GradientCalculat
}

NegativeG2LineSearch ng2ls;
std::cout << dgrad.Vec() << std::endl;
std::cout << prec << std::endl;
std::cout << ng2ls.HasNegativeG2(dgrad, prec) << std::endl;
std::cout << dgrad.Vec().size() << std::endl;
// std::cout << dgrad.Vec() << std::endl;
// std::cout << prec << std::endl;
// std::cout << ng2ls.HasNegativeG2(dgrad, prec) << std::endl;
// std::cout << dgrad.Vec().size() << std::endl;

if(ng2ls.HasNegativeG2(dgrad, prec)) {
#ifdef DEBUG
Expand Down Expand Up @@ -160,12 +160,12 @@ MinimumSeed MnSeedGenerator::operator()(const MnFcn& fcn, const AnalyticalGradie
// FunctionGradient tmp = igc(pa);
// std::cout << "tmp.G2: " << tmp.G2() << std::endl;
FunctionGradient grd = gc(pa);
std::cout << "grd.G2: " << grd.G2() << std::endl;
// std::cout << "grd.G2: " << grd.G2() << std::endl;

// FunctionGradient dgrad(grd.Grad(), tmp.G2(), tmp.Gstep());
FunctionGradient dgrad(grd.Grad(), grd.G2(), grd.Gstep());
std::cout << "dgrad.Vec: " << dgrad.Vec() << std::endl;
std::cout << "dgrad.G2: " << dgrad.G2() << std::endl;
// std::cout << "dgrad.Vec: " << dgrad.Vec() << std::endl;
// std::cout << "dgrad.G2: " << dgrad.G2() << std::endl;

if(gc.CheckGradient()) {
bool good = true;
Expand Down Expand Up @@ -195,12 +195,12 @@ MinimumSeed MnSeedGenerator::operator()(const MnFcn& fcn, const AnalyticalGradie
MnAlgebraicSymMatrix mat(n);
double dcovar = 1.;
if(st.HasCovariance()) {
std::cout << "has covariance" << std::endl;
// std::cout << "has covariance" << std::endl;
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 {
std::cout << "has no covariance" << std::endl;
// std::cout << "has no covariance" << std::endl;
for(unsigned int i = 0; i < n; i++)
mat(i,i) = (fabs(dgrad.G2()(i)) > prec.Eps2() ? 1./dgrad.G2()(i) : 1.);
}
Expand All @@ -209,10 +209,10 @@ MinimumSeed MnSeedGenerator::operator()(const MnFcn& fcn, const AnalyticalGradie
MinimumState state(pa, err, dgrad, edm, fcn.NumOfCalls());

NegativeG2LineSearch ng2ls;
std::cout << dgrad.Vec() << std::endl;
std::cout << prec << std::endl;
std::cout << ng2ls.HasNegativeG2(dgrad, prec) << std::endl;
std::cout << dgrad.Vec().size() << std::endl;
// std::cout << dgrad.Vec() << std::endl;
// std::cout << prec << std::endl;
// std::cout << ng2ls.HasNegativeG2(dgrad, prec) << std::endl;
// std::cout << dgrad.Vec().size() << std::endl;

if(ng2ls.HasNegativeG2(dgrad, prec)) {
// Numerical2PGradientCalculator ngc(fcn, st.Trafo(), stra);
Expand Down
2 changes: 2 additions & 0 deletions roofit/roofitcore/inc/NumericalDerivatorMinuit2.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -85,6 +85,8 @@ namespace RooFit {
// mutable std::vector <double> fG2;
// mutable std::vector <double> fGstep;
mutable ROOT::Minuit2::FunctionGradient fG;
mutable ROOT::Minuit2::FunctionGradient fG_internal;

// same story for SetParameterHasLimits
mutable std::vector <bool> _parameter_has_limits;

Expand Down
61 changes: 38 additions & 23 deletions roofit/roofitcore/src/NumericalDerivatorMinuit2.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -45,7 +45,8 @@ namespace RooFit {
Up(1),
fVal(0),
fN(0),
fG(0)//,
fG(0),
fG_internal(0)
// eps(std::numeric_limits<double>::epsilon()),
// eps2(2 * std::sqrt(eps))
{}
Expand All @@ -59,7 +60,8 @@ namespace RooFit {
Up(error_level),
fVal(0),
fN(f.NDim()),
fG(f.NDim())//,
fG(f.NDim()),
fG_internal(f.NDim())
// eps(precision),
// eps2(2 * std::sqrt(eps))
{
Expand Down Expand Up @@ -95,6 +97,7 @@ namespace RooFit {
fVal(other.fVal),
fN(other.fN),
fG(other.fG),
fG_internal(other.fG_internal),
_parameter_has_limits(other._parameter_has_limits),
precision(other.precision)
// eps(other.eps),
Expand All @@ -107,6 +110,7 @@ namespace RooFit {
// fG2 = other.fG2;
// fGstep = other.fGstep;
fG = other.fG;
fG_internal = other.fG_internal;
_parameter_has_limits = other._parameter_has_limits;
fFunction = other.fFunction;
fStepTolerance = other.fStepTolerance;
Expand Down Expand Up @@ -189,6 +193,9 @@ namespace RooFit {
ROOT::Minuit2::MnAlgebraicVector grad_vec(fG.Grad()),
gr2_vec(fG.G2()),
gstep_vec(fG.Gstep());
ROOT::Minuit2::MnAlgebraicVector grad_vec_internal(fG_internal.Grad()),
gr2_vec_internal(fG_internal.G2()),
gstep_vec_internal(fG_internal.Gstep());

// MODIFIED: Up
// In Minuit2, this depends on the type of function to minimize, e.g.
Expand All @@ -214,24 +221,24 @@ namespace RooFit {
std::cout << "fVal = " << fVal << "\t";
std::cout << std::endl;

grad_vec(i) *= DInt2Ext(parameters[i], vx[i]);
gr2_vec(i) *= D2Int2Ext(parameters[i], vx[i]);
gstep_vec(i) *= GStepInt2Ext(parameters[i], vx[i]);
// grad_vec(i) *= DInt2Ext(parameters[i], vx[i]);
// gr2_vec(i) *= D2Int2Ext(parameters[i], vx[i]);
// gstep_vec(i) *= GStepInt2Ext(parameters[i], vx[i]);

std::cout << "BEFORE, INTERNAL: fGrd[" << i <<"] = " << grad_vec(i) << "\t";
std::cout << "fG2[" << i <<"] = " << gr2_vec(i) << "\t";
std::cout << "fGstep[" << i <<"] = " << gstep_vec(i) << "\t";
std::cout << "BEFORE, INTERNAL: fGrd[" << i <<"] = " << grad_vec_internal(i) << "\t";
std::cout << "fG2[" << i <<"] = " << gr2_vec_internal(i) << "\t";
std::cout << "fGstep[" << i <<"] = " << gstep_vec_internal(i) << "\t";
std::cout << "x[" << i << "] = " << vx[i] << "\t";
std::cout << "fVal = " << fVal << "\t";
std::cout << std::endl;

double xtf = vx[i];
double epspri = eps2 + std::abs(grad_vec(i) * eps2);
double epspri = eps2 + std::abs(grad_vec_internal(i) * eps2);
double step_old = 0.;
for (unsigned int j = 0; j < ncycle; ++ j) {

double optstp = std::sqrt(dfmin/(std::abs(gr2_vec(i))+epspri));
double step = std::max(optstp, std::abs(0.1*gstep_vec(i)));
double optstp = std::sqrt(dfmin/(std::abs(gr2_vec_internal(i))+epspri));
double step = std::max(optstp, std::abs(0.1*gstep_vec_internal(i)));

// MODIFIED: in Minuit2 we have here the following condition:
// if(Trafo().Parameter(Trafo().ExtOfInt(i)).HasLimits()) {
Expand All @@ -242,7 +249,7 @@ namespace RooFit {
// See the discussion above NumericalDerivatorMinuit2::SetInitialGradient
// below on how to pass parameter information to this derivator.

double stpmax = 10.*std::abs(gstep_vec(i));
double stpmax = 10.*std::abs(gstep_vec_internal(i));
if (step > stpmax) step = stpmax;

double stpmin = std::max(vrysml, 8.*std::abs(eps2*vx[i])); //8.*std::abs(double(eps2*x[i]))
Expand All @@ -251,7 +258,7 @@ namespace RooFit {
//answer = fGrd[i];
break;
}
gstep_vec(i) = step;
gstep_vec_internal(i) = step;
step_old = step;
// std::cout << "step = " << step << std::endl;
vx[i] = xtf + step;
Expand All @@ -266,37 +273,37 @@ namespace RooFit {
vx[i] = xtf;
vx_external[i] = Int2ext(parameters[i], vx[i]);

double fGrd_old = grad_vec(i);
grad_vec(i) = 0.5*(fs1-fs2)/step;
double fGrd_old = grad_vec_internal(i);
grad_vec_internal(i) = 0.5*(fs1-fs2)/step;
// std::cout << "int i = " << i << std::endl;
// std::cout << "fs1 = " << fs1 << std::endl;
// std::cout << "fs2 = " << fs2 << std::endl;
// std::cout << "fVal = " << fVal << std::endl;
// std::cout << "step^2 = " << (step*step) << std::endl;
// std::cout << std::endl;
gr2_vec(i) = (fs1 + fs2 -2.*fVal)/step/step;
gr2_vec_internal(i) = (fs1 + fs2 -2.*fVal)/step/step;

// MODIFIED:
// The condition below had a closing parenthesis differently than
// Minuit. Fixed in this version.
// if (std::abs(fGrd_old-fGrd[i])/(std::abs(fGrd[i]+dfmin/step)) < grad_tolerance)
if (std::abs(fGrd_old - grad_vec(i))/(std::abs(grad_vec(i)) + dfmin/step) < grad_tolerance)
if (std::abs(fGrd_old - grad_vec_internal(i))/(std::abs(grad_vec_internal(i)) + dfmin/step) < grad_tolerance)
{
//answer = fGrd[i];
break;
}
}

std::cout << "AFTER, INTERNAL: fGrd[" << i <<"] = " << grad_vec(i) << "\t";
std::cout << "fG2[" << i <<"] = " << gr2_vec(i) << "\t";
std::cout << "fGstep[" << i <<"] = " << gstep_vec(i) << "\t";
std::cout << "AFTER, INTERNAL: fGrd[" << i <<"] = " << grad_vec_internal(i) << "\t";
std::cout << "fG2[" << i <<"] = " << gr2_vec_internal(i) << "\t";
std::cout << "fGstep[" << i <<"] = " << gstep_vec_internal(i) << "\t";
std::cout << "x[" << i << "] = " << xtf << "\t";
std::cout << "fVal = " << fVal << "\t";
std::cout << std::endl;

grad_vec(i) /= DInt2Ext(parameters[i], xtf);
gr2_vec(i) /= D2Int2Ext(parameters[i], xtf);
gstep_vec(i) /= GStepInt2Ext(parameters[i], xtf);
grad_vec(i) = grad_vec_internal(i) / DInt2Ext(parameters[i], xtf);
gr2_vec(i) = gr2_vec_internal(i) / D2Int2Ext(parameters[i], xtf);
gstep_vec(i) = gstep_vec_internal(i) / GStepInt2Ext(parameters[i], xtf);

std::cout << "AFTER, EXTERNAL: fGrd[" << i <<"] = " << grad_vec(i) << "\t";
std::cout << "fG2[" << i <<"] = " << gr2_vec(i) << "\t";
Expand All @@ -315,6 +322,7 @@ namespace RooFit {
// }

fG = ROOT::Minuit2::FunctionGradient(grad_vec, gr2_vec, gstep_vec);
fG_internal = ROOT::Minuit2::FunctionGradient(grad_vec_internal, gr2_vec_internal, gstep_vec_internal);
return fG;
}

Expand Down Expand Up @@ -437,6 +445,9 @@ namespace RooFit {
ROOT::Minuit2::MnAlgebraicVector grad_vec(fFunction->NDim()),
gr2_vec(fFunction->NDim()),
gstep_vec(fFunction->NDim());
ROOT::Minuit2::MnAlgebraicVector grad_vec_internal(fG_internal.Grad()),
gr2_vec_internal(fG_internal.G2()),
gstep_vec_internal(fG_internal.Gstep());

unsigned ix = 0;
for (auto parameter = parameters.begin(); parameter != parameters.end(); ++parameter, ++ix) {
Expand Down Expand Up @@ -522,6 +533,9 @@ namespace RooFit {
gr2_vec(ix) = g2 / D2Int2Ext(*parameter, var);
gstep_vec(ix) = gstep / GStepInt2Ext(*parameter, var);

grad_vec_internal(ix) = grd;
gr2_vec_internal(ix) = g2;
gstep_vec_internal(ix) = gstep;

std::cout << "INTERNAL: fGrd[" << ix <<"] = " << grd << "\t";
std::cout << "fG2[" << ix <<"] = " << g2 << "\t";
Expand All @@ -533,6 +547,7 @@ namespace RooFit {
}

fG = ROOT::Minuit2::FunctionGradient(grad_vec, gr2_vec, gstep_vec);
fG_internal = ROOT::Minuit2::FunctionGradient(grad_vec_internal, gr2_vec_internal, gstep_vec_internal);
}


Expand Down

0 comments on commit 4fccb6c

Please sign in to comment.