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    double operator()(std::vector<double> const &v) const override { return fFunc.operator()(&v[0]); }
0052    double operator()(const double *v) const { return fFunc.operator()(v); }
0053    double Up() const override { return fUp; }
0054
0055    void SetErrorDef(double up) override { fUp = up; }
0056
0057
0058    /**
0059        evaluate gradient hessian and function value needed by Fumili
0060      */
0061    void EvaluateAll(std::vector<double> const &v) override;
0062
0063 private:
0064    const Function &fFunc;
0065    double fUp;
0066 };
0067
0068 template <class Function>
0069 void FumiliFCNAdapter<Function>::EvaluateAll(std::vector<double> const &v)
0070 {
0071    MnPrint print("FumiliFCNAdapter");
0072
0073    // evaluate all elements
0074    unsigned int npar = Dimension();
0075    if (npar != v.size())
0076       print.Error("npar", npar, "v.size()", v.size());
0077    assert(npar == v.size());
0078    // must distinguish case of likelihood or LS
0079
0080    std::vector<double> &grad = Gradient();
0081    std::vector<double> &hess = Hessian();
0082    // reset
0083    assert(grad.size() == npar);
0084    grad.assign(npar, 0.0);
0085    hess.assign(hess.size(), 0.0);
0086
0087    unsigned int ndata = fFunc.NPoints();
0088
0089    std::vector<double> gf(npar);
0090    std::vector<double> h(hess.size());
0091
0092    // loop on the data points
0093
0094    // if FCN is of type least-square
0095    if (fFunc.Type() == Function::kLeastSquare) {
0096       print.Debug("Chi2 FCN: Evaluate gradient and Hessian");
0097
0098       for (unsigned int i = 0; i < ndata; ++i) {
0099          // calculate data element and gradient (no need to compute Hessian)
0100          double fval = fFunc.DataElement(&v.front(), i, &gf[0]);
0101
0102          for (unsigned int j = 0; j < npar; ++j) {
0103             grad[j] += 2. * fval * gf[j];
0104             for (unsigned int k = j; k < npar; ++k) {
0105                int idx = j + k * (k + 1) / 2;
0106                hess[idx] += 2.0 * gf[j] * gf[k];
0107             }
0108          }
0109       }
0110    } else if (fFunc.Type() == Function::kLogLikelihood) {
0111       print.Debug("LogLikelihood FCN: Evaluate gradient and Hessian");
0112       for (unsigned int i = 0; i < ndata; ++i) {
0113
0114          // calculate data element and gradient: returns derivative of log(pdf)
0115          fFunc.DataElement(&v.front(), i, &gf[0]);
0116
0117          for (unsigned int j = 0; j < npar; ++j) {
0118             double gfj = gf[j];
0119             grad[j] -= gfj; // need a minus sign since is a NLL
0120             for (unsigned int k = j; k < npar; ++k) {
0121                int idx = j + k * (k + 1) / 2;
0122                hess[idx] += gfj * gf[k];
0123             }
0124          }
0125       }
0126    } else if (fFunc.Type() == Function::kPoissonLikelihood) {
0127       print.Debug("Poisson Likelihood FCN: Evaluate gradient and Hessian");
0128       // for Poisson need Hessian computed in DataElement since one needs the bin expected value ad bin observed value
0129       for (unsigned int i = 0; i < ndata; ++i) {
0130          // calculate data element and gradient
0131          fFunc.DataElement(&v.front(), i, gf.data(), h.data());
0132          for (size_t j = 0; j < npar; ++j) {
0133             grad[j] += gf[j];
0134             for (unsigned int k = j; k < npar; ++k) {
0135                int idx = j + k * (k + 1) / 2;
0136                hess[idx] += h[idx];
0137             }
0138          }
0139       }
0140    } else {
0141       print.Error("Type of fit method is not supported, it must be chi2 or log-likelihood or Poisson Likelihood");
0142    }
0143 }
0144
0145 } // end namespace Minuit2
0146
0147 } // end namespace ROOT
0148
0149 #endif // ROOT_Minuit2_FCNAdapter