root/Minuit2/FumiliFCNAdapter.h

0001 // @(#)root/minuit2:$Id$
0002 // Author: L. Moneta    10/2006
0003
0004 /**********************************************************************
0005  *                                                                    *
0006  * Copyright (c) 2006 ROOT Foundation,  CERN/PH-SFT                   *
0007  *                                                                    *
0008  **********************************************************************/
0009
0010 #ifndef ROOT_Minuit2_FumiliFCNAdapter
0011 #define ROOT_Minuit2_FumiliFCNAdapter
0012
0013 #include "Minuit2/FumiliFCNBase.h"
0014
0015 #include "Math/FitMethodFunction.h"
0016
0017 #include "Minuit2/MnPrint.h"
0018
0019 // #ifndef ROOT_Math_Util
0020 // #include "Math/Util.h"
0021 // #endif
0022
0023 #include <cmath>
0024 #include <cassert>
0025 #include <vector>
0026
0027 namespace ROOT {
0028
0029 namespace Minuit2 {
0030
0031 /**
0032
0033
0034 template wrapped class for adapting to FumiliFCNBase signature
0035
0036 @author Lorenzo Moneta
0037
0038 @ingroup Minuit
0039
0040 */
0041
0042 template <class Function>
0043 class FumiliFCNAdapter : public FumiliFCNBase {
0044
0045 public:
0046    //   typedef ROOT::Math::FitMethodFunction Function;
0047    typedef typename Function::Type_t Type_t;
0048
0049    FumiliFCNAdapter(const Function &f, unsigned int ndim, double up = 1.) : FumiliFCNBase(ndim), fFunc(f), fUp(up) {}
0050
0051    ~FumiliFCNAdapter() override {}
0052
0053    double operator()(const std::vector<double> &v) const override { return fFunc.operator()(&v[0]); }
0054    double operator()(const double *v) const { return fFunc.operator()(v); }
0055    double Up() const override { return fUp; }
0056
0057    void SetErrorDef(double up) override { fUp = up; }
0058
0059    // virtual std::vector<double> Gradient(const std::vector<double>&) const;
0060
0061    // forward interface
0062    // virtual double operator()(int npar, double* params,int iflag = 4) const;
0063
0064    /**
0065        evaluate gradient hessian and function value needed by fumili
0066      */
0067    void EvaluateAll(const std::vector<double> &v) override;
0068
0069 private:
0070    // data member
0071
0072    const Function &fFunc;
0073    double fUp;
0074 };
0075
0076 template <class Function>
0077 void FumiliFCNAdapter<Function>::EvaluateAll(const std::vector<double> &v)
0078 {
0079    MnPrint print("FumiliFCNAdapter");
0080
0081    // typedef FumiliFCNAdapter::Function Function;
0082
0083    // evaluate all elements
0084    unsigned int npar = Dimension();
0085    if (npar != v.size())
0086       print.Error("npar", npar, "v.size()", v.size());
0087    assert(npar == v.size());
0088    // must distinguish case of likelihood or LS
0089
0090    std::vector<double> &grad = Gradient();
0091    std::vector<double> &hess = Hessian();
0092    // reset
0093    assert(grad.size() == npar);
0094    grad.assign(npar, 0.0);
0095    hess.assign(hess.size(), 0.0);
0096
0097    unsigned int ndata = fFunc.NPoints();
0098
0099    std::vector<double> gf(npar);
0100    std::vector<double> h(hess.size());
0101
0102    // loop on the data points
0103
0104    // if FCN is of type least-square
0105    if (fFunc.Type() == Function::kLeastSquare) {
0106       print.Debug("Chi2 FCN: Evaluate gradient and Hessian");
0107
0108       for (unsigned int i = 0; i < ndata; ++i) {
0109          // calculate data element and gradient (no need to compute Hessian)
0110          double fval = fFunc.DataElement(&v.front(), i, &gf[0]);
0111
0112          for (unsigned int j = 0; j < npar; ++j) {
0113             grad[j] += 2. * fval * gf[j];
0114             for (unsigned int k = j; k < npar; ++k) {
0115                int idx = j + k * (k + 1) / 2;
0116                hess[idx] += 2.0 * gf[j] * gf[k];
0117             }
0118          }
0119       }
0120    } else if (fFunc.Type() == Function::kLogLikelihood) {
0121       print.Debug("LogLikelihood FCN: Evaluate gradient and Hessian");
0122       for (unsigned int i = 0; i < ndata; ++i) {
0123
0124          // calculate data element and gradient: returns derivative of log(pdf)
0125          fFunc.DataElement(&v.front(), i, &gf[0]);
0126
0127          for (unsigned int j = 0; j < npar; ++j) {
0128             double gfj = gf[j];
0129             grad[j] -= gfj; // need a minus sign since is a NLL
0130             for (unsigned int k = j; k < npar; ++k) {
0131                int idx = j + k * (k + 1) / 2;
0132                hess[idx] += gfj * gf[k];
0133             }
0134          }
0135       }
0136    } else if (fFunc.Type() == Function::kPoissonLikelihood) {
0137       print.Debug("Poisson Likelihood FCN: Evaluate gradient and Hessian");
0138       // for Poisson need Hessian computed in DataElement since one needs the bin expected value ad bin observed value
0139      for (unsigned int i = 0; i < ndata; ++i) {
0140          // calculate data element and gradient
0141          fFunc.DataElement(&v.front(), i, gf.data(), h.data());
0142          for (size_t j = 0; j < npar; ++j) {
0143             grad[j] += gf[j];
0144             for (unsigned int k = j; k < npar; ++k) {
0145                int idx = j + k * (k + 1) / 2;
0146                hess[idx] += h[idx];
0147             }
0148          }
0149      }
0150    } else {
0151       print.Error("Type of fit method is not supported, it must be chi2 or log-likelihood or Poisson Likelihood");
0152    }
0153 }
0154
0155 } // end namespace Minuit2
0156
0157 } // end namespace ROOT
0158
0159 #endif // ROOT_Minuit2_FCNAdapter