From 4fccb6c7876efca83b8d285ea5c5dee4c031f9ca Mon Sep 17 00:00:00 2001 From: "E. G. Patrick Bos" Date: Tue, 21 Nov 2017 08:16:42 +0100 Subject: [PATCH] NumericalDerivatorMinuit2 now also uses exact same internal gradients Another important step in #11. Also removed many debug prints. --- math/minuit2/src/MnSeedGenerator.cxx | 36 +++++------ .../inc/NumericalDerivatorMinuit2.h | 2 + .../src/NumericalDerivatorMinuit2.cxx | 61 ++++++++++++------- 3 files changed, 58 insertions(+), 41 deletions(-) diff --git a/math/minuit2/src/MnSeedGenerator.cxx b/math/minuit2/src/MnSeedGenerator.cxx index ac65bff7d957e..168a72bb6fbfd 100644 --- a/math/minuit2/src/MnSeedGenerator.cxx +++ b/math/minuit2/src/MnSeedGenerator.cxx @@ -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.); } @@ -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 @@ -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; @@ -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.); } @@ -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); diff --git a/roofit/roofitcore/inc/NumericalDerivatorMinuit2.h b/roofit/roofitcore/inc/NumericalDerivatorMinuit2.h index aac7eb87cc042..a2067a2b60fa1 100644 --- a/roofit/roofitcore/inc/NumericalDerivatorMinuit2.h +++ b/roofit/roofitcore/inc/NumericalDerivatorMinuit2.h @@ -85,6 +85,8 @@ namespace RooFit { // mutable std::vector fG2; // mutable std::vector fGstep; mutable ROOT::Minuit2::FunctionGradient fG; + mutable ROOT::Minuit2::FunctionGradient fG_internal; + // same story for SetParameterHasLimits mutable std::vector _parameter_has_limits; diff --git a/roofit/roofitcore/src/NumericalDerivatorMinuit2.cxx b/roofit/roofitcore/src/NumericalDerivatorMinuit2.cxx index bd11abf3a3bc4..6bda8cbbab228 100644 --- a/roofit/roofitcore/src/NumericalDerivatorMinuit2.cxx +++ b/roofit/roofitcore/src/NumericalDerivatorMinuit2.cxx @@ -45,7 +45,8 @@ namespace RooFit { Up(1), fVal(0), fN(0), - fG(0)//, + fG(0), + fG_internal(0) // eps(std::numeric_limits::epsilon()), // eps2(2 * std::sqrt(eps)) {} @@ -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)) { @@ -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), @@ -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; @@ -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. @@ -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()) { @@ -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])) @@ -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; @@ -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"; @@ -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; } @@ -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) { @@ -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"; @@ -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); }