LOrPE1DSymbetaKernel.hh

Go to the documentation of this file.

#ifndef NPSTAT_LORPE1DSYMBETAKERNEL_HH_
#define NPSTAT_LORPE1DSYMBETAKERNEL_HH_
 
/*!
// \file LOrPE1DSymbetaKernel.hh
//
// \brief Gaussian or symmetric beta kernels for unbinned LOrPE
//
// Author: I. Volobouev
//
// May 2022
*/
 
#include <map>
#include <vector>
#include <utility>
#include <memory>
 
#include "npstat/nm/AbsClassicalOrthoPoly1D.hh"
#include "npstat/nm/SimpleFunctors.hh"
#include "npstat/nm/fcnOrConst.hh"
 
#include "npstat/stat/BoundaryHandling.hh"
#include "npstat/stat/AbsDistribution1D.hh"
 
namespace npstat {
    class LOrPE1DSymbetaKernel
    {
    public:
        LOrPE1DSymbetaKernel(int symbetaPower, double filterDegree,
                             double leftBoundary, double rightBoundary,
                             const BoundaryHandling& bh,
                             bool memoizeBoundaryPolys = false);
 
        LOrPE1DSymbetaKernel(const LOrPE1DSymbetaKernel& r);
        LOrPE1DSymbetaKernel& operator=(const LOrPE1DSymbetaKernel& r);
 
        inline ~LOrPE1DSymbetaKernel() {deallocateThings();}
 
        inline double filterDegree() const {return filterDegree_;}
        inline double leftBoundary() const {return leftBoundary_;}
        inline double rightBoundary() const {return rightBoundary_;}
 
        inline bool isMemoizing() const {return memoizingBoundaryPolys_;}
        inline void memoize(const bool b) {memoizingBoundaryPolys_ = b;}
        inline void clearMemoizedInfo() {bpMap_.clear();}
 
        /**
        // Check if we need a boundary correction for this
        // particular combination of x and bandwidth
        */
        bool isBoundaryCorrectionNeeded(double x, double bandwidth) const;
 
        /**
        // Build boundary corrected polynomial system for this
        // combination of x and bandwidth. The code will throw
        // an exception if a boundary correction is not needed.
        */
        inline std::shared_ptr<AbsClassicalOrthoPoly1D> buildBoundaryPoly(
            const double x, const double bandwidth) const
        {
            return std::shared_ptr<AbsClassicalOrthoPoly1D>(
                buildBoundaryPolyOnTheHeap(x, bandwidth));
        } 
 
        /**
        // Normalization factor should typically be close to
        // 1/(sum of point weights).
        */
        void setNormFactor(double normfactor);
        inline double normFactor() const {return norm_;}
 
        /**
        // "Numeric" can be a pair. In this case the first element
        // of the pair must be the coordinate and the second element
        // must be the weight. The weights must be non-negative.
        */
        template <typename Numeric>
        double lorpe(double x, double bandwidth,
                     const Numeric* coords, unsigned long nCoords,
                     bool coordinatesSorted = false) const;
 
        /**
        // This function is not applying any weight. This function will
        // throw an exception if no calls to "lorpe" method have been made
        // of if a previous call to "lorpe" threw an exception.
        */
        double getLastKernelAt0() const;
 
        /**
        // Did we use boundary correction during the last "lorpe" call?
        // This function will throw an exception if no calls to "lorpe"
        // method have been made of if a previous call to "lorpe" threw
        // an exception.
        */
        bool usedBoundaryCorrection() const;
 
        bool operator==(const LOrPE1DSymbetaKernel& r) const;
        inline bool operator!=(const LOrPE1DSymbetaKernel& r) const
            {return !(*this == r);}
 
    private:
        struct BoundaryPolyInfo
        {
            std::shared_ptr<AbsClassicalOrthoPoly1D> boundaryPoly;
            std::vector<long double> polyAt0;
        };
        typedef std::map<std::pair<double,double>, BoundaryPolyInfo> BoundaryPolyMap;
 
        void deallocateThings();
        void copyInternals(const LOrPE1DSymbetaKernel& r);
        double kernelValue(const AbsClassicalOrthoPoly1D& poly,
                           const long double* poly0, double x_i) const;
        AbsClassicalOrthoPoly1D* buildBoundaryPolyOnTheHeap(
            double x, double bandwidth) const;
 
        int power_;
        unsigned maxdeg_;
        double filterDegree_;
        double leftBoundary_;
        double rightBoundary_;
        double lastShrinkage_;
        double norm_;
        double kernelRange_;
        mutable double lastCenterValue_;
        BoundaryHandling bh_;
        AbsClassicalOrthoPoly1D* poly_;
        std::vector<long double> polyAt0_;
        bool memoizingBoundaryPolys_;
 
        mutable std::vector<long double> polyBuf_;
        mutable bool usingBoundaryCorr_;
        mutable BoundaryPolyMap bpMap_;
 
#ifdef SWIG
    public:
        inline double lorpe_2(double x, double bandwidth,
                              const double* coordsd, unsigned long nCoords,
                              bool coordsSorted = false) const
        {
            return this->lorpe(x, bandwidth, coordsd, nCoords, coordsSorted);
        }
 
        inline double lorpe_2(double x, double bandwidth,
                              const std::vector<std::pair<double, double> >& vec,
                              bool coordsSorted = false) const
        {
            const unsigned long n = vec.size();
            const std::pair<double, double>* points = 0;
            if (n)
                points = &vec[0];
            return this->lorpe(x, bandwidth, points, n, coordsSorted);
        }
 
        inline AbsClassicalOrthoPoly1D* buildBoundaryPoly_2(
            double x, double bandwidth) const
        {
            return buildBoundaryPolyOnTheHeap(x, bandwidth);
        }
#endif
    };
 
    /**
    // A lightweight functor for getting LOrPE values. It will not copy
    // the kernel and will not own the data. It is the responsibility
    // of the user of this class to ensure that the lifetimes of the
    // kernel and the data exceed the lifetime of the functor.
    //
    // Intended usage is via the "LOrPE1DFunctor" utility function.
    */
    template <typename Numeric, class Fcn>
    class LOrPE1DFunctorHelper : public Functor1<double, double>
    {
    public:
        inline LOrPE1DFunctorHelper(const LOrPE1DSymbetaKernel& kernel,
                                    Fcn bandwidthFunctor,
                                    const Numeric* coords,
                                    const unsigned long nCoords,
                                    const bool coordinatesSorted)
            : kernel_(kernel), bandwidthFunctor_(bandwidthFunctor),
              coords_(coords), nCoords_(nCoords),
              coordinatesSorted_(coordinatesSorted) {}
 
        inline virtual ~LOrPE1DFunctorHelper() {}
 
        inline double operator()(const double& x) const
        {
            const double h = fcnOrConst(bandwidthFunctor_, x);
            return kernel_.lorpe(x, h, coords_, nCoords_, coordinatesSorted_);
        }
 
    private:
        const LOrPE1DSymbetaKernel& kernel_;
        Fcn bandwidthFunctor_;
        const Numeric* coords_;
        unsigned long nCoords_;
        bool coordinatesSorted_;
    };
 
    /** A convenience function for making lightweight LOrPE functors */
    template <typename Numeric, class Fcn>
    inline LOrPE1DFunctorHelper<Numeric, Fcn> LOrPE1DFunctor(
        const LOrPE1DSymbetaKernel& kernel, Fcn bandwidthFunctor,
        const Numeric* coords, const unsigned long nCoords,
        const bool coordinatesSorted = false)
    {
        return LOrPE1DFunctorHelper<Numeric,Fcn>(
            kernel, bandwidthFunctor, coords, nCoords, coordinatesSorted);
    }
 
    /** A functor for calculating LOrPE ISE */
    template <typename Numeric, class Fcn>
    class LOrPE1DKernelISEFunctor : public Functor1<double, double>
    {
    public:
        inline LOrPE1DKernelISEFunctor(
            const AbsDistribution1D& distro,
            const LOrPE1DFunctorHelper<Numeric,Fcn>& fcn)
            : distro_(distro), densityFcn_(fcn) {}
 
        inline virtual ~LOrPE1DKernelISEFunctor() {}
 
        inline double operator()(const double& x) const
        {
            const double delta = densityFcn_(x) - distro_.density(x);
            return delta*delta;
        }
 
    private:
        const AbsDistribution1D& distro_;
        const LOrPE1DFunctorHelper<Numeric,Fcn> densityFcn_;
    };
}
 
#include "npstat/stat/LOrPE1DSymbetaKernel.icc"
 
#endif // NPSTAT_LORPE1DSYMBETAKERNEL_HH_