root/Minuit2/FCNGradAdapter.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_FCNGradAdapter
0011 #define ROOT_Minuit2_FCNGradAdapter
0012
0013 #include "Minuit2/FCNGradientBase.h"
0014 #include "Minuit2/MnPrint.h"
0015
0016 #include <vector>
0017 #include <functional>
0018
0019 namespace ROOT {
0020
0021 namespace Minuit2 {
0022
0023 /**
0024
0025
0026 template wrapped class for adapting to FCNBase signature a IGradFunction
0027
0028 @author Lorenzo Moneta
0029
0030 @ingroup Minuit
0031
0032 */
0033
0034 template <class Function>
0035 class FCNGradAdapter : public FCNGradientBase {
0036
0037 public:
0038    FCNGradAdapter(const Function &f, double up = 1.) : fFunc(f), fUp(up), fGrad(std::vector<double>(fFunc.NDim())) {}
0039
0040    ~FCNGradAdapter() override {}
0041
0042    double operator()(const std::vector<double> &v) const override { return fFunc.operator()(&v[0]); }
0043    double operator()(const double *v) const { return fFunc.operator()(v); }
0044
0045    double Up() const override { return fUp; }
0046
0047    std::vector<double> Gradient(const std::vector<double> &v) const override
0048    {
0049       fFunc.Gradient(&v[0], &fGrad[0]);
0050       return fGrad;
0051    }
0052    std::vector<double> GradientWithPrevResult(const std::vector<double> &v, double *previous_grad, double *previous_g2,
0053                                               double *previous_gstep) const override
0054    {
0055       fFunc.GradientWithPrevResult(&v[0], &fGrad[0], previous_grad, previous_g2, previous_gstep);
0056       return fGrad;
0057    }
0058    // forward interface
0059    // virtual double operator()(int npar, double* params,int iflag = 4) const;
0060    bool CheckGradient() const override { return false; }
0061
0062    GradientParameterSpace gradParameterSpace() const override {
0063       if (fFunc.returnsInMinuit2ParameterSpace()) {
0064          return GradientParameterSpace::Internal;
0065       } else {
0066          return GradientParameterSpace::External;
0067       }
0068    }
0069
0070    /// return second derivatives (diagonal of the Hessian matrix)
0071    std::vector<double> G2(const std::vector<double> & x) const override {
0072       if (fG2Func)
0073          return fG2Func(x);
0074       if (fHessianFunc) {
0075          unsigned int n = fFunc.NDim();
0076          if (fG2Vec.empty() ) fG2Vec.resize(n);
0077          if (fHessian.empty() ) fHessian.resize(n*n);
0078          fHessianFunc(x,fHessian.data());
0079          if (!fHessian.empty()) {
0080             // get diagonal element of h
0081             for (unsigned int i = 0; i < n; i++)
0082                fG2Vec[i] = fHessian[i*n+i];
0083          }
0084          else fG2Vec.clear();
0085       }
0086       else
0087          if (!fG2Vec.empty()) fG2Vec.clear();
0088       return fG2Vec;
0089    }
0090
0091    /// compute Hessian. Return Hessian as a std::vector of size(n*n)
0092    std::vector<double> Hessian(const std::vector<double> & x ) const override {
0093       unsigned int n = fFunc.NDim();
0094       if (fHessianFunc) {
0095          if (fHessian.empty() ) fHessian.resize(n * n);
0096          bool ret = fHessianFunc(x,fHessian.data());
0097          if (!ret) {
0098             fHessian.clear();
0099             fHessianFunc = nullptr;
0100          }
0101       } else {
0102          fHessian.clear();
0103       }
0104
0105       return fHessian;
0106    }
0107
0108    bool HasG2() const override {
0109       return bool(fG2Func);
0110    }
0111    bool HasHessian() const override {
0112       return bool(fHessianFunc);
0113    }
0114
0115    template<class Func>
0116    void SetG2Function(Func f) { fG2Func = f;}
0117
0118    template<class Func>
0119    void SetHessianFunction(Func f) { fHessianFunc = f;}
0120
0121    void SetErrorDef(double up) override { fUp = up; }
0122
0123 private:
0124    const Function &fFunc;
0125    double fUp;
0126    mutable std::vector<double> fGrad;
0127    mutable std::vector<double> fHessian;
0128    mutable std::vector<double> fG2Vec;
0129
0130    std::function<std::vector<double>(const std::vector<double> &)> fG2Func;
0131    mutable std::function<bool(const std::vector<double> &, double *)> fHessianFunc;
0132 };
0133
0134 } // end namespace Minuit2
0135
0136 } // end namespace ROOT
0137
0138 #endif // ROOT_Minuit2_FCNGradAdapter