Back to home page

EIC code displayed by LXR
 
 

    

Warning, file /include/root/Fit/FitUtil.h was not indexed or was modified since last indexation (in which case cross-reference links may be missing, inaccurate or erroneous).

0001 // @(#)root/mathcore:$Id$
0002 // Author: L. Moneta Tue Nov 28 10:52:47 2006
0003 
0004 /**********************************************************************
0005  *                                                                    *
0006  * Copyright (c) 2006  LCG ROOT Math Team, CERN/PH-SFT                *
0007  *                                                                    *
0008  *                                                                    *
0009  **********************************************************************/
0010 
0011 // Header file for class FitUtil
0012 
0013 #ifndef ROOT_Fit_FitUtil
0014 #define ROOT_Fit_FitUtil
0015 
0016 #include "Math/IParamFunctionfwd.h"
0017 #include "Math/IParamFunction.h"
0018 
0019 #include "Fit/BinData.h"
0020 #include "Fit/UnBinData.h"
0021 #include "ROOT/EExecutionPolicy.hxx"
0022 
0023 #include "Math/Integrator.h"
0024 #include "Math/IntegratorMultiDim.h"
0025 
0026 #include "TError.h"
0027 #include <vector>
0028 
0029 // using parameter cache is not thread safe but needed for normalizing the functions
0030 #define USE_PARAMCACHE
0031 
0032 //#define DEBUG_FITUTIL
0033 
0034 namespace ROOT {
0035 
0036 namespace Fit {
0037 
0038 /**
0039    namespace defining utility free functions using in Fit for evaluating the various fit method
0040    functions (chi2, likelihood, etc..)  given the data and the model function
0041 
0042    @ingroup FitMain
0043 */
0044 namespace FitUtil {
0045 
0046   typedef  ROOT::Math::IParamMultiFunction IModelFunction;
0047   typedef  ROOT::Math::IParamMultiGradFunction IGradModelFunction;
0048 
0049   template <class T>
0050   using IGradModelFunctionTempl = ROOT::Math::IParamMultiGradFunctionTempl<T>;
0051 
0052   template <class T>
0053   using IModelFunctionTempl = ROOT::Math::IParamMultiFunctionTempl<T>;
0054 
0055   // internal class defining
0056   template <class T>
0057   class LikelihoodAux {
0058   public:
0059      LikelihoodAux(T logv = {}, T w = {}, T w2 = {}) : logvalue(logv), weight(w), weight2(w2) {}
0060 
0061      LikelihoodAux operator+(const LikelihoodAux &l) const
0062      {
0063         return LikelihoodAux<T>(logvalue + l.logvalue, weight + l.weight, weight2 + l.weight2);
0064      }
0065 
0066      LikelihoodAux &operator+=(const LikelihoodAux &l)
0067      {
0068         logvalue += l.logvalue;
0069         weight += l.weight;
0070         weight2 += l.weight2;
0071         return *this;
0072      }
0073 
0074      T logvalue;
0075      T weight;
0076      T weight2;
0077   };
0078 
0079   template <>
0080   class LikelihoodAux<double> {
0081   public:
0082      LikelihoodAux(double logv = 0.0, double w = 0.0, double w2 = 0.0) : logvalue(logv), weight(w), weight2(w2){};
0083 
0084      LikelihoodAux operator+(const LikelihoodAux &l) const
0085      {
0086         return LikelihoodAux<double>(logvalue + l.logvalue, weight + l.weight, weight2 + l.weight2);
0087      }
0088 
0089      LikelihoodAux &operator+=(const LikelihoodAux &l)
0090      {
0091         logvalue += l.logvalue;
0092         weight += l.weight;
0093         weight2 += l.weight2;
0094         return *this;
0095      }
0096 
0097      double logvalue;
0098      double weight;
0099      double weight2;
0100   };
0101 
0102   // internal class to evaluate the function or the integral
0103   // and cached internal integration details
0104   // if useIntegral is false no allocation is done
0105   // and this is a dummy class
0106   // class is templated on any parametric functor implementing operator()(x,p) and NDim()
0107   // contains a constant pointer to the function
0108 
0109   template <class ParamFunc = ROOT::Math::IParamMultiFunctionTempl<double>>
0110   class IntegralEvaluator {
0111 
0112   public:
0113      IntegralEvaluator(const ParamFunc &func, const double *p, bool useIntegral = true,
0114                        ROOT::Math::IntegrationOneDim::Type igType = ROOT::Math::IntegrationOneDim::kDEFAULT)
0115         : fDim(0), fParams(nullptr), fFunc(nullptr), fIg1Dim(nullptr), fIgNDim(nullptr), fFunc1Dim(nullptr), fFuncNDim(nullptr)
0116      {
0117         if (useIntegral) {
0118            SetFunction(func, p, igType);
0119         }
0120      }
0121 
0122      void SetFunction(const ParamFunc &func, const double *p = nullptr,
0123                       ROOT::Math::IntegrationOneDim::Type igType = ROOT::Math::IntegrationOneDim::kDEFAULT)
0124      {
0125         // set the integrand function and create required wrapper
0126         // to perform integral in (x) of a generic  f(x,p)
0127         fParams = p;
0128         fDim = func.NDim();
0129         // copy the function object to be able to modify the parameters
0130         // fFunc = dynamic_cast<ROOT::Math::IParamMultiFunction *>( func.Clone() );
0131         fFunc = &func;
0132         assert(fFunc != nullptr);
0133         // set parameters in function
0134         // fFunc->SetParameters(p);
0135         if (fDim == 1) {
0136            fFunc1Dim =
0137               new ROOT::Math::WrappedMemFunction<IntegralEvaluator, double (IntegralEvaluator::*)(double) const>(
0138                  *this, &IntegralEvaluator::F1);
0139            fIg1Dim = new ROOT::Math::IntegratorOneDim(igType);
0140            // fIg1Dim->SetFunction( static_cast<const ROOT::Math::IMultiGenFunction & >(*fFunc),false);
0141            fIg1Dim->SetFunction(static_cast<const ROOT::Math::IGenFunction &>(*fFunc1Dim));
0142         } else if (fDim > 1) {
0143            fFuncNDim =
0144               new ROOT::Math::WrappedMemMultiFunction<IntegralEvaluator, double (IntegralEvaluator::*)(const double *)
0145                                                                             const>(*this, &IntegralEvaluator::FN, fDim);
0146            fIgNDim = new ROOT::Math::IntegratorMultiDim();
0147            fIgNDim->SetFunction(*fFuncNDim);
0148         } else
0149            assert(fDim > 0);
0150      }
0151 
0152      void SetParameters(const double *p)
0153      {
0154         // copy just the pointer
0155         fParams = p;
0156      }
0157 
0158      ~IntegralEvaluator()
0159      {
0160         if (fIg1Dim)
0161            delete fIg1Dim;
0162         if (fIgNDim)
0163            delete fIgNDim;
0164         if (fFunc1Dim)
0165            delete fFunc1Dim;
0166         if (fFuncNDim)
0167            delete fFuncNDim;
0168         // if (fFunc) delete fFunc;
0169      }
0170 
0171      // evaluation of integrand function (one-dim)
0172      double F1(double x) const
0173      {
0174         double xx = x;
0175         return ExecFunc(fFunc, &xx, fParams);
0176      }
0177      // evaluation of integrand function (multi-dim)
0178      double FN(const double *x) const { return ExecFunc(fFunc, x, fParams); }
0179 
0180      double Integral(const double *x1, const double *x2)
0181      {
0182         // return unnormalized integral
0183         return (fIg1Dim) ? fIg1Dim->Integral(*x1, *x2) : fIgNDim->Integral(x1, x2);
0184      }
0185 
0186      double operator()(const double *x1, const double *x2)
0187      {
0188         // return normalized integral, divided by bin volume (dx1*dx...*dxn)
0189         if (fIg1Dim) {
0190            double dV = *x2 - *x1;
0191            return fIg1Dim->Integral(*x1, *x2) / dV;
0192         } else if (fIgNDim) {
0193            double dV = 1;
0194            for (unsigned int i = 0; i < fDim; ++i)
0195               dV *= (x2[i] - x1[i]);
0196            return fIgNDim->Integral(x1, x2) / dV;
0197            //                   std::cout << " do integral btw x " << x1[0] << "  " << x2[0] << " y " << x1[1] << "  "
0198            //                   << x2[1] << " dV = " << dV << " result = " << result << std::endl; return result;
0199         } else
0200            assert(1.); // should never be here
0201         return 0;
0202      }
0203 
0204   private:
0205      template <class T>
0206      inline double ExecFunc(T *f, const double *x, const double *p) const
0207      {
0208         return (*f)(x, p);
0209      }
0210 
0211 #ifdef R__HAS_STD_EXPERIMENTAL_SIMD
0212 
0213 #if __clang_major__ > 16
0214 #pragma clang diagnostic push
0215 #pragma clang diagnostic ignored "-Wvla-cxx-extension"
0216 #endif
0217 
0218      inline double ExecFunc(const IModelFunctionTempl<ROOT::Double_v> *f, const double *x, const double *p) const
0219      {
0220         ROOT::Double_v xx[fDim];
0221         for (unsigned int i = 0; i < fDim; ++i) {
0222            xx[i][0] = x[i];
0223            for (std::size_t j = 1; j < ROOT::Double_v::size(); ++j) {
0224               xx[i][j] = 0.0;
0225            }
0226         }
0227         auto res = (*f)(xx, p);
0228         return res[0];
0229      }
0230 
0231 #if __clang_major__ > 16
0232 #pragma clang diagnostic pop
0233 #endif
0234 
0235 #endif
0236 
0237      // objects of this class are not meant to be copied / assigned
0238      IntegralEvaluator(const IntegralEvaluator &rhs) = delete;
0239      IntegralEvaluator &operator=(const IntegralEvaluator &rhs) = delete;
0240 
0241      unsigned int fDim;
0242      const double *fParams;
0243      // ROOT::Math::IParamMultiFunction * fFunc;  // copy of function in order to be able to change parameters
0244      // const ParamFunc * fFunc;       //  reference to a generic parametric function
0245      const ParamFunc *fFunc;
0246      ROOT::Math::IntegratorOneDim *fIg1Dim;
0247      ROOT::Math::IntegratorMultiDim *fIgNDim;
0248      ROOT::Math::IGenFunction *fFunc1Dim;
0249      ROOT::Math::IMultiGenFunction *fFuncNDim;
0250   };
0251 
0252   /** Chi2 Functions */
0253 
0254   /**
0255       evaluate the Chi2 given a model function and the data at the point x.
0256       return also nPoints as the effective number of used points in the Chi2 evaluation
0257   */
0258   double EvaluateChi2(const IModelFunction &func, const BinData &data, const double *p, unsigned int &nPoints,
0259                       ::ROOT::EExecutionPolicy executionPolicy, unsigned nChunks = 0);
0260 
0261   /**
0262       evaluate the effective Chi2 given a model function and the data at the point x.
0263       The effective chi2 uses the errors on the coordinates : W = 1/(sigma_y**2 + ( sigma_x_i * df/dx_i )**2 )
0264       return also nPoints as the effective number of used points in the Chi2 evaluation
0265   */
0266   double EvaluateChi2Effective(const IModelFunction &func, const BinData &data, const double *x, unsigned int &nPoints);
0267 
0268   /**
0269       evaluate the Chi2 gradient given a model function and the data at the point p.
0270       return also nPoints as the effective number of used points in the Chi2 evaluation
0271   */
0272   void EvaluateChi2Gradient(const IModelFunction &func, const BinData &data, const double *p, double *grad,
0273                             unsigned int &nPoints,
0274                             ::ROOT::EExecutionPolicy executionPolicy = ::ROOT::EExecutionPolicy::kSequential,
0275                             unsigned nChunks = 0);
0276 
0277   /**
0278       evaluate the LogL given a model function and the data at the point x.
0279       return also nPoints as the effective number of used points in the LogL evaluation
0280   */
0281   double EvaluateLogL(const IModelFunction &func, const UnBinData &data, const double *p, int iWeight, bool extended,
0282                       unsigned int &nPoints, ::ROOT::EExecutionPolicy executionPolicy, unsigned nChunks = 0);
0283 
0284   /**
0285       evaluate the LogL gradient given a model function and the data at the point p.
0286       return also nPoints as the effective number of used points in the LogL evaluation
0287   */
0288   void EvaluateLogLGradient(const IModelFunction &func, const UnBinData &data, const double *p, double *grad,
0289                             unsigned int &nPoints,
0290                             ::ROOT::EExecutionPolicy executionPolicy = ::ROOT::EExecutionPolicy::kSequential,
0291                             unsigned nChunks = 0);
0292 
0293   // #ifdef R__HAS_STD_EXPERIMENTAL_SIMD
0294   //    template <class NotCompileIfScalarBackend = std::enable_if<!(std::is_same<double, ROOT::Double_v>::value)>>
0295   //    void EvaluateLogLGradient(const IModelFunctionTempl<ROOT::Double_v> &, const UnBinData &, const double *, double
0296   //    *, unsigned int & ) {}
0297   // #endif
0298 
0299   /**
0300       evaluate the Poisson LogL given a model function and the data at the point p.
0301       return also nPoints as the effective number of used points in the LogL evaluation
0302       By default is extended, pass extend to false if want to be not extended (MultiNomial)
0303   */
0304   double EvaluatePoissonLogL(const IModelFunction &func, const BinData &data, const double *p, int iWeight,
0305                              bool extended, unsigned int &nPoints, ::ROOT::EExecutionPolicy executionPolicy,
0306                              unsigned nChunks = 0);
0307 
0308   /**
0309       evaluate the Poisson LogL given a model function and the data at the point p.
0310       return also nPoints as the effective number of used points in the LogL evaluation
0311   */
0312   void EvaluatePoissonLogLGradient(const IModelFunction &func, const BinData &data, const double *p, double *grad,
0313                                    unsigned int &nPoints,
0314                                    ::ROOT::EExecutionPolicy executionPolicy = ::ROOT::EExecutionPolicy::kSequential,
0315                                    unsigned nChunks = 0);
0316 
0317   // methods required by dedicate minimizer like Fumili
0318 
0319   /**
0320       evaluate the residual contribution to the Chi2 given a model function and the BinPoint data
0321       and if the pointer g is not null evaluate also the gradient of the residual.
0322       If the function provides parameter derivatives they are used otherwise a simple derivative calculation
0323       is used
0324   */
0325   double EvaluateChi2Residual(const IModelFunction &func, const BinData &data, const double *p, unsigned int ipoint,
0326                               double *g = nullptr, double * h = nullptr, bool hasGrad = false, bool fullHessian = false);
0327 
0328   /**
0329       evaluate the pdf contribution to the LogL given a model function and the BinPoint data.
0330       If the pointer g is not null evaluate also the gradient of the pdf.
0331       If the function provides parameter derivatives they are used otherwise a simple derivative calculation
0332       is used
0333   */
0334   double
0335   EvaluatePdf(const IModelFunction &func, const UnBinData &data, const double *p, unsigned int ipoint, double *g = nullptr, double * h = nullptr, bool hasGrad = false, bool fullHessian = false);
0336 
0337 
0338    /**
0339        evaluate the pdf contribution to the Poisson LogL given a model function and the BinPoint data.
0340        If the pointer g is not null evaluate also the gradient of the Poisson pdf.
0341        If the function provides parameter derivatives they are used otherwise a simple derivative calculation
0342        is used
0343    */
0344    double EvaluatePoissonBinPdf(const IModelFunction & func, const BinData & data, const double * x, unsigned int ipoint, double * g = nullptr, double * h = nullptr, bool hasGrad = false, bool fullHessian = false);
0345 
0346    unsigned setAutomaticChunking(unsigned nEvents);
0347 
0348    template <class T>
0349    struct Evaluate {};
0350 
0351 #ifdef R__HAS_STD_EXPERIMENTAL_SIMD
0352    template <>
0353    struct Evaluate<Double_v> {
0354       static double EvalChi2(const IModelFunctionTempl<Double_v> &func, const BinData &data, const double *p,
0355                              unsigned int &nPoints, ::ROOT::EExecutionPolicy executionPolicy, unsigned nChunks = 0);
0356 
0357       static double EvalLogL(const IModelFunctionTempl<Double_v> &func, const UnBinData &data, const double *const p,
0358                              int iWeight, bool extended, unsigned int &nPoints,
0359                              ::ROOT::EExecutionPolicy executionPolicy, unsigned nChunks = 0);
0360 
0361       static double EvalPoissonLogL(const IModelFunctionTempl<Double_v> &func, const BinData &data, const double *p,
0362                                     int iWeight, bool extended, unsigned int, ::ROOT::EExecutionPolicy executionPolicy,
0363                                     unsigned nChunks = 0);
0364 
0365       static double
0366       EvalChi2Effective(const IModelFunctionTempl<Double_v> &, const BinData &, const double *, unsigned int &)
0367       {
0368          Error("FitUtil::Evaluate<T>::EvalChi2Effective", "The vectorized evaluation of the Chi2 with coordinate errors is still not supported");
0369          return -1.;
0370       }
0371 
0372       static void EvalChi2Gradient(const IModelFunctionTempl<Double_v> &f, const BinData &data, const double *p,
0373                                    double *grad, unsigned int &nPoints,
0374                                    ::ROOT::EExecutionPolicy executionPolicy = ::ROOT::EExecutionPolicy::kSequential,
0375                                    unsigned nChunks = 0);
0376 
0377       static double EvalChi2Residual(const IModelFunctionTempl<Double_v> &, const BinData &, const double *,
0378                                      unsigned int, double *, double *, bool, bool)
0379       {
0380          Error("FitUtil::Evaluate<T>::EvalChi2Residual", "The vectorized evaluation of the Chi2 with the ith residual is still not supported");
0381          return -1.;
0382       }
0383 
0384       /// evaluate the pdf (Poisson) contribution to the logl (return actually log of pdf)
0385       /// and its gradient
0386       static double EvalPoissonBinPdf(const IModelFunctionTempl<Double_v> &, const BinData &, const double *,
0387                                       unsigned int, double *, double *, bool, bool)
0388       {
0389          Error("FitUtil::Evaluate<T>::EvaluatePoissonBinPdf", "The vectorized evaluation of the BinnedLikelihood fit evaluated point by point is still not supported");
0390          return -1.;
0391       }
0392 
0393       static double EvalPdf(const IModelFunctionTempl<Double_v> &, const UnBinData &, const double *, unsigned int,
0394                             double *, double *, bool, bool)
0395       {
0396          Error("FitUtil::Evaluate<T>::EvalPdf", "The vectorized evaluation of the LogLikelihood fit evaluated point by point is still not supported");
0397          return -1.;
0398       }
0399 
0400       static void
0401       EvalPoissonLogLGradient(const IModelFunctionTempl<Double_v> &f, const BinData &data, const double *p,
0402                               double *grad, unsigned int &,
0403                               ::ROOT::EExecutionPolicy executionPolicy = ::ROOT::EExecutionPolicy::kSequential,
0404                               unsigned nChunks = 0);
0405 
0406       static void EvalLogLGradient(const IModelFunctionTempl<Double_v> &f, const UnBinData &data, const double *p,
0407                                    double *grad, unsigned int &,
0408                                    ::ROOT::EExecutionPolicy executionPolicy = ::ROOT::EExecutionPolicy::kSequential,
0409                                    unsigned nChunks = 0);
0410    };
0411 #endif // R__HAS_STD_EXPERIMENTAL_SIMD
0412 
0413    template <>
0414    struct Evaluate<double> {
0415 
0416       static double EvalChi2(const IModelFunction &func, const BinData &data, const double *p, unsigned int &nPoints,
0417                              ::ROOT::EExecutionPolicy executionPolicy, unsigned nChunks = 0)
0418       {
0419          // evaluate the chi2 given a  function reference, the data and returns the value and also in nPoints
0420          // the actual number of used points
0421          // normal chi2 using only error on values (from fitting histogram)
0422          // optionally the integral of function in the bin is used
0423 
0424 
0425          //Info("EvalChi2","Using non-vectorized implementation %d",(int) data.Opt().fIntegral);
0426 
0427          return FitUtil::EvaluateChi2(func, data, p, nPoints, executionPolicy, nChunks);
0428       }
0429 
0430       static double EvalLogL(const IModelFunctionTempl<double> &func, const UnBinData &data, const double *p,
0431                              int iWeight, bool extended, unsigned int &nPoints,
0432                              ::ROOT::EExecutionPolicy executionPolicy, unsigned nChunks = 0)
0433       {
0434          return FitUtil::EvaluateLogL(func, data, p, iWeight, extended, nPoints, executionPolicy, nChunks);
0435       }
0436 
0437       static double EvalPoissonLogL(const IModelFunctionTempl<double> &func, const BinData &data, const double *p,
0438                                     int iWeight, bool extended, unsigned int &nPoints,
0439                                     ::ROOT::EExecutionPolicy executionPolicy, unsigned nChunks = 0)
0440       {
0441          return FitUtil::EvaluatePoissonLogL(func, data, p, iWeight, extended, nPoints, executionPolicy, nChunks);
0442       }
0443 
0444       static double EvalChi2Effective(const IModelFunctionTempl<double> &func, const BinData & data, const double * p, unsigned int &nPoints)
0445       {
0446          return FitUtil::EvaluateChi2Effective(func, data, p, nPoints);
0447       }
0448       static void EvalChi2Gradient(const IModelFunctionTempl<double> &func, const BinData &data, const double *p,
0449                                    double *g, unsigned int &nPoints,
0450                                    ::ROOT::EExecutionPolicy executionPolicy = ::ROOT::EExecutionPolicy::kSequential,
0451                                    unsigned nChunks = 0)
0452       {
0453          FitUtil::EvaluateChi2Gradient(func, data, p, g, nPoints, executionPolicy, nChunks);
0454       }
0455 
0456       static double EvalChi2Residual(const IModelFunctionTempl<double> &func, const BinData & data, const double * p, unsigned int i, double *g, double * h,
0457                                     bool hasGrad, bool fullHessian)
0458       {
0459          return FitUtil::EvaluateChi2Residual(func, data, p, i, g, h, hasGrad, fullHessian);
0460       }
0461 
0462       /// evaluate the pdf (Poisson) contribution to the logl (return actually log of pdf)
0463       /// and its gradient
0464       static double EvalPoissonBinPdf(const IModelFunctionTempl<double> &func, const BinData & data, const double *p, unsigned int i, double *g, double * h, bool hasGrad, bool fullHessian) {
0465          return FitUtil::EvaluatePoissonBinPdf(func, data, p, i, g, h, hasGrad, fullHessian);
0466       }
0467 
0468       static double EvalPdf(const IModelFunctionTempl<double> &func, const UnBinData & data, const double *p, unsigned int i, double *g, double * h, bool hasGrad, bool fullHessian) {
0469          return FitUtil::EvaluatePdf(func, data, p, i, g, h, hasGrad, fullHessian);
0470       }
0471 
0472       static void
0473       EvalPoissonLogLGradient(const IModelFunctionTempl<double> &func, const BinData &data, const double *p, double *g,
0474                               unsigned int &nPoints,
0475                               ::ROOT::EExecutionPolicy executionPolicy = ::ROOT::EExecutionPolicy::kSequential,
0476                               unsigned nChunks = 0)
0477       {
0478          FitUtil::EvaluatePoissonLogLGradient(func, data, p, g, nPoints, executionPolicy, nChunks);
0479       }
0480 
0481       static void EvalLogLGradient(const IModelFunctionTempl<double> &func, const UnBinData &data, const double *p,
0482                                    double *g, unsigned int &nPoints,
0483                                    ::ROOT::EExecutionPolicy executionPolicy = ::ROOT::EExecutionPolicy::kSequential,
0484                                    unsigned nChunks = 0)
0485       {
0486          FitUtil::EvaluateLogLGradient(func, data, p, g, nPoints, executionPolicy, nChunks);
0487       }
0488    };
0489 
0490    } // end namespace FitUtil
0491 
0492    } // end namespace Fit
0493 
0494    } // end namespace ROOT
0495 
0496 #endif /* ROOT_Fit_FitUtil */
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!