MinuitPLogliOSDE1DFcn.hh

Go to the documentation of this file.

#ifndef NPSI_MINUITPLOGLIOSDE1DFCN_HH_
#define NPSI_MINUITPLOGLIOSDE1DFCN_HH_
 
/*!
// \file MinuitPLogliOSDE1DFcn.hh
//
// \brief Minuit function to minimize for pseudo-log-likelihood OSDE
//
// Author: I. Volobouev
//
// March 2023
*/
 
#include <cassert>
 
#include "Minuit2/FCNGradientBase.h"
#include "npstat/stat/PLogliOSDE1D.hh"
 
namespace npsi {
    /**
    // Target minimization function adapter class for running PLogli OSDE
    // fits by Minuit2. Note that PLogliOSDE1D object is not copied,
    // only a reference is held, with the corresponding consequences
    // for the object lifetimes.
    */
    class MinuitPLogliOSDE1DFcn : public ROOT::Minuit2::FCNGradientBase
    {
    public:
        inline MinuitPLogliOSDE1DFcn(
            const npstat::PLogliOSDE1D& i_osde,
            const unsigned nCoeffsExpected,
            const double up=0.05)
            : osde_(i_osde),
              lastGradient_(nCoeffsExpected),
              lastResult_(0.0),
              up_(up),
              nCoeffs_(nCoeffsExpected),
              first_(true)
        {
            assert(nCoeffs_);
        }
 
        inline virtual ~MinuitPLogliOSDE1DFcn() {}
 
        inline const npstat::PLogliOSDE1D& getOsde() const {return osde_;}
 
        // Make Minuit look for a better minimum by reducing
        // the recommended "Up" value by an order of magnitude
        inline double Up() const {return up_;}
 
        // Do not waste time checking gradient calculations. They are correct.
        inline bool CheckGradient() const {return false;}
 
        inline virtual double operator()(const std::vector<double>& x) const
        {
            calculateMemoized(x);
            return lastResult_;
        }
 
        inline std::vector<double> Gradient(const std::vector<double>& x) const
        {
            // Returning std::vector<double> on the stack is
            // a pretty inefficient thing to do. However,
            // this is how Minuit2 interface is designed.
            calculateMemoized(x);
            return lastGradient_;
        }
 
    private:
        MinuitPLogliOSDE1DFcn();
 
        // Minuit2 has separate interfaces for calculating the function
        // value and the gradient. I have yet to see a non-trivial realistic
        // problem in which the function value can not be easily obtained
        // along with the gradient. It would be much better to have them
        // calculated together. We will alleviate this design problem
        // by memoizing the argument / result pairs so that we do not
        // have to repeat the calculation when the gradient is wanted
        // at the same point where the function was just evaluated.
        inline void calculateMemoized(const std::vector<double>& x) const
        {
            if (x != lastX_ || first_)
            {
                assert(x.size() == nCoeffs_);
                lastResult_ = -osde_.pseudoLogli(
                    &x[0], nCoeffs_, &lastGradient_[0]);
                lastX_ = x;
                for (unsigned i=0; i<nCoeffs_; ++i)
                    lastGradient_[i] *= -1.0;
                first_ = false;
            }
        }
 
        const npstat::PLogliOSDE1D& osde_;
        mutable std::vector<double> lastX_;
        mutable std::vector<double> lastGradient_;
        mutable double lastResult_;
        double up_;
        unsigned nCoeffs_;
        mutable bool first_;
    };
}
 
#endif // MINUITPLOGLIOSDE1DFCN_HH_