Merge pull request #7 from fermi-lat/optimizers-1

Optimizers 1
fermi-lat · Jan 30, 2020 · 2d58f58 · 2d58f58
2 parents 6e5ff5d + a46ece5
commit 2d58f58
Show file tree

Hide file tree

Showing 2 changed files with 100 additions and 0 deletions.
diff --git a/optimizers/NewMinuit.h b/optimizers/NewMinuit.h
@@ -15,6 +15,7 @@
 #include "Minuit2/MnStrategy.h"
 #include "Minuit2/MnUserParameterState.h"
 #include "Minuit2/FunctionMinimum.h"
+#include "Minuit2/MinimumError.h"
 
 namespace optimizers {
 
@@ -72,9 +73,43 @@ namespace optimizers {
         return m_strategy_value;
      }
 
+
+    const ROOT::Minuit2::MnUserParameterState& userState() const;
+    const ROOT::Minuit2::MinimumError& minuitError() const;
+
+    std::vector<std::vector<double> >  userCovariance() const;
+    const std::vector<double>& userGlobalCC() const;
+
+    // hessian (inverse of covariance matrix)
+    std::vector<std::vector<double> > userHessian() const;
+
+    // covariance matrix status (0 = not valid, 1 approximate, 2, full but made pos def, 3 accurate and not pos def
+    int userCovarianceStatus() const { return userState().CovarianceStatus(); } 
+
+    bool userIsValid() const {return userState().IsValid();}
+    bool userHasCovariance() const {return userState().HasCovariance();}
+    bool userHasGlobalCC() const {return userState().HasGlobalCC();}    
+    double userFval() const {return userState().Fval();}
+    double userEdm() const {return userState().Edm();}
+    unsigned int userNFcn() const {return userState().NFcn();}  
+
+
+    double minuitDcovar() const {return minuitError().Dcovar();}
+    bool minuitIsAccurate() const {return minuitError().IsAccurate();}
+    bool minuitIsValid() const {return minuitError().IsValid();}
+    bool minuitIsPosDef() const {return minuitError().IsPosDef();}
+    bool minuitIsMadePosDef() const {return minuitError().IsMadePosDef();}
+    bool minuitHesseFailed() const {return minuitError().HesseFailed();}
+    bool minuitInvertFailed() const {return minuitError().InvertFailed();}
+    bool minuitIsAvailable() const {return minuitError().IsAvailable();}
+
+    std::vector<std::vector<double> > minuitInvHessian() const;
+
     double getDistance(void) const {return m_distance;};
     virtual const std::vector<double> & getUncertainty(bool useBase = false);
     virtual std::vector<std::vector<double> > covarianceMatrix() const;
+
+
     virtual std::ostream& put (std::ostream& s) const;
     std::pair<double,double> Minos(unsigned int n, double level=1., bool numericDeriv=false);
 

diff --git a/src/NewMinuit.cxx b/src/NewMinuit.cxx
@@ -18,6 +18,7 @@
 #include "Minuit2/FunctionMinimum.h"
 #include "Minuit2/MnHesse.h"
 #include "Minuit2/MnPrint.h"
+#include "Minuit2/MnMatrix.h"
 #include "optimizers/Exception.h"
 #include "optimizers/OutOfBounds.h"
 #include "StMnMinos.h"
@@ -300,6 +301,18 @@ namespace optimizers {
     return grad;
   }
 
+  const ROOT::Minuit2::MnUserParameterState& NewMinuit::userState() const {
+    static const ROOT::Minuit2::MnUserParameterState null_state;
+    return m_min == 0 ? null_state : m_min->UserState();
+  }
+
+
+  const ROOT::Minuit2::MinimumError& NewMinuit::minuitError() const {
+    static const ROOT::Minuit2::MinimumError null_error(1);
+    return m_min == 0 ? null_error : m_min->Error();
+  }
+
+
   // Get the uncertainty values from covariance matrix
    const std::vector<double> & NewMinuit::getUncertainty(bool useBase) {
       std::vector<double> parValues;
@@ -324,6 +337,58 @@ namespace optimizers {
     return s;
   }
 
+  std::vector<std::vector<double> >  NewMinuit::userCovariance() const {
+    std::vector<std::vector<double> > cov;
+    if ( m_min == 0 ) return cov;
+    const ROOT::Minuit2::MnUserCovariance& userCov = m_min->UserState().Covariance();
+
+    for (unsigned int x = 0; x < userCov.Nrow(); ++x) {
+      std::vector<double> vec;
+      for (unsigned int y = 0; y < userCov.Nrow(); ++y) {
+	vec.push_back(userCov(x,y));
+      }
+      cov.push_back(vec);
+    }
+    return cov;
+  }
+
+  const std::vector<double>& NewMinuit::userGlobalCC() const {
+    static const std::vector<double> null_globalcc;
+    return m_min == 0 ? null_globalcc : m_min->UserState().GlobalCC().GlobalCC();    
+  }
+
+
+  std::vector<std::vector<double> >  NewMinuit::userHessian() const {
+    std::vector<std::vector<double> > hesse;
+    if ( m_min == 0 ) return hesse;
+    ROOT::Minuit2::MnUserCovariance mn_hesse(m_min->UserState().Hessian());
+
+    for (unsigned int x = 0; x < mn_hesse.Nrow(); ++x) {
+      std::vector<double> vec;
+      for (unsigned int y = 0; y < mn_hesse.Nrow(); ++y) {
+	vec.push_back(mn_hesse(x,y));
+      }
+      hesse.push_back(vec);
+    }
+    return hesse;
+  }
+
+
+  std::vector<std::vector<double> > NewMinuit::minuitInvHessian() const {
+    std::vector<std::vector<double> > inv_hesse;
+    if ( m_min == 0 ) return inv_hesse;
+    const ROOT::Minuit2::MnAlgebraicSymMatrix& inv_mat = m_min->Error().InvHessian();
+
+    for (unsigned int x = 0; x < inv_mat.Nrow(); ++x) {
+      std::vector<double> vec;
+      for (unsigned int y = 0; y < inv_mat.Nrow(); ++y) {
+	vec.push_back(inv_mat(x,y));
+      }
+      inv_hesse.push_back(vec);
+    }
+    return inv_hesse;
+  }
+
    std::vector<std::vector<double> > NewMinuit::covarianceMatrix() const {
       std::vector<double> parValues;
       m_stat->getFreeParamValues(parValues);