From 59a417dca514ae3618cd3e773c45a4daad2fafbd Mon Sep 17 00:00:00 2001 From: "E. G. Patrick Bos" Date: Wed, 1 Nov 2017 17:03:40 +0100 Subject: [PATCH] Debugging initial gradient (#11) + adding Int2ext/Ext2int --- .../inc/Math/NumericalDerivatorMinuit2.h | 11 ++++ .../src/NumericalDerivatorMinuit2.cxx | 59 ++++++++++++++++++- math/minuit2/src/DavidonErrorUpdator.cxx | 4 ++ .../minuit2/src/InitialGradientCalculator.cxx | 21 +++++++ math/minuit2/src/VariableMetricBuilder.cxx | 1 + 5 files changed, 94 insertions(+), 2 deletions(-) diff --git a/math/mathcore/inc/Math/NumericalDerivatorMinuit2.h b/math/mathcore/inc/Math/NumericalDerivatorMinuit2.h index d62910404c2fc..0380cb44f2c01 100644 --- a/math/mathcore/inc/Math/NumericalDerivatorMinuit2.h +++ b/math/mathcore/inc/Math/NumericalDerivatorMinuit2.h @@ -25,6 +25,10 @@ #include #include "Fit/ParameterSettings.h" +#include "Minuit2/SinParameterTransformation.h" +#include "Minuit2/SqrtUpParameterTransformation.h" +#include "Minuit2/SqrtLowParameterTransformation.h" + namespace ROOT { namespace Math { @@ -55,6 +59,9 @@ class NumericalDerivatorMinuit2 { void SetInitialValues(const double* g, const double* g2, const double* s); + double Int2ext(const ROOT::Fit::ParameterSettings& parameter, double val) const; + double Ext2int(const ROOT::Fit::ParameterSettings& parameter, double val) const; + void SetInitialGradient(std::vector& parameters) const; void SetParameterHasLimits(std::vector& parameters) const; @@ -83,6 +90,10 @@ class NumericalDerivatorMinuit2 { double eps; double eps2; + SinParameterTransformation fDoubleLimTrafo; + SqrtUpParameterTransformation fUpperLimTrafo; + SqrtLowParameterTransformation fLowerLimTrafo; + }; diff --git a/math/mathcore/src/NumericalDerivatorMinuit2.cxx b/math/mathcore/src/NumericalDerivatorMinuit2.cxx index 24e94db4c9c49..0f876862fb138 100644 --- a/math/mathcore/src/NumericalDerivatorMinuit2.cxx +++ b/math/mathcore/src/NumericalDerivatorMinuit2.cxx @@ -33,8 +33,6 @@ #include "Fit/ParameterSettings.h" #include // needed here because in Fitter is only a forward declaration -//#include "Minuit2/MnStrategy.h" -//#include "Minuit2/MnMachinePrecision.h" namespace ROOT { @@ -277,6 +275,36 @@ void NumericalDerivatorMinuit2::SetParameterHasLimits(std::vectorValue(); + std::cout << "var: " << var << std::endl; // Judging by the ParameterSettings.h constructor argument name "err", // I'm guessing what MINUIT calls "Error" is stepsize on the ROOT side. // double werr = parameter->Error(); double werr = parameter->StepSize(); + std::cout << "werr: " << werr << std::endl; + // double sav = Trafo().Int2ext(i, var); // Int2Ext is not necessary, we're doing everything externally here // double sav2 = sav + werr; double sav2 = var + werr; + std::cout << "sav2: " << sav2 << std::endl; + // if(parameter->HasLimits()) { // this if statement in MINUIT is superfluous if(parameter->HasUpperLimit() && sav2 > parameter->UpperLimit()) { sav2 = parameter->UpperLimit(); } + std::cout << "sav2: " << sav2 << std::endl; + // double var2 = Trafo().Ext2int(exOfIn, sav2); // Ext2int is not necessary, we're doing everything externally here // double vplu = var2 - var; double vplu = sav2 - var; + std::cout << "vplu: " << vplu << std::endl; + // sav2 = sav - werr; sav2 = var - werr; + std::cout << "sav2: " << sav2 << std::endl; + // if(parameter->HasLimits()) { // this if statement in MINUIT is superfluous if(parameter->HasLowerLimit() && sav2 < parameter->LowerLimit()) { sav2 = parameter->LowerLimit(); } + std::cout << "sav2: " << sav2 << std::endl; + // var2 = Trafo().Ext2int(exOfIn, sav2); // Ext2int is not necessary, we're doing everything externally here // double vmin = var2 - var; double vmin = sav2 - var; + std::cout << "vmin: " << vmin << std::endl; + double gsmin = 8. * eps2 * (fabs(var) + eps2); // protect against very small step sizes which can cause dirin to zero and then nan values in grd double dirin = std::max(0.5*(fabs(vplu) + fabs(vmin)), gsmin ); + std::cout << "gsmin: " << gsmin << std::endl; + std::cout << "dirin: " << dirin << std::endl; + // double g2 = 2.0*fFcn.ErrorDef()/(dirin*dirin); // ErrorDef is the same as Up, which we already have in here double g2 = 2.0*Up/(dirin*dirin); double gstep = std::max(gsmin, 0.1*dirin); double grd = g2*dirin; + + std::cout << "gstep: " << gstep << std::endl; + if(parameter->IsBound()) { if(gstep > 0.5) gstep = 0.5; } fGrd[ix] = grd; fG2[ix] = g2; fGstep[ix] = gstep; + + std::cout << "fGrd[" << ix <<"] = " << fGrd[ix] << "\t"; + std::cout << "fG2[" << ix <<"] = " << fG2[ix] << "\t"; + std::cout << "fGstep[" << ix <<"] = " << fGstep[ix] << std::endl; + } } diff --git a/math/minuit2/src/DavidonErrorUpdator.cxx b/math/minuit2/src/DavidonErrorUpdator.cxx index 2d7302e36434d..e0c7b3d850d0c 100644 --- a/math/minuit2/src/DavidonErrorUpdator.cxx +++ b/math/minuit2/src/DavidonErrorUpdator.cxx @@ -41,6 +41,10 @@ MinimumError DavidonErrorUpdator::Update(const MinimumState& s0, MnAlgebraicVector dx = p1.Vec() - s0.Vec(); MnAlgebraicVector dg = g1.Vec() - s0.Gradient().Vec(); +// std::cout << "v0: " << v0 << std::endl; +// std::cout << "s0.Vec(): " << s0.Vec() << std::endl; +// std::cout << "s0.Gradient().Vec(): " << s0.Gradient().Vec() << std::endl; + double delgam = inner_product(dx, dg); double gvg = similarity(dg, v0); diff --git a/math/minuit2/src/InitialGradientCalculator.cxx b/math/minuit2/src/InitialGradientCalculator.cxx index 3292e417f8e75..31fb658102aef 100644 --- a/math/minuit2/src/InitialGradientCalculator.cxx +++ b/math/minuit2/src/InitialGradientCalculator.cxx @@ -43,34 +43,51 @@ FunctionGradient InitialGradientCalculator::operator()(const MinimumParameters& #endif MnAlgebraicVector gr(n), gr2(n), gst(n); + std::cout << "initial gradient minuit: " << std::endl; for(unsigned int i = 0; i < n; i++) { unsigned int exOfIn = Trafo().ExtOfInt(i); double var = par.Vec()(i); + std::cout << "var: " << var << std::endl; + double werr = Trafo().Parameter(exOfIn).Error(); + std::cout << "werr: " << werr << std::endl; double sav = Trafo().Int2ext(i, var); + std::cout << "sav: " << sav << std::endl; + double sav2 = sav + werr; + std::cout << "sav2: " << sav2 << std::endl; if(Trafo().Parameter(exOfIn).HasLimits()) { if(Trafo().Parameter(exOfIn).HasUpperLimit() && sav2 > Trafo().Parameter(exOfIn).UpperLimit()) sav2 = Trafo().Parameter(exOfIn).UpperLimit(); } + std::cout << "sav2: " << sav2 << std::endl; double var2 = Trafo().Ext2int(exOfIn, sav2); + std::cout << "var2: " << var2 << std::endl; double vplu = var2 - var; + std::cout << "vplu: " << vplu << std::endl; sav2 = sav - werr; + std::cout << "sav2: " << sav2 << std::endl; if(Trafo().Parameter(exOfIn).HasLimits()) { if(Trafo().Parameter(exOfIn).HasLowerLimit() && sav2 < Trafo().Parameter(exOfIn).LowerLimit()) sav2 = Trafo().Parameter(exOfIn).LowerLimit(); } + std::cout << "sav2: " << sav2 << std::endl; var2 = Trafo().Ext2int(exOfIn, sav2); + std::cout << "var2: " << var2 << std::endl; double vmin = var2 - var; + std::cout << "vmin: " << vmin << std::endl; double gsmin = 8.*Precision().Eps2()*(fabs(var) + Precision().Eps2()); + std::cout << "gsmin: " << gsmin << std::endl; // protect against very small step sizes which can cause dirin to zero and then nan values in grd double dirin = std::max(0.5*(fabs(vplu) + fabs(vmin)), gsmin ); + std::cout << "dirin: " << dirin << std::endl; double g2 = 2.0*fFcn.ErrorDef()/(dirin*dirin); double gstep = std::max(gsmin, 0.1*dirin); + std::cout << "gstep: " << gstep << std::endl; double grd = g2*dirin; if(Trafo().Parameter(exOfIn).HasLimits()) { if(gstep > 0.5) gstep = 0.5; @@ -79,6 +96,10 @@ FunctionGradient InitialGradientCalculator::operator()(const MinimumParameters& gr2(i) = g2; gst(i) = gstep; + std::cout << "fGrd[" << i <<"] = " << gr(i) << "\t"; + std::cout << "fG2[" << i <<"] = " << gr2(i) << "\t"; + std::cout << "fGstep[" << i <<"] = " << gst(i) << std::endl; + #ifdef DEBUG std::cout << "computing initial gradient for parameter " << Trafo().Name(exOfIn) << " value = " << var << " [ " << vmin << " , " << vplu << " ] " << "dirin " << dirin << " grd " << grd << " g2 " << g2 << std::endl; diff --git a/math/minuit2/src/VariableMetricBuilder.cxx b/math/minuit2/src/VariableMetricBuilder.cxx index c3097dcce91ef..13da29e7250e3 100644 --- a/math/minuit2/src/VariableMetricBuilder.cxx +++ b/math/minuit2/src/VariableMetricBuilder.cxx @@ -380,6 +380,7 @@ FunctionMinimum VariableMetricBuilder::Minimum(const MnFcn& fcn, const GradientC // update the state s0 = MinimumState(p, e, g, edm, fcn.NumOfCalls()); + std::cout << "s0: " << s0 << std::endl; if (StorageLevel() || result.size() <= 1) AddResult(result, s0); else