root/Fit/PoissonLikelihoodFCN.h

0001 // @(#)root/mathcore:$Id$
0002 // Author: L. Moneta Fri Aug 17 14:29:24 2007
0003
0004 /**********************************************************************
0005  *                                                                    *
0006  * Copyright (c) 2007  LCG ROOT Math Team, CERN/PH-SFT                *
0007  *                                                                    *
0008  *                                                                    *
0009  **********************************************************************/
0010
0011 // Header file for class PoissonLikelihoodFCN
0012
0013 #ifndef ROOT_Fit_PoissonLikelihoodFCN
0014 #define ROOT_Fit_PoissonLikelihoodFCN
0015
0016 #include "ROOT/EExecutionPolicy.hxx"
0017 #include "Fit/BasicFCN.h"
0018 #include "Fit/BinData.h"
0019 #include "Fit/FitUtil.h"
0020 #include "Math/IParamFunction.h"
0021
0022 #include <memory>
0023 #include <vector>
0024
0025 //#define PARALLEL
0026 // #ifdef PARALLEL
0027 // #ifndef ROOT_Fit_FitUtilParallel
0028 // #include "Fit/FitUtilParallel.h"
0029 // #endif
0030 // #endif
0031
0032 namespace ROOT {
0033
0034    namespace Fit {
0035
0036
0037 //___________________________________________________________________________________
0038 /**
0039    class evaluating the log likelihood
0040    for binned Poisson likelihood fits
0041    it is template to distinguish gradient and non-gradient case
0042
0043    @ingroup  FitMethodFunc
0044 */
0045 template<class DerivFunType, class ModelFunType = ROOT::Math::IParamMultiFunction>
0046 class PoissonLikelihoodFCN : public BasicFCN<DerivFunType,ModelFunType,BinData>  {
0047
0048 public:
0049    typedef typename ModelFunType::BackendType T;
0050    typedef  BasicFCN<DerivFunType,ModelFunType,BinData> BaseFCN;
0051
0052    typedef  ::ROOT::Math::BasicFitMethodFunction<DerivFunType> BaseObjFunction;
0053    typedef typename  BaseObjFunction::BaseFunction BaseFunction;
0054
0055    typedef ::ROOT::Math::IParamMultiFunctionTempl<T> IModelFunction;
0056    typedef typename BaseObjFunction::Type_t Type_t;
0057
0058    /**
0059       Constructor from unbin data set and model function (pdf)
0060    */
0061    PoissonLikelihoodFCN (const std::shared_ptr<BinData> & data, const std::shared_ptr<IModelFunction> & func, int weight = 0, bool extended = true, const ::ROOT::EExecutionPolicy &executionPolicy = ::ROOT::EExecutionPolicy::kSequential ) :
0062       BaseFCN( data, func),
0063       fIsExtended(extended),
0064       fWeight(weight),
0065       fNEffPoints(0),
0066       fGrad ( std::vector<double> ( func->NPar() ) ),
0067       fExecutionPolicy(executionPolicy)
0068    { }
0069
0070    /**
0071       Constructor from unbin data set and model function (pdf) managed by the users
0072    */
0073    PoissonLikelihoodFCN (const BinData & data, const IModelFunction & func, int weight = 0, bool extended = true, const ::ROOT::EExecutionPolicy &executionPolicy = ::ROOT::EExecutionPolicy::kSequential ) :
0074       BaseFCN(std::make_shared<BinData>(data), std::shared_ptr<IModelFunction>(dynamic_cast<IModelFunction*>(func.Clone() ) ) ),
0075       fIsExtended(extended),
0076       fWeight(weight),
0077       fNEffPoints(0),
0078       fGrad ( std::vector<double> ( func.NPar() ) ),
0079       fExecutionPolicy(executionPolicy)
0080    { }
0081
0082
0083    /**
0084       Destructor (no operations)
0085    */
0086    virtual ~PoissonLikelihoodFCN () {}
0087
0088    /**
0089       Copy constructor
0090    */
0091    PoissonLikelihoodFCN(const PoissonLikelihoodFCN & f) :
0092       BaseFCN(f.DataPtr(), f.ModelFunctionPtr() ),
0093       fIsExtended(f.fIsExtended ),
0094       fWeight( f.fWeight ),
0095       fNEffPoints( f.fNEffPoints ),
0096       fGrad( f.fGrad),
0097       fExecutionPolicy(f.fExecutionPolicy)
0098    {  }
0099
0100    /**
0101       Assignment operator
0102    */
0103    PoissonLikelihoodFCN & operator = (const PoissonLikelihoodFCN & rhs) {
0104       SetData(rhs.DataPtr() );
0105       SetModelFunction(rhs.ModelFunctionPtr() );
0106       fNEffPoints = rhs.fNEffPoints;
0107       fGrad = rhs.fGrad;
0108       fIsExtended = rhs.fIsExtended;
0109       fWeight = rhs.fWeight;
0110       fExecutionPolicy = rhs.fExecutionPolicy;
0111    }
0112
0113
0114    /// clone the function (need to return Base for Windows)
0115    virtual BaseFunction * Clone() const { return new  PoissonLikelihoodFCN(*this); }
0116
0117    // effective points used in the fit
0118    virtual unsigned int NFitPoints() const { return fNEffPoints; }
0119
0120    /// i-th likelihood element and its gradient
0121    virtual double DataElement(const double * x, unsigned int i, double * g, double * h, bool fullHessian) const {
0122       if (i==0) this->UpdateNCalls();
0123       return FitUtil::Evaluate<T>::EvalPoissonBinPdf(BaseFCN::ModelFunction(), BaseFCN::Data(), x, i, g, h, BaseFCN::IsAGradFCN(), fullHessian);
0124    }
0125
0126    /// evaluate gradient
0127    virtual void Gradient(const double *x, double *g) const
0128    {
0129       // evaluate the Poisson gradient
0130       FitUtil::Evaluate<typename BaseFCN::T>::EvalPoissonLogLGradient(BaseFCN::ModelFunction(), BaseFCN::Data(), x, g,
0131                                                                       fNEffPoints, fExecutionPolicy);
0132    }
0133
0134    /**
0135     * Computes the full Hessian. Return false if Hessian is not supported
0136     */
0137    // virtual bool Hessian(const double * x, double * hess) const {
0138    //    //return full Hessian
0139    //    unsigned int np = NPoints();
0140    //    unsigned int ndim = NDim();
0141    //    unsigned int nh = ndim*(ndim+1)/2;
0142    //    for (unsigned int k = 0; k < nh;  ++k) {
0143    //       hess[k] = 0;
0144    //    }
0145    //    std::vector<double> g(np);  // gradient of the residual (y-f)/sigma
0146    //    std::vector<double> h(nh);  // hessian of residual
0147    //    for (unsigned int i = 0; i < np; i++) {
0148    //       double f = DataElement(x,i,g,h,true);
0149    //       if (f == std::numeric_limits<double>::quiet_NaN() ) return false;
0150    //       for (unsigned int j = 0; j < np; j++) {
0151    //          for (unsigned int k = 0; k <=j; k++) {
0152    //             unsigned int index = j + k * (k + 1) / 2;
0153    //             h[index] += 2. * ( g[j]*g[k] + f * h[index]*h[index] );
0154    //          }
0155    //       }
0156    //    }
0157    //    return true;
0158    // }
0159
0160    /// get type of fit method function
0161    virtual  typename BaseObjFunction::Type_t Type() const { return BaseObjFunction::kPoissonLikelihood; }
0162
0163    bool IsWeighted() const { return (fWeight != 0); }
0164
0165    // Use the weights in evaluating the likelihood
0166    void UseSumOfWeights() {
0167       if (fWeight == 0) return; // do nothing if it was not weighted
0168       fWeight = 1;
0169    }
0170
0171    // Use sum of the weight squared in evaluating the likelihood
0172    // (this is needed for calculating the errors)
0173    void UseSumOfWeightSquare(bool on = true) {
0174       if (fWeight == 0) return; // do nothing if it was not weighted
0175       if (on) fWeight = 2;
0176       else fWeight = 1;
0177    }
0178
0179
0180 protected:
0181
0182
0183 private:
0184
0185    /**
0186       Evaluation of the  function (required by interface)
0187     */
0188    virtual double DoEval (const double * x) const {
0189       this->UpdateNCalls();
0190       return FitUtil::Evaluate<T>::EvalPoissonLogL(BaseFCN::ModelFunction(), BaseFCN::Data(), x, fWeight, fIsExtended,
0191                                                    fNEffPoints, fExecutionPolicy);
0192    }
0193
0194    // for derivatives
0195    virtual double  DoDerivative(const double * x, unsigned int icoord ) const {
0196       Gradient(x, &fGrad[0]);
0197       return fGrad[icoord];
0198    }
0199
0200
0201       //data member
0202
0203    bool fIsExtended; ///< flag to indicate if is extended (when false is a Multinomial likelihood), default is true
0204    int fWeight;  ///< flag to indicate if needs to evaluate using weight or weight squared (default weight = 0)
0205
0206    mutable unsigned int fNEffPoints;  ///< number of effective points used in the fit
0207
0208    mutable std::vector<double> fGrad; ///< for derivatives
0209
0210    ::ROOT::EExecutionPolicy fExecutionPolicy; ///< Execution policy
0211 };
0212
0213       // define useful typedef's
0214       typedef PoissonLikelihoodFCN<ROOT::Math::IMultiGenFunction, ROOT::Math::IParamMultiFunction> PoissonLLFunction;
0215       typedef PoissonLikelihoodFCN<ROOT::Math::IMultiGradFunction, ROOT::Math::IParamMultiFunction> PoissonLLGradFunction;
0216
0217
0218    } // end namespace Fit
0219
0220 } // end namespace ROOT
0221
0222
0223 #endif /* ROOT_Fit_PoissonLikelihoodFCN */