Skip to content

Commit

Permalink
Got NumericalDerivatorMinuit2 to find exact same minimum
Browse files Browse the repository at this point in the history 
Important step in #11.
  • Loading branch information
@egpbos
egpbos committed Nov 20, 2017
1 parent caac087 commit 52c5fcc
Show file tree
Hide file tree
Showing 3 changed files with 64 additions and 20 deletions.
16 changes: 15 additions & 1 deletion math/minuit2/src/Numerical2PGradientCalculator.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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++) {
Expand Down Expand Up @@ -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";
Expand Down
65 changes: 46 additions & 19 deletions roofit/roofitcore/src/NumericalDerivatorMinuit2.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -167,17 +167,24 @@ namespace RooFit {
std::cout << "########### NumericalDerivatorMinuit2::Differentiate()" <<std::endl;

assert(fFunction != 0);
std::vector<double> vx(fFunction->NDim());
std::vector<double> 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()),
Expand All @@ -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) {
Expand All @@ -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];
Expand All @@ -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] <<std::endl;
double fs1 = f(x);
double fs1 = f(vx_external.data());
//std::cout << "xtf + step = " << x[i] << ", fs1 = " << fs1 << std::endl;
x[i] = xtf - step;
double fs2 = f(x);
vx[i] = xtf - step;
vx_external[i] = Int2ext(parameters[i], vx[i]);
double fs2 = f(vx_external.data());
//std::cout << "xtf - step = " << x[i] << ", fs2 = " << fs2 << std::endl;
x[i] = xtf;

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;
Expand All @@ -260,21 +287,21 @@ namespace RooFit {
}
}

std::cout << "INTERNAL: fGrd[" << i <<"] = " << grad_vec(i) << "\t";
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 << "x[" << i << "] = " << xtf << "\t";
std::cout << "fVal = " << fVal << "\t";
std::cout << std::endl;

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

std::cout << "EXTERNAL: fGrd[" << i <<"] = " << grad_vec(i) << "\t";
std::cout << "AFTER, 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 << "] = " << xtf << "\t";
std::cout << "x[" << i << "] = " << vx_external[i] << "\t";
std::cout << "fVal = " << fVal << "\t";
std::cout << std::endl;

Expand Down
3 changes: 3 additions & 0 deletions roofit/roofitcore/src/RooGradMinimizerFcn.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -649,6 +649,9 @@ void RooGradMinimizerFcn::InitGradient() const {


void RooGradMinimizerFcn::run_derivator(const double *x) const {
for (int i = 0; i < _nDim; ++i) {
std::cout << "x[" << i << "] = " << x[i] << std::endl;
}
if (!_grad_initialized) {
InitGradient();
}
Expand Down

0 comments on commit 52c5fcc

Please sign in to comment.