Minou.h

Go to the documentation of this file.

#ifndef _utl_Minou_h_
#define _utl_Minou_h_

#include <utl/Math.h>
#include <TMinuit.h>
#include <TMatrixD.h>
//#include <TGraph.h>
#include <vector>
#include <exception>
#include <iostream>


namespace utl {

  namespace Minou {

    struct ParameterDef {

      ParameterDef() = default;

      ParameterDef(const std::string name,
                   const double value = 0, const double step = 1,
                   const double low = 0, const double high = 0,
                   const bool fixed = false) :
        fName(name),
        fValue(value),
        fStep(step),
        fLow(low),
        fHigh(high),
        fFixed(fixed)
      { }

      bool IsFree() const { return !fFixed && fStep; }
      bool IsFixed() const { return !IsFree(); }

      std::string fName;
      double fValue = 0;
      double fStep = 1;
      double fLow = 0;
      double fHigh = 0;
      bool fFixed = false;

    };


    struct MinosParameterResult {
      int fParNumber = 0;
      double fValue = 0;
      double fParabolicError = 0;
      double fLowerError = 0;
      double fUpperError = 0;
      double fGlobalCorrelation = 0;
    };


    /*
      Your functor should inherit from Base.

      The minimizer will take number of parameters, initial values, step size
      and free/fixed information from the Base::GetParameterDefs() so you should
      either initialize your class and its Base with an appropriate vector or set it
      (after creation through your ctor) with Base::SetParameterDefs(vector).

      Minimizer will call operator()(const vector<double>& p) in your
      Base-derived class where "p" are the current values of parameters, which are
      also available (to other methods) with Base::GetCurrentParameters().
    */

    class Base {
    public:
      Base(const unsigned int nPar = 0, const double up = 1)
        : fParameterDefs(nPar), fCurrentParameters(nPar), fErrorDefinition(up) { }

      Base(const std::vector<ParameterDef>& pars, const double up = 1)
        : fParameterDefs(pars), fCurrentParameters(pars.size()), fErrorDefinition(up) { }

      unsigned int GetNParameters() const { return fParameterDefs.size(); }

      unsigned int
      GetNFreeParameters()
        const
      {
        unsigned int n = 0;
        for (const auto& p : fParameterDefs)
          if (p.IsFree())
            ++n;
        return n;
      }

      // parameter definitions including initial values
      std::vector<ParameterDef>& GetParameterDefs() { return fParameterDefs; }
      const std::vector<ParameterDef>& GetParameterDefs() const { return fParameterDefs; }
      void SetParameterDefs(const std::vector<ParameterDef>& defs) { fParameterDefs = defs; }

      std::vector<double>
      GetParameterDefValues()
        const
      {
        std::vector<double> r;
        r.reserve(GetNParameters());
        for (const auto& p : fParameterDefs)
          r.push_back(p.fValue);
        return r;
      }

      void
      SetParameterDefValues(const std::vector<double>& vals)
      {
        const unsigned int n = GetNParameters();
        if (vals.size() != n)
          throw std::runtime_error("SetParameterDefValues: requires the same length");
        for (unsigned int i = 0; i < n; ++i)
          fParameterDefs[i].fValue = vals[i];
      }

      std::vector<int>
      GetParameterDefFixed()
        const
      {
        std::vector<int> r;
        r.reserve(GetNParameters());
        for (const auto& p : fParameterDefs)
          r.push_back(p.fFixed);
        return r;
      }

      void
      SetParameterDefFixed(const std::vector<int>& fixed)
      {
        const unsigned int n = GetNParameters();
        if (fixed.size() != n)
          throw std::runtime_error("SetParameterDefFixed: requires the same length");
        for (unsigned int i = 0; i < n; ++i)
          fParameterDefs[i].fFixed = fixed[i];
      }

      // current values of parameters during the minimization
      std::vector<double>& GetCurrentParameters() { return fCurrentParameters; }
      const std::vector<double>& GetCurrentParameters() const { return fCurrentParameters; }

      void
      SetCurrentParameters(const std::vector<double>& p)
      {
        if (p.size() != GetNParameters())
          throw std::runtime_error("SetCurrentParameters: requires the same number");
        fCurrentParameters = p;
      }

      void SetErrorDefinition(const double up) { fErrorDefinition = up; }

      double GetErrorDefinition() const { return fErrorDefinition; }

      std::string GetName(const int i) const
      { return fParameterDefs.at(i).fName; }

      std::vector<std::string>
      GetNames()
        const
      {
        const int n = GetNParameters();
        std::vector<std::string> r;
        r.reserve(n);
        for (int i = 0; i < n; ++i)
          r.push_back(GetName(i));
        return r;
      }

      unsigned int
      GetIndex(const std::string& parameterName)
        const
      {
        for (unsigned int i = 0, n = fParameterDefs.size(); i < n; ++i) {
          if (parameterName == fParameterDefs[i].fName)
            return i;
        }
        throw std::out_of_range("GetIndex: parameter " + parameterName + " not found!");
      }

    private:
      std::vector<ParameterDef> fParameterDefs;
      std::vector<double> fCurrentParameters;
      double fErrorDefinition = 0;
    };


    /*
      F = functor derived from Base and implements the operator()(const vector<double>& par)
      to reduce copying
    */
    template<class F>
    class Minimizer {
    public:
      Minimizer(F& func) :
        fFunction(&func),
        fMinuit(func.GetNParameters())
      {
        if (fgMinimizer)
          throw std::runtime_error("only one instance at the time allowed");
        fgMinimizer = this;
        SetVerbose(false);
        fMinuit.SetFCN(Function);
        SetErrorDefinition(func.GetErrorDefinition());
        DefineParameters();
      }

      ~Minimizer() { fgMinimizer = nullptr; }

      void
      DefineParameters()
      {
        const auto& ps = fFunction->GetParameterDefs();
        for (unsigned int i = 0, n = ps.size(); i < n; ++i) {
          const auto& p = ps[i];
          if (p.fStep <= 0 && !p.fFixed) {
            const std::string err = "In Minou::Minimizer::DefineParameters(): "
                                    "Parameter " + p.fName + ": "
                                    "step should be > 0! (or fix the parameter)";
            std::cerr << err << std::endl;
            throw std::runtime_error(err);
          }
          fMinuit.DefineParameter(i, p.fName.c_str(), p.fValue, p.fStep, p.fLow, p.fHigh);
          if (p.fFixed)
            fMinuit.FixParameter(i);
        }
      }

      void
      SetVerbose(const bool verbose)
      {
        fArgumentList[0] = (verbose ? 1 : -1);
        fMinuit.mnexcm("SET PRINTOUT", fArgumentList, 1, fErrorFlag);
        if (verbose)
          fMinuit.mnexcm("SET WARNINGS", fArgumentList, 0, fErrorFlag);
        else
          fMinuit.mnexcm("SET NOWARNINGS", fArgumentList, 0, fErrorFlag);
        fMinuit.SetPrintLevel(fArgumentList[0]);
      }

      F& GetFunction() { return *fFunction; }

      TMinuit& GetMinuit() { return fMinuit; }

      int
      Minimize(const int n = 500)
      {
        fArgumentList[0] = n;
        fMinuit.mnexcm("MINIMIZE", fArgumentList, 1, fErrorFlag);
        return fErrorFlag;
      }

      int
      Minos(const int n = 500)
      {
        fArgumentList[0] = n;
        fMinuit.mnexcm("MINOS", fArgumentList, 1, fErrorFlag);
        return fErrorFlag;
      }

      int
      Migrad(const int n = 500)
      {
        fArgumentList[0] = n;
        fErrorFlag = 0;
        fMinuit.mnexcm("MIGRAD", fArgumentList, 1, fErrorFlag);
        return fErrorFlag;
      }

      MinosParameterResult
      GetMinosError(const int parIndex)
        const
      {
        MinosParameterResult eRes;
        eRes.fParNumber = parIndex;
        fMinuit.GetParameter(parIndex, eRes.fValue, eRes.fParabolicError);
        fMinuit.mnerrs(parIndex, eRes.fUpperError, eRes.fLowerError,
                       eRes.fParabolicError, eRes.fGlobalCorrelation);
        return eRes;
      }

      std::vector<MinosParameterResult>
      GetMinosErrors()
        const
      {
        const int n = fFunction->GetNParameters();
        std::vector<MinosParameterResult> r;
        r.reserve(n);
        for (int i = 0; i < n; ++i) {
          const auto mRes = GetMinosError(i);
          r.push_back(mRes);
        }
        return r;
      }

      std::pair<double, double>
      GetParameterAndError(const int i)
        const
      {
        std::pair<double, double> ve;
        fMinuit.GetParameter(i, ve.first, ve.second);
        return ve;
      }

      std::vector<std::pair<double, double>>
      GetParametersAndErrors()
        const
      {
        const int n = fFunction->GetNParameters();
        std::vector<std::pair<double, double>> r;
        r.reserve(n);
        for (int i = 0; i < n; ++i) {
          const auto p = GetParameterAndError(i);
          r.push_back(p);
        }
        return r;
      }

      double GetParameter(const int i) const
      { return GetParameterAndError(i).first; }

      std::vector<double>
      GetParameters()
        const
      {
        const int n = fFunction->GetNParameters();
        std::vector<double> r;
        r.reserve(n);
        for (int i = 0; i < n; ++i)
          r.push_back(GetParameter(i));
        return r;
      }

      void WriteBack() const { WriteBack(*fFunction); }

      void WriteBack(Base& base) const
      { base.GetCurrentParameters() = GetParameters(); }

      double GetMinimum() const { return fMinuit.fAmin; }

      double EvaluateFunction(const std::vector<double>& p) const
      { return (*fFunction)(p); }

      TMatrixD
      GetCovarianceMatrix()
        const
      {
        const int n = fMinuit.GetNumFreePars();
        double matrix[n][n];
        for (int i = 0; i < n; i++)
          for (int j = 0; j < n; j++)
            matrix[i][j] = 0;
        fMinuit.mnemat(&matrix[0][0], n);

        const int nTotal = fFunction->GetNParameters();
        int inverseIndex[nTotal];
        const auto& ps = fFunction->GetParameterDefs();

        {
          int iFree = -1;
          for (int i = 0; i < nTotal; ++i)
            if (ps[i].IsFree())
              inverseIndex[i] = ++iFree;
            else
              inverseIndex[i] = 0;
        }

        TMatrixD fullMatrix(nTotal, nTotal);

        for (int i = 0; i < nTotal; ++i) {
          const bool iIsFree = ps[i].IsFree();
          const auto& mii = matrix[inverseIndex[i]];
          auto fmi = fullMatrix[i];
          for (int j = i; j < nTotal; ++j) {
            const double covariance = (iIsFree && ps[j].IsFree()) ? mii[inverseIndex[j]] : 0;
            const auto el = fmi[j] = covariance;
            if (i != j)
              fullMatrix[j][i] = el;
          }
        }

        return fullMatrix;
      }

      /*TGraph
      ParameterScan(const int parNumber, const double minValue,
                    const double maxValue, const int nPoints = 100)
      {
        fArgumentList[0] = parNumber;
        fArgumentList[1] = nPoints;
        fArgumentList[2] = minValue;
        fArgumentList[3] = maxValue;
        fMinuit.mnexcm("SCAN", fArgumentList, 4, fErrorFlag);
        return *(TGraph*)(fMinuit.GetPlot());
      }*/

    private:
      void
      SetErrorDefinition(const double up)
      {
        fArgumentList[0] = up;
        fMinuit.mnexcm("SET ERRORDEF", fArgumentList, 1, fErrorFlag);
      }

      static
      void
      Function(int& /*npar*/, double* /*gin*/, double& f, double* par, int /*iflag*/)
      {
        auto& func = *fgMinimizer->fFunction;
        auto& fpars = func.GetCurrentParameters();
        fpars.assign(par, par + func.GetNParameters());
        f = func(fpars);
      }

      static Minimizer* fgMinimizer;
      F* fFunction = nullptr;
      mutable TMinuit fMinuit;
      int fErrorFlag = 0;
      double fArgumentList[10] = { 0 };
    };

    template<class Functor> Minimizer<Functor>* Minimizer<Functor>::fgMinimizer = nullptr;

  }

}


#endif