From 52c5fcc397ecb9e1e5275e5496610640f33ae4fb Mon Sep 17 00:00:00 2001 From: "E. G. Patrick Bos" Date: Mon, 20 Nov 2017 14:57:27 +0100 Subject: [PATCH] Got NumericalDerivatorMinuit2 to find exact same minimum Important step in #11. --- .../src/Numerical2PGradientCalculator.cxx | 16 ++++- .../src/NumericalDerivatorMinuit2.cxx | 65 +++++++++++++------ roofit/roofitcore/src/RooGradMinimizerFcn.cxx | 3 + 3 files changed, 64 insertions(+), 20 deletions(-) diff --git a/math/minuit2/src/Numerical2PGradientCalculator.cxx b/math/minuit2/src/Numerical2PGradientCalculator.cxx index f97332d803570..c8af46e597862 100644 --- a/math/minuit2/src/Numerical2PGradientCalculator.cxx +++ b/math/minuit2/src/Numerical2PGradientCalculator.cxx @@ -106,6 +106,11 @@ FunctionGradient Numerical2PGradientCalculator::operator()(const MinimumParamete MnAlgebraicVector g2 = Gradient.G2(); MnAlgebraicVector gstep = Gradient.Gstep(); + for (int i = 0; i < n; ++i) { + std::cout << par.Vec()(i) << "\t"; + } + std::cout << std::endl; + #ifndef _OPENMP MPIProcess mpiproc(n,0); #endif @@ -148,7 +153,15 @@ FunctionGradient Numerical2PGradientCalculator::operator()(const MinimumParamete MnAlgebraicVector x = par.Vec(); #endif - double xtf = x(i); + std::cout << "BEFORE: "; + std::cout << "fGrd[" << i <<"] = " << grd(i) << "\t"; + std::cout << "fG2[" << i <<"] = " << g2(i) << "\t"; + std::cout << "fGstep[" << i <<"] = " << gstep(i) << "\t"; + std::cout << "x[" << i << "] = " << x(i) << "\t"; + std::cout << "fVal = " << fcnmin << "\t"; + std::cout << std::endl; + + double xtf = x(i); double epspri = eps2 + fabs(grd(i)*eps2); double stepb4 = 0.; for(unsigned int j = 0; j < ncycle; j++) { @@ -194,6 +207,7 @@ FunctionGradient Numerical2PGradientCalculator::operator()(const MinimumParamete << " grd " << grd(i) << " g2 " << g2(i) << std::endl; std::cout.precision(pr); #endif + std::cout << "AFTER: "; std::cout << "fGrd[" << i <<"] = " << grd(i) << "\t"; std::cout << "fG2[" << i <<"] = " << g2(i) << "\t"; std::cout << "fGstep[" << i <<"] = " << gstep(i) << "\t"; diff --git a/roofit/roofitcore/src/NumericalDerivatorMinuit2.cxx b/roofit/roofitcore/src/NumericalDerivatorMinuit2.cxx index 95f11ded1896b..bd11abf3a3bc4 100644 --- a/roofit/roofitcore/src/NumericalDerivatorMinuit2.cxx +++ b/roofit/roofitcore/src/NumericalDerivatorMinuit2.cxx @@ -167,17 +167,24 @@ namespace RooFit { std::cout << "########### NumericalDerivatorMinuit2::Differentiate()" < vx(fFunction->NDim()); + std::vector vx(fFunction->NDim()), vx_external(fFunction->NDim()); assert (vx.size() > 0); - TODO: x is anders in Numerical2PGradient (zie huidige test output), fix dat! +// double *x = &vx[0]; + std::copy (cx, cx+fFunction->NDim(), vx.data()); + std::copy (cx, cx+fFunction->NDim(), vx_external.data()); + + // convert to Minuit internal parameters + for (int i = 0; i < int(fN); i++) { + std::cout << vx[i] << "\t"; + vx[i] = Ext2int(parameters[i], vx[i]); + std::cout << vx[i] << std::endl; + } - double *x = &vx[0]; - std::copy (cx, cx+fFunction->NDim(), x); double step_tolerance = fStepTolerance; double grad_tolerance = fGradTolerance; const ROOT::Math::IBaseFunctionMultiDim &f = *fFunction; - fVal = f(x); //value of function at given points + fVal = f(vx_external.data()); //value of function at given points ROOT::Minuit2::MnAlgebraicVector grad_vec(fG.Grad()), gr2_vec(fG.G2()), @@ -200,7 +207,25 @@ namespace RooFit { for (int i = 0; i < int(fN); i++) { - double xtf = x[i]; + std::cout << "BEFORE, EXTERNAL: fGrd[" << i <<"] = " << grad_vec(i) << "\t"; + std::cout << "fG2[" << i <<"] = " << gr2_vec(i) << "\t"; + std::cout << "fGstep[" << i <<"] = " << gstep_vec(i) << "\t"; + std::cout << "x[" << i << "] = " << vx_external[i] << "\t"; + 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]); + + 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 << "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 step_old = 0.; for (unsigned int j = 0; j < ncycle; ++ j) { @@ -220,7 +245,7 @@ namespace RooFit { double stpmax = 10.*std::abs(gstep_vec(i)); if (step > stpmax) step = stpmax; - double stpmin = std::max(vrysml, 8.*std::abs(eps2*x[i])); //8.*std::abs(double(eps2*x[i])) + double stpmin = std::max(vrysml, 8.*std::abs(eps2*vx[i])); //8.*std::abs(double(eps2*x[i])) if (step < stpmin) step = stpmin; if (std::abs((step-step_old)/step) < step_tolerance) { //answer = fGrd[i]; @@ -229,15 +254,17 @@ namespace RooFit { gstep_vec(i) = step; step_old = step; // std::cout << "step = " << step << std::endl; - x[i] = xtf + step; + vx[i] = xtf + step; + vx_external[i] = Int2ext(parameters[i], vx[i]); //std::cout << "x[" << i << "] = " << x[i] <