Debugging initial gradient (#11) + adding Int2ext/Ext2int

math/mathcore/inc/Math/NumericalDerivatorMinuit2.h

@@ -25,6 +25,10 @@
 
 #include <vector>
 #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<ROOT::Fit::ParameterSettings>& parameters) const;
   void SetParameterHasLimits(std::vector<ROOT::Fit::ParameterSettings>& parameters) const;
 
@@ -83,6 +90,10 @@ class NumericalDerivatorMinuit2 {
   double eps;
   double eps2;
 
+  SinParameterTransformation fDoubleLimTrafo;
+  SqrtUpParameterTransformation fUpperLimTrafo;
+  SqrtLowParameterTransformation fLowerLimTrafo;
+
 };
 
 

math/mathcore/src/NumericalDerivatorMinuit2.cxx

@@ -33,8 +33,6 @@
 #include "Fit/ParameterSettings.h"
 
 #include <Math/Minimizer.h>  // 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::vector<ROOT::Fit::Par
   }
 }
 
+double NumericalDerivatorMinuit2::Int2ext(const ROOT::Fit::ParameterSettings& parameter, double val) const {
+  // return external value from internal value for parameter i
+  if(parameter.IsBound()) {
+    if(parameter.IsDoubleBound()) {
+      return fDoubleLimTrafo.Int2ext(val, parameter.UpperLimit(), parameter.LowerLimit());
+    } else if (parameter.HasUpperLimit() && !parameter.HasLowerLimit()) {
+      return fUpperLimTrafo.Int2ext(val, parameter.UpperLimit());
+    } else {
+      return fLowerLimTrafo.Int2ext(val, parameter.LowerLimit());
+    }
+  }
+
+  return val;
+}
+
+double NumericalDerivatorMinuit2::Ext2int(const ROOT::Fit::ParameterSettings& parameter, double val) const {
+  // return the internal value for parameter i with external value val
+
+  if(parameter.IsBound()) {
+    if(parameter.IsDoubleBound()) {
+      return fDoubleLimTrafo.Ext2int(val, parameter.UpperLimit(), parameter.LowerLimit(), Precision());
+    } else if (parameter.HasUpperLimit() && !parameter.HasLowerLimit()) {
+      return fUpperLimTrafo.Ext2int(val, parameter.UpperLimit(), Precision());
+    } else {
+      return fLowerLimTrafo.Ext2int(val, parameter.LowerLimit(), Precision());
+    }
+  }
+
+  return val;
+}
 
 // MODIFIED:
 // This function was not implemented as in Minuit2. Now it copies the behavior
@@ -290,6 +318,7 @@ void NumericalDerivatorMinuit2::SetInitialGradient(std::vector<ROOT::Fit::Parame
 //    Double_t oldVerr = parameters[index].StepSize();
 //    Double_t oldVlo = parameters[index].LowerLimit();
 //    Double_t oldVhi = parameters[index].UpperLimit();
+      std::cout << "initial gradient numericalderivatorminuit2: " << std::endl;
 
   unsigned ix = 0;
   for (auto parameter = parameters.begin(); parameter != parameters.end(); ++parameter, ++ix) {
@@ -314,56 +343,82 @@ void NumericalDerivatorMinuit2::SetInitialGradient(std::vector<ROOT::Fit::Parame
 //    double var = par.Vec()(i);
     // I'm guessing par.Vec()(i) should give the value of variable i...
     double var = parameter->Value();
+    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;
+
   }
 }
 

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);
 

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;

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