MinuitNeymanOSDE1DFcn.hh

Go to the documentation of this file.

#ifndef NPSI_MINUITNEYMANOSDE1DFCN_HH_
#define NPSI_MINUITNEYMANOSDE1DFCN_HH_
 
/*!
// \file MinuitNeymanOSDE1DFcn.hh
//
// \brief Minuit function to minimize for Neyman OSDE
//
// Author: I. Volobouev
//
// March 2023
*/
 
#include <cassert>
 
#include "Minuit2/FCNGradientBase.h"
#include "npstat/stat/NeymanOSDE1D.hh"
 
namespace npsi {
    /**
    // Target minimization function adapter class for running Neyman OSDE
    // fits by Minuit2. Note that NeymanOSDE1D object is not copied,
    // only a reference is held, with the corresponding consequences
    // for the object lifetimes.
    */
    class MinuitNeymanOSDE1DFcn : public ROOT::Minuit2::FCNGradientBase
    {
    public:
        inline MinuitNeymanOSDE1DFcn(
            const npstat::NeymanOSDE1D& i_osde,
            const unsigned nCoeffsExpected,
            const double* shrinkages, const unsigned nShrinkages,
            const double up=0.1)
            : osde_(i_osde),
              shrinkages_(shrinkages, shrinkages+nShrinkages),
              lastGradient_(nCoeffsExpected),
              polyStatistics_(nShrinkages),
              lastResult_(-1.0),
              up_(up),
              nCoeffs_(nCoeffsExpected)
        {
            assert(nCoeffs_);
            assert(nShrinkages);
            assert(shrinkages);
        }
 
        inline virtual ~MinuitNeymanOSDE1DFcn() {}
 
        inline const npstat::NeymanOSDE1D& 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:
        MinuitNeymanOSDE1DFcn();
 
        // 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_ || lastResult_ < 0.0)
            {
                assert(x.size() == nCoeffs_);
                lastResult_ = osde_.chiSquare(
                    &x[0], nCoeffs_,
                    &shrinkages_[0], shrinkages_.size(),
                    &polyStatistics_[0], &lastGradient_[0]);
                lastX_ = x;
            }
        }
 
        const npstat::NeymanOSDE1D& osde_;
        std::vector<double> shrinkages_;
        mutable std::vector<double> lastX_;
        mutable std::vector<double> lastGradient_;
        mutable std::vector<double> polyStatistics_;
        mutable double lastResult_;
        double up_;
        unsigned nCoeffs_;
    };
}
 
#endif // MINUITNEYMANOSDE1DFCN_HH_