MnUserParameterState.cxx

Go to the documentation of this file.

// @(#)root/minuit2:$Id$
// Authors: M. Winkler, F. James, L. Moneta, A. Zsenei   2003-2005
 
/**********************************************************************
 *                                                                    *
 * Copyright (c) 2005 LCG ROOT Math team,  CERN/PH-SFT                *
 *                                                                    *
 **********************************************************************/
 
#include "Minuit2/MnUserParameterState.h"
#include "Minuit2/MnCovarianceSqueeze.h"
#include "Minuit2/MinimumState.h"
#include "Minuit2/MnPrint.h"
 
namespace ROOT {
 
namespace Minuit2 {
 
//
// construct from user parameters (before minimization)
//
MnUserParameterState::MnUserParameterState(std::span<const double> par, std::span<const double> err)
   : fValid(true), fCovStatus(-1), fParameters(MnUserParameters(par, err)), fIntParameters(par.begin(), par.end())
{
}
 
MnUserParameterState::MnUserParameterState(const MnUserParameters &par) : fValid(true), fCovStatus(-1), fParameters(par)
{
   // construct from user parameters (before minimization)
 
   for (auto const &ipar : MinuitParameters()) {
      if (ipar.IsConst() || ipar.IsFixed())
         continue;
      if (ipar.HasLimits())
         fIntParameters.push_back(Ext2int(ipar.Number(), ipar.Value()));
      else
         fIntParameters.push_back(ipar.Value());
   }
}
 
//
// construct from user parameters + errors (before minimization)
//
MnUserParameterState::MnUserParameterState(std::span<const double> par, std::span<const double> cov, unsigned int nrow)
   : fValid(true), fCovStatus(-1), fIntParameters(par.begin(), par.end())
{
   // construct from user parameters + errors (before minimization) using std::vector for parameter error and    // an
   // std::vector of size n*(n+1)/2 for the covariance matrix  and n (rank of cov matrix)
 
   std::vector<double> err;
   err.reserve(par.size());
   for (unsigned int i = 0; i < par.size(); i++) {
      assert(fCovariance(i, i) > 0.);
      err.push_back(std::sqrt(fCovariance(i, i)));
   }
   fParameters = MnUserParameters(par, err);
   MnUserParameterState::AddCovariance({cov, nrow});
}
 
MnUserParameterState::MnUserParameterState(std::span<const double> par, const MnUserCovariance &cov)
   : fValid(true), fCovStatus(-1), fIntParameters(par.begin(), par.end())
{
   // construct from user parameters + errors (before minimization) using std::vector (params) and MnUserCovariance
   // class
 
   std::vector<double> err;
   err.reserve(par.size());
   for (unsigned int i = 0; i < par.size(); i++) {
      assert(fCovariance(i, i) > 0.);
      err.push_back(std::sqrt(fCovariance(i, i)));
   }
   fParameters = MnUserParameters(par, err);
   MnUserParameterState::AddCovariance(cov);
}
 
MnUserParameterState::MnUserParameterState(const MnUserParameters &par, const MnUserCovariance &cov)
   : fValid(true), fCovStatus(-1), fParameters(par)
{
   // construct from user parameters + errors (before minimization) using
   // MnUserParameters and MnUserCovariance objects
 
   for (auto const &ipar : MinuitParameters()) {
      if (ipar.IsConst() || ipar.IsFixed())
         continue;
      if (ipar.HasLimits())
         fIntParameters.push_back(Ext2int(ipar.Number(), ipar.Value()));
      else
         fIntParameters.push_back(ipar.Value());
   }
 
   MnUserParameterState::AddCovariance(cov);
}
 
//
//
MnUserParameterState::MnUserParameterState(const MinimumState &st, double up, const MnUserTransformation &trafo)
   : fValid(st.IsValid()), fCovStatus(-1), fFVal(st.Fval()), fEDM(st.Edm()), fNFcn(st.NFcn())
{
   //
   // construct from internal parameters (after minimization)
   //
   // std::cout << "build a MnUserParameterState after minimization.." << std::endl;
 
   for (auto const &ipar : trafo.Parameters()) {
      if (ipar.IsConst()) {
         Add(ipar.GetName(), ipar.Value());
      } else if (ipar.IsFixed()) {
         Add(ipar.GetName(), ipar.Value(), ipar.Error());
         if (ipar.HasLimits()) {
            if (ipar.HasLowerLimit() && ipar.HasUpperLimit())
               SetLimits(ipar.GetName(), ipar.LowerLimit(), ipar.UpperLimit());
            else if (ipar.HasLowerLimit() && !ipar.HasUpperLimit())
               SetLowerLimit(ipar.GetName(), ipar.LowerLimit());
            else
               SetUpperLimit(ipar.GetName(), ipar.UpperLimit());
         }
         Fix(ipar.GetName());
      } else if (ipar.HasLimits()) {
         unsigned int i = trafo.IntOfExt(ipar.Number());
         double err =
            st.Error().IsValid() ? std::sqrt(2. * up * st.Error().InvHessian()(i, i)) : st.Parameters().Dirin()(i);
         Add(ipar.GetName(), trafo.Int2ext(i, st.Vec()(i)), trafo.Int2extError(i, st.Vec()(i), err));
         if (ipar.HasLowerLimit() && ipar.HasUpperLimit())
            SetLimits(ipar.GetName(), ipar.LowerLimit(), ipar.UpperLimit());
         else if (ipar.HasLowerLimit() && !ipar.HasUpperLimit())
            SetLowerLimit(ipar.GetName(), ipar.LowerLimit());
         else
            SetUpperLimit(ipar.GetName(), ipar.UpperLimit());
      } else {
         unsigned int i = trafo.IntOfExt(ipar.Number());
         double err =
            st.Error().IsValid() ? std::sqrt(2. * up * st.Error().InvHessian()(i, i)) : st.Parameters().Dirin()(i);
         Add(ipar.GetName(), st.Vec()(i), err);
      }
   }
 
   // need to be set afterwards because becore the ::Add method set fCovarianceValid to false
   fCovarianceValid = st.Error().IsValid();
 
   fCovStatus = -1; // when not available
   // if (st.Error().HesseFailed() || st.Error().InvertFailed() ) fCovStatus = -1;
   // when available
   if (st.Error().IsAvailable())
      fCovStatus = 0;
 
   if (fCovarianceValid) {
      fCovariance = trafo.Int2extCovariance(st.Vec(), st.Error().InvHessian());
      fIntCovariance =
         MnUserCovariance({st.Error().InvHessian().Data(), st.Error().InvHessian().size()}, st.Error().InvHessian().Nrow());
      fCovariance.Scale(2. * up);
 
      assert(fCovariance.Nrow() == VariableParameters());
 
      if (!st.Error().IsNotPosDef())
         fCovStatus = 1; // when is valid and not NotPosDef
   }
   if (st.Error().IsMadePosDef())
      fCovStatus = 2;
   if (st.Error().IsAccurate())
      fCovStatus = 3;
 
}
 
MnUserCovariance MnUserParameterState::Hessian() const
{
   // invert covariance matrix and return Hessian
   // need to copy in a MnSymMatrix
   MnPrint print("MnUserParameterState::Hessian");
 
   MnAlgebraicSymMatrix mat(fCovariance.Nrow());
   std::copy(fCovariance.Data().begin(), fCovariance.Data().end(), mat.Data());
   int ifail = Invert(mat);
   if (ifail != 0) {
      print.Warn("Inversion failed; return diagonal matrix");
 
      MnUserCovariance tmp(fCovariance.Nrow());
      for (unsigned int i = 0; i < fCovariance.Nrow(); i++) {
         tmp(i, i) = 1. / fCovariance(i, i);
      }
      return tmp;
   }
 
   MnUserCovariance hessian(mat.Data(), fCovariance.Nrow());
   return hessian;
}
 
// facade: forward interface of MnUserParameters and MnUserTransformation
// via MnUserParameterState
 
const std::vector<MinuitParameter> &MnUserParameterState::MinuitParameters() const
{
   // access to parameters (row-wise)
   return fParameters.Parameters();
}
 
std::vector<double> MnUserParameterState::Params() const
{
   // access to parameters in column-wise representation
   return fParameters.Params();
}
std::vector<double> MnUserParameterState::Errors() const
{
   // access to errors in column-wise representation
   return fParameters.Errors();
}
 
const MinuitParameter &MnUserParameterState::Parameter(unsigned int i) const
{
   // access to single Parameter i
   return fParameters.Parameter(i);
}
 
void MnUserParameterState::Add(const std::string &name, double val, double err)
{
   MnPrint print("MnUserParameterState::Add");
 
   // add free Parameter
   if (fParameters.Add(name, val, err)) {
      fIntParameters.push_back(val);
      fCovarianceValid = false;
      fValid = true;
   } else {
      // redefine an existing parameter
      int i = Index(name);
      SetValue(i, val);
      if (Parameter(i).IsConst()) {
         print.Warn("Cannot modify status of constant parameter", name);
         return;
      }
      SetError(i, err);
      // release if it was fixed
      if (Parameter(i).IsFixed())
         Release(i);
   }
}
 
void MnUserParameterState::Add(const std::string &name, double val, double err, double low, double up)
{
   MnPrint print("MnUserParameterState::Add");
 
   // add limited Parameter
   if (fParameters.Add(name, val, err, low, up)) {
      fCovarianceValid = false;
      fIntParameters.push_back(Ext2int(Index(name), val));
      fValid = true;
   } else { // Parameter already exist - just set values
      int i = Index(name);
      SetValue(i, val);
      if (Parameter(i).IsConst()) {
         print.Warn("Cannot modify status of constant parameter", name);
         return;
      }
      SetError(i, err);
      SetLimits(i, low, up);
      // release if it was fixed
      if (Parameter(i).IsFixed())
         Release(i);
   }
}
 
void MnUserParameterState::Add(const std::string &name, double val)
{
   // add const Parameter
   if (fParameters.Add(name, val))
      fValid = true;
   else
      SetValue(name, val);
}
 
void MnUserParameterState::AddCovariance(const MnUserCovariance &cov)
{
 
   unsigned int nrow = VariableParameters();
   assert(cov.Nrow() >= nrow);
 
   // add external covariance matrix
   fCovariance = cov;
 
   // compute and add internal covariance matrix
   MnUserCovariance covsqueezed;
   if (cov.Nrow() > nrow)
      covsqueezed = MnCovarianceSqueeze()(cov, nrow);
   else if (cov.Nrow() == nrow)
      covsqueezed = cov;
 
   MnAlgebraicVector params(nrow);
   for (unsigned int i = 0; i < nrow; i++)
      params(i) = fParameters.Params()[i];
 
   MnAlgebraicSymMatrix covmat(nrow);
   for (unsigned int i = 0; i < nrow; i++)
      for (unsigned int j = i; j < nrow; j++)
         covmat(i, j) = covsqueezed(i, j);
 
   fIntCovariance = fParameters.Trafo().Ext2intCovariance(params, covmat);
   fIntCovariance.Scale(0.5); //?
   fCovarianceValid = true;
   fCovStatus = 0; //?
}
 
// interaction via external number of Parameter
 
void MnUserParameterState::Fix(unsigned int e)
{
   // fix parameter e (external index)
   if (!Parameter(e).IsFixed() && !Parameter(e).IsConst()) {
      unsigned int i = IntOfExt(e);
      if (fCovarianceValid) {
         // when covariance is valid remove fixed parameter row/column in covariance matrix
         // use squeeze to remove first from Hessian and obtain then new covariance
         fCovariance = MnCovarianceSqueeze()(fCovariance, i);
         fIntCovariance = MnCovarianceSqueeze()(fIntCovariance, i);
      }
      fIntParameters.erase(fIntParameters.begin() + i, fIntParameters.begin() + i + 1);
   }
   fParameters.Fix(e);
}
 
void MnUserParameterState::Release(unsigned int e)
{
   // release parameter e (external index)
   // no-op if parameter is const or if it is not fixed
   if (Parameter(e).IsConst() || !Parameter(e).IsFixed())
      return;
   fParameters.Release(e);
   fCovarianceValid = false;
   unsigned int i = IntOfExt(e);
   if (Parameter(e).HasLimits())
      fIntParameters.insert(fIntParameters.begin() + i, Ext2int(e, Parameter(e).Value()));
   else
      fIntParameters.insert(fIntParameters.begin() + i, Parameter(e).Value());
}
 
void MnUserParameterState::SetValue(unsigned int e, double val)
{
   // set value for parameter e ( external index )
   fParameters.SetValue(e, val);
   if (!Parameter(e).IsFixed() && !Parameter(e).IsConst()) {
      unsigned int i = IntOfExt(e);
      if (Parameter(e).HasLimits())
         fIntParameters[i] = Ext2int(e, val);
      else
         fIntParameters[i] = val;
   }
}
 
void MnUserParameterState::SetError(unsigned int e, double val)
{
   // set error for  parameter e (external index)
   fParameters.SetError(e, val);
}
 
void MnUserParameterState::SetLimits(unsigned int e, double low, double up)
{
   // set limits for parameter e (external index)
   fParameters.SetLimits(e, low, up);
   fCovarianceValid = false;
   if (!Parameter(e).IsFixed() && !Parameter(e).IsConst()) {
      unsigned int i = IntOfExt(e);
      if (low < fIntParameters[i] && fIntParameters[i] < up)
         fIntParameters[i] = Ext2int(e, fIntParameters[i]);
      else if (low >= fIntParameters[i])
         fIntParameters[i] = Ext2int(e, low + 0.1 * Parameter(e).Error());
      else
         fIntParameters[i] = Ext2int(e, up - 0.1 * Parameter(e).Error());
   }
}
 
void MnUserParameterState::SetUpperLimit(unsigned int e, double up)
{
   // set upper limit for parameter e (external index)
   fParameters.SetUpperLimit(e, up);
   fCovarianceValid = false;
   if (!Parameter(e).IsFixed() && !Parameter(e).IsConst()) {
      unsigned int i = IntOfExt(e);
      if (fIntParameters[i] < up)
         fIntParameters[i] = Ext2int(e, fIntParameters[i]);
      else
         fIntParameters[i] = Ext2int(e, up - 0.1 * Parameter(e).Error());
   }
}
 
void MnUserParameterState::SetLowerLimit(unsigned int e, double low)
{
   // set lower limit for parameter e (external index)
   fParameters.SetLowerLimit(e, low);
   fCovarianceValid = false;
   if (!Parameter(e).IsFixed() && !Parameter(e).IsConst()) {
      unsigned int i = IntOfExt(e);
      if (low < fIntParameters[i])
         fIntParameters[i] = Ext2int(e, fIntParameters[i]);
      else
         fIntParameters[i] = Ext2int(e, low + 0.1 * Parameter(e).Error());
   }
}
 
void MnUserParameterState::RemoveLimits(unsigned int e)
{
   // remove limit for parameter e (external index)
   fParameters.RemoveLimits(e);
   fCovarianceValid = false;
   if (!Parameter(e).IsFixed() && !Parameter(e).IsConst())
      fIntParameters[IntOfExt(e)] = Value(e);
}
 
double MnUserParameterState::Value(unsigned int i) const
{
   // get value for parameter e (external index)
   return fParameters.Value(i);
}
double MnUserParameterState::Error(unsigned int i) const
{
   // get error for parameter e (external index)
   return fParameters.Error(i);
}
 
// interaction via name of Parameter
 
void MnUserParameterState::Fix(const std::string &name)
{
   Fix(Index(name));
}
 
void MnUserParameterState::Release(const std::string &name)
{
   Release(Index(name));
}
 
void MnUserParameterState::SetValue(const std::string &name, double val)
{
   SetValue(Index(name), val);
}
 
void MnUserParameterState::SetError(const std::string &name, double val)
{
   SetError(Index(name), val);
}
 
void MnUserParameterState::SetLimits(const std::string &name, double low, double up)
{
   SetLimits(Index(name), low, up);
}
 
void MnUserParameterState::SetUpperLimit(const std::string &name, double up)
{
   SetUpperLimit(Index(name), up);
}
 
void MnUserParameterState::SetLowerLimit(const std::string &name, double low)
{
   SetLowerLimit(Index(name), low);
}
 
void MnUserParameterState::RemoveLimits(const std::string &name)
{
   RemoveLimits(Index(name));
}
 
double MnUserParameterState::Value(const std::string &name) const
{
   return Value(Index(name));
}
 
double MnUserParameterState::Error(const std::string &name) const
{
   return Error(Index(name));
}
 
unsigned int MnUserParameterState::Index(const std::string &name) const
{
   // convert name into external number of Parameter
   return fParameters.Index(name);
}
 
const char *MnUserParameterState::Name(unsigned int i) const
{
   // convert external number into name of Parameter (API returning a const char *)
   return fParameters.Name(i);
}
const std::string &MnUserParameterState::GetName(unsigned int i) const
{
   // convert external number into name of Parameter (new interface returning a string)
   return fParameters.GetName(i);
}
 
// transformation internal <-> external (forward to transformation class)
 
double MnUserParameterState::Int2ext(unsigned int i, double val) const
{
   // internal to external value
   return fParameters.Trafo().Int2ext(i, val);
}
double MnUserParameterState::Ext2int(unsigned int e, double val) const
{
   // external  to internal value
   return fParameters.Trafo().Ext2int(e, val);
}
unsigned int MnUserParameterState::IntOfExt(unsigned int ext) const
{
   // return internal index for external index ext
   return fParameters.Trafo().IntOfExt(ext);
}
unsigned int MnUserParameterState::ExtOfInt(unsigned int internal) const
{
   // return external index for internal index internal
   return fParameters.Trafo().ExtOfInt(internal);
}
unsigned int MnUserParameterState::VariableParameters() const
{
   // return number of variable parameters
   return fParameters.Trafo().VariableParameters();
}
const MnMachinePrecision &MnUserParameterState::Precision() const
{
   // return global parameter precision
   return fParameters.Precision();
}
 
void MnUserParameterState::SetPrecision(double eps)
{
   // set global parameter precision
   fParameters.SetPrecision(eps);
}
 
MnGlobalCorrelationCoeff MnUserParameterState::GlobalCC() const
{
   if (!fCovarianceValid) {
      return MnGlobalCorrelationCoeff{}; // invalidate
   }
   int n = fIntCovariance.Nrow();
   MnAlgebraicSymMatrix invHessian{static_cast<unsigned int>(n)};
   for (int i = 0; i < n; ++i) {
      for (int j = 0; j <= i; ++j) {
         invHessian(i, j) = fIntCovariance(i, j);
      }
   }
   return MnGlobalCorrelationCoeff{invHessian};
}
 
} // namespace Minuit2
 
} // namespace ROOT