LOrPE1DCV.hh

Go to the documentation of this file.

#ifndef NPSTAT_LORPE1DCV_HH_
#define NPSTAT_LORPE1DCV_HH_
 
/*!
// \file LOrPE1DCV.hh
//
// \brief Unbinned density estimation by LOrPE with cross-validation
//
// Author: I. Volobouev
//
// June 2022
*/
 
#include <map>
#include <cmath>
#include <vector>
#include <climits>
#include <cassert>
#include <memory>
#include <utility>
#include <stdexcept>
#include <algorithm>
#include <sstream>
 
#include "npstat/stat/LOrPE1D.hh"
 
namespace npstat {
    template <class Helper>
    class LOrPE1DCVDensityFunctor : public Functor1<double, double>
    {
    public:
        inline LOrPE1DCVDensityFunctor(const Helper& helper,
                                       const double filterDegree,
                                       const double bwFactor)
            : helper_(helper),
              filterDegree_(filterDegree),
              bwFactor_(bwFactor)
        {
            if (filterDegree < 0.0) throw std::invalid_argument(
                "In npstat::LOrPE1DCVDensityFunctor constructor: "
                "filter degree can not be negative");
            if (bwFactor <= 0.0) throw std::invalid_argument(
                "In npstat::LOrPE1DCVDensityFunctor constructor: "
                "bandwidth factor must be positive");
        }
 
        inline virtual ~LOrPE1DCVDensityFunctor() {}
 
        inline double operator()(const double& x) const
            {return helper_.density(x, filterDegree_, bwFactor_);}
 
    private:
        const Helper& helper_;
        double filterDegree_;
        double bwFactor_;
    };
 
    template <class Helper>
    class LOrPE1DCVLocalizingWeightFunctor : public Functor1<double, double>
    {
    public:
        inline LOrPE1DCVLocalizingWeightFunctor(const Helper& helper)
            : helper_(helper) {}
 
        inline virtual ~LOrPE1DCVLocalizingWeightFunctor() {}
 
        inline double operator()(const double& x) const
            {return helper_.cvLocalizingWeight(x);}
 
    private:
        const Helper& helper_;
    };
 
 
    template <typename Numeric, class BwFcn, class WFcn>
    class LOrPE1DCVFunctorHelper : public Functor1<double, double>,
                                   public Functor2<double, double, double>
    {
    public:
        typedef Numeric point_type;
 
        inline LOrPE1DCVFunctorHelper(
            const int i_symbetaPower,
            const double i_leftBoundary,
            const double i_rightBoundary,
            const BoundaryHandling& i_bh,
            const BwFcn& i_bandwidthFunctor,
            const WFcn& i_localizingWeight,
            const double i_localizingWeightXmin,
            const double i_localizingWeightXmax,
            const unsigned i_nIntegIntervals,
            const unsigned i_nIntegPoints,
            const bool i_normalizeLOrPEEstimate,
            const std::vector<Numeric>& i_coords)
            : coords_(i_coords),
              kernel_(0),
              lastBwFactor_(-1.0),
              lastNorm_(-1.0),
              boundFilterDegree_(-1.0),
              symbetaPower_(i_symbetaPower),
              leftBoundary_(i_leftBoundary),
              rightBoundary_(i_rightBoundary),
              bh_(i_bh),
              bandwidthFunctor_(i_bandwidthFunctor),
              localizingWeight_(i_localizingWeight),
              xminForWeight_(i_localizingWeightXmin),
              xmaxForWeight_(i_localizingWeightXmax),
              nIntegIntervals_(i_nIntegIntervals),
              nIntegPoints_(i_nIntegPoints),
              normalize_(i_normalizeLOrPEEstimate),
              memoizeKernels_(false),
              memoizeNorms_(false),
              useMemoisingKernels_(false)
        {
            if (coords_.empty()) throw std::invalid_argument(
                "In npstat::LOrPE1DCVFunctorHelper constructor: empty point sample");
            if (leftBoundary_ >= rightBoundary_) throw std::invalid_argument(
                "In npstat::LOrPE1DCVFunctorHelper constructor: invalid interval boundaries");
            if (xminForWeight_ >= xmaxForWeight_) throw std::invalid_argument(
                "In npstat::LOrPE1DCVFunctorHelper constructor: invalid weight boundaries");
            if (!nIntegIntervals_) throw std::invalid_argument(
                "In npstat::LOrPE1DCVFunctorHelper constructor: "
                "number of integration intervals must be positive");
            if (!GaussLegendreQuadrature::isAllowed(nIntegPoints_)) throw std::invalid_argument(
                "In npstat::LOrPE1DCVFunctorHelper constructor: "
                "unsupported number of integration points");
            std::sort(coords_.begin(), coords_.end());
            calculateEffectiveSampleSize();
        }
 
        inline virtual ~LOrPE1DCVFunctorHelper() {}
 
        inline void bindFilterDegree(const double deg)
        {
            if (deg < 0.0) throw std::invalid_argument(
                "In npstat::LOrPE1DCVFunctorHelper::bindFilterDegree: "
                "filter degree can not be negative");
            boundFilterDegree_ = deg;
        }
        inline void unbindFilterDegree()
            {boundFilterDegree_ = -1.0;}
        inline double boundFilterDegree() const
            {return boundFilterDegree_;}
 
        /**
        // This method with throw std::runtime_error if Numeric
        // is not an std::pair
        */
        inline virtual void setWeights(const double* weights, unsigned long nWeights)
        {
            if (nWeights != coords_.size()) throw std::invalid_argument(
                "In npstat::LOrPE1DCVFunctorHelper::setWeights: "
                "incompatible length of the weight array");
            Private::SetPairSeconds<Numeric>::set(&coords_[0], weights, nWeights);
            calculateEffectiveSampleSize();
            mapOfNorms_.clear();
            lastBwFactor_ = -1.0;
        }
 
        inline int symbetaPower() const {return symbetaPower_;}
        inline double leftBoundary() const {return leftBoundary_;}
        inline double rightBoundary() const {return rightBoundary_;}
        inline BoundaryHandling bh() const {return bh_;}
        inline double localizingWeightXmin() const {return xminForWeight_;}
        inline double localizingWeighXmax() const {return xmaxForWeight_;}
        inline unsigned nIntegIntervals() const {return nIntegIntervals_;}
        inline unsigned nIntegPoints() const {return nIntegPoints_;}
 
        /** Is the density estimate auto-normalizing? */
        inline bool normalizingLOrPEEstimate() const {return normalize_;}
 
        inline const std::vector<Numeric>& coords() const {return coords_;}
        inline unsigned long nCoords() const {return coords_.size();}
        inline Numeric minCoordinate() const {return coords_[0];}
        inline Numeric maxCoordinate() const {return coords_.back();}
        inline double totalSampleWeight() const {return sampleWeight_;}
        inline double effectiveSampleSize() const {return effSampleSize_;}
 
        inline bool isMemoizingKernels() const {return memoizeKernels_;}
        inline bool isMemoizingNorms() const {return memoizeNorms_;}
        inline bool isUsingMemoisingKernels() const {return useMemoisingKernels_;}
 
        // The following modifiers should be called immediately
        // after object creation
        inline void memoizeKernels(const bool b) {memoizeKernels_ = b;}
        inline void memoizeNorms(const bool b) {memoizeNorms_ = b;}
        inline void useMemoisingKernels(const bool b)
        {
            assert(!kernel_);
            useMemoisingKernels_ = b;
        }
 
        // Clear some of the memoized information. Note that
        // this will not clear the info memoized by the current
        // kernel. Call "clearKernel" function after this one
        // if you must clear everything.
        inline void clearMemoizedInfo()
        {
            mapOfNorms_.clear();
            mapOfKernels_.clear();
        }
 
        inline void clearKernel()
        {
            kernel_ = std::shared_ptr<LOrPE1DSymbetaKernel>(0);
            lastBwFactor_ = -1.0;
        }
 
        inline double pointCoordinate(const unsigned long i) const
        {
            double x, w;
            Private::coordAndWeight(coords_[i], &x, &w);
            return x;
        }
 
        inline double pointWeight(const unsigned long i) const
        {
            double x, w;
            Private::coordAndWeight(coords_[i], &x, &w);
            return w;
        }
 
        inline double bandwidthCurve(const double x) const
            {return fcnOrConst(bandwidthFunctor_, x);}
 
        inline double cvLocalizingWeight(const double x) const
        {
            const double xlolim = std::max(leftBoundary_, xminForWeight_);
            const double xuplim = std::min(rightBoundary_, xmaxForWeight_);
            if (x >= xlolim && x <= xuplim)
                return fcnOrConst(localizingWeight_, x);
            else
                return 0.0;
        }
 
        inline LOrPE1DCVLocalizingWeightFunctor<LOrPE1DCVFunctorHelper<Numeric, BwFcn, WFcn> >
        cvLocalizingWeightFunctor() const
        {
            return LOrPE1DCVLocalizingWeightFunctor<LOrPE1DCVFunctorHelper<Numeric, BwFcn, WFcn> >(*this);
        }
 
        inline double density(const double x, const double filterDegree,
                              const double bwFactor) const
        {
            setupKernel(filterDegree, bwFactor);
            const double bw = bwFactor*fcnOrConst(bandwidthFunctor_, x);
            return kernel_->lorpe(x, bw, &coords_[0], coords_.size(), true);
        }
 
        inline double density(const double x, const double bwFactor) const
        {
            validateBoundDegree("density");
            return density(x, boundFilterDegree_, bwFactor);
        }
 
        inline LOrPE1DCVDensityFunctor<LOrPE1DCVFunctorHelper<Numeric, BwFcn, WFcn> >
        densityFunctor(const double filterDegree, const double bwFactor) const
        {
            return LOrPE1DCVDensityFunctor<LOrPE1DCVFunctorHelper<Numeric, BwFcn, WFcn> >(
                *this, filterDegree, bwFactor);
        }
 
        inline LOrPE1DCVDensityFunctor<LOrPE1DCVFunctorHelper<Numeric, BwFcn, WFcn> >
        densityFunctor(const double bwFactor) const
        {
            validateBoundDegree("densityFunctor");
            return densityFunctor(boundFilterDegree_, bwFactor);
        }
 
        inline LOrPE1D<Numeric> getLOrPE1D(
            const double filterDegree, const double bwFactor) const
        {
            setupKernel(filterDegree, bwFactor);
            LOrPE1D<Numeric> lorpe(*kernel_, coords_);
            lorpe.setNormFactor(kernel_->normFactor());
            return lorpe;
        }
 
        inline LOrPE1D<Numeric> getLOrPE1D(const double bwFactor) const
        {
            validateBoundDegree("getLOrPE1D");
            return getLOrPE1D(boundFilterDegree_, bwFactor);
        }
 
        inline double operator()(
            const double& filterDegree, const double& bwFactor) const
        {
            setupKernel(filterDegree, bwFactor);
 
            const unsigned long ncoords = coords_.size();
            const Numeric* first = &coords_[0];
            const Numeric* last = first + ncoords;
            const double xlolim = std::max(leftBoundary_, xminForWeight_);
            const double xuplim = std::min(rightBoundary_, xmaxForWeight_);
            if (ncoords > 10UL)
            {
                const Numeric xlow = Private::WeightedCoord<Numeric>(xlolim, 0.0).value;
                if (xlow > *first)
                    first = std::lower_bound(first, last, xlow);
 
                const Numeric xup = Private::WeightedCoord<Numeric>(xuplim, DBL_MAX).value;
                if (xup < *(last - 1))
                    last = std::upper_bound(first, last, xup);
            }
 
            long double msum = 0.0L;
            double xi, w;
            for (const Numeric* ptr = first; ptr != last; ++ptr)
            {
                Private::coordAndWeight(*ptr, &xi, &w);
                assert(xi >= leftBoundary_ && xi <= rightBoundary_);
                assert(w >= 0.0);
                if (w > 0.0 && xi >= xlolim && xi <= xuplim)
                {
                    const double localWeight = fcnOrConst(localizingWeight_, xi);
                    assert(localWeight >= 0.0);
                    if (localWeight > 0.0)
                    {
                        const double bw = bwFactor*fcnOrConst(bandwidthFunctor_, xi);
                        const double estimate = kernel_->lorpe(
                            xi, bw, &coords_[0], ncoords, true);
                        const double selfC = kernel_->getLastKernelAt0();
                        const double loo = (estimate - w*selfC)*sampleWeight_/(sampleWeight_ - w);
                        msum += pointContribution(selfC, loo, w, localWeight);
                    }
                }
            }
            return valueToMaximize(msum);
        }
 
        inline double operator()(const double& bwFactor) const
        {
            validateBoundDegree("operator()");
            return operator()(boundFilterDegree_, bwFactor);
        }
 
        /**
        // This method will return meaningful results only
        // if the automatic normalization is not turned on
        */
        inline double densityIntegral(
            const double filterDegree, const double bwFactor) const
        {
            setupKernel(filterDegree, bwFactor);
            GaussLegendreQuadrature glq(nIntegPoints_);
            return glq.integrate(this->bwDensityFunctor(bwFactor),
                                 leftBoundary_, rightBoundary_, nIntegIntervals_);
        }
 
        inline double densityIntegral(const double bwFactor) const
        {
            validateBoundDegree("densityIntegral");
            return densityIntegral(boundFilterDegree_, bwFactor);
        }
 
        inline double integratedSquaredError(
            const AbsDistribution1D& distro,
            const double filterDegree, const double bwFactor) const
        {
            setupKernel(filterDegree, bwFactor);
            LOrPE1DKernelISEFunctor<Numeric,BwFunctor> iseFcn(
                distro, this->bwDensityFunctor(bwFactor));
            GaussLegendreQuadrature glq(nIntegPoints_);
            return glq.integrate(iseFcn, leftBoundary_,
                                 rightBoundary_, nIntegIntervals_);
        }
 
        inline double integratedSquaredError(
            const AbsDistribution1D& distro, const double bwFactor) const
        {
            validateBoundDegree("integratedSquaredError");
            return integratedSquaredError(distro, boundFilterDegree_, bwFactor);
        }
 
        /** Is this class capable of cross-validation? */
        inline static bool cvCapable() {return true;}
 
    protected:
        typedef MultFunctor<BwFcn> BwFunctor;
        typedef std::map<std::pair<double,double>,double> MapOfNorms;
        typedef std::map<double,std::shared_ptr<LOrPE1DSymbetaKernel> > MapOfKernels;
 
        virtual double pointContribution(
            double selfC, double leaveOneOutEstimate,
            double pointWeight, double localWeight) const = 0;
        virtual double valueToMaximize(double sum) const = 0;
 
        std::vector<Numeric> coords_;
        double sampleWeight_;
        double effSampleSize_;
 
        mutable std::shared_ptr<LOrPE1DSymbetaKernel> kernel_;
        mutable double lastBwFactor_;
        mutable double lastNorm_;
 
        double boundFilterDegree_;
 
        int symbetaPower_;
        double leftBoundary_;
        double rightBoundary_;
        BoundaryHandling bh_;
 
        BwFcn bandwidthFunctor_;
        WFcn localizingWeight_;
        double xminForWeight_;
        double xmaxForWeight_;
        unsigned nIntegIntervals_;
        unsigned nIntegPoints_;
        bool normalize_;
 
        bool memoizeKernels_;
        bool memoizeNorms_;
        bool useMemoisingKernels_;
 
        mutable MapOfNorms mapOfNorms_;
        mutable MapOfKernels mapOfKernels_;
 
        inline LOrPE1DFunctorHelper<Numeric,BwFunctor> bwDensityFunctor(
            const double bwFactor) const
        {
            return LOrPE1DFunctorHelper<Numeric,BwFunctor>(
                *kernel_, BwFunctor(bandwidthFunctor_, bwFactor),
                &coords_[0], coords_.size(), true);
        }
 
    private:
        inline void validateBoundDegree(const char* where) const
        {
            if (boundFilterDegree_ < 0.0)
            {
                std::ostringstream os;
                os << "In npstat::LOrPE1DCVFunctorHelper::" << where
                   << ": can't use the bandwidth factor alone, "
                   << "please bind the filter degree first";
                throw std::runtime_error(os.str());
            }
        }
 
        inline void setupKernel(const double filterDegree,
                                const double bwFactor) const
        {
            if (filterDegree < 0.0) throw std::invalid_argument(
                "In npstat::LOrPE1DCVFunctorHelper::setupKernel: "
                "filter degree can not be negative");
            if (bwFactor <= 0.0) throw std::invalid_argument(
                "In npstat::LOrPE1DCVFunctorHelper::setupKernel: "
                "bandwidth factor must be positive");
 
            const bool needAnotherKernel =
                !(kernel_ && kernel_->filterDegree() == filterDegree);
            if (needAnotherKernel)
                makeKernel(filterDegree);
 
            // Check if precise auto normalization is requested
            if (normalize_)
            {
                if (bwFactor == lastBwFactor_ && !needAnotherKernel)
                {
                    assert(lastNorm_ > 0.0);
                    kernel_->setNormFactor(lastNorm_);
                }
                else
                {
                    normalizeEstimate(bwFactor);
                    lastBwFactor_ = bwFactor;
                    lastNorm_ = kernel_->normFactor();
                }
            }
            else
            {
                // Set up reasonable default normalization
                kernel_->setNormFactor(1.0/sampleWeight_);
                lastBwFactor_ = bwFactor;
                lastNorm_ = kernel_->normFactor();
            }
        }
 
        inline void makeKernel(const double filterDegree) const
        {
            bool reallyMakeNewKernel = true;
            if (memoizeKernels_)
            {
                kernel_ = mapOfKernels_[filterDegree];
                reallyMakeNewKernel = !kernel_;
            }
            if (reallyMakeNewKernel)
            {
                LOrPE1DSymbetaKernel* k = new LOrPE1DSymbetaKernel(
                    symbetaPower_, filterDegree, leftBoundary_,
                    rightBoundary_, bh_, useMemoisingKernels_);
                kernel_ = std::shared_ptr<LOrPE1DSymbetaKernel>(k);
 
                if (memoizeKernels_)
                    mapOfKernels_[filterDegree] = kernel_;
            }
        }
 
        inline void normalizeEstimate(const double bwFactor) const
        {
            assert(kernel_);
 
            bool reallyRecalculateNorm = true;
            if (memoizeNorms_)
            {
                const double deg = kernel_->filterDegree();
                MapOfNorms::iterator it = mapOfNorms_.find(std::make_pair(deg, bwFactor));
                if (it != mapOfNorms_.end())
                {
                    reallyRecalculateNorm = false;
                    kernel_->setNormFactor(it->second);
                }
            }
 
            if (reallyRecalculateNorm)
            {
                const double defaultNorm = 1.0/sampleWeight_;
                kernel_->setNormFactor(defaultNorm);
                GaussLegendreQuadrature glq(nIntegPoints_);
                const double integ = glq.integrate(this->bwDensityFunctor(bwFactor),
                                                   leftBoundary_, rightBoundary_, nIntegIntervals_);
                assert(integ > 0.0);
                const double newNorm = defaultNorm/integ;
                kernel_->setNormFactor(newNorm);
 
                if (memoizeNorms_)
                {
                    const double deg = kernel_->filterDegree();
                    mapOfNorms_[std::make_pair(deg, bwFactor)] = newNorm;
                }
            }
        }
 
        inline void calculateEffectiveSampleSize()
        {
            sampleWeight_ = Private::SampleWeight<Numeric>(coords_).value;
            if (sampleWeight_ <= 0.0) throw std::invalid_argument(
                "In npstat::LOrPE1DCV::calculateEffectiveSampleSize: "
                "total sample weight must be positive");
            const double sumWsq = Private::SampleSquaredWeight<Numeric>(coords_).value;
            effSampleSize_ = sampleWeight_/sumWsq*sampleWeight_;
        }
 
#ifdef SWIG
    public:
        inline std::vector<double> sortedCoords() const
        {
            const unsigned long sz = coords_.size();
            std::vector<double> out;
            out.reserve(sz);
            double x, w;
            for (unsigned long i = 0; i < sz; ++i)
            {
                Private::coordAndWeight(coords_[i], &x, &w);
                out.push_back(x);
            }
            return out;
        }
 
        inline std::vector<double> sortedCoordWeights() const
        {
            const unsigned long sz = coords_.size();
            std::vector<double> out;
            out.reserve(sz);
            double x, w;
            for (unsigned long i = 0; i < sz; ++i)
            {
                Private::coordAndWeight(coords_[i], &x, &w);
                out.push_back(w);
            }
            return out;
        }
#endif
    };
 
    template <typename Numeric, class BwFcn, class WFcn>
    class LOrPE1DLSCVFunctorHelper : public LOrPE1DCVFunctorHelper<Numeric, BwFcn, WFcn>
    {
    public:
        typedef LOrPE1DCVFunctorHelper<Numeric, BwFcn, WFcn> Base;
 
        using Base::LOrPE1DCVFunctorHelper;
 
        inline virtual ~LOrPE1DLSCVFunctorHelper() {}
 
    protected:
        inline virtual double pointContribution(
            const double /* selfC */, const double leaveOneOutEstimate,
            const double pointWeight, const double localWeight) const
        {
            return pointWeight*localWeight*leaveOneOutEstimate;
        }
 
        inline virtual double valueToMaximize(const double sum) const
        {
            const double xlolim = std::max(Base::leftBoundary_, Base::xminForWeight_);
            const double xuplim = std::min(Base::rightBoundary_, Base::xmaxForWeight_);
            GaussLegendreQuadrature glq(Base::nIntegPoints_);
            const double integ = glq.integrate(
                make_DensitySquaredTimesWeight(
                    this->bwDensityFunctor(Base::lastBwFactor_), Base::localizingWeight_),
                xlolim, xuplim, Base::nIntegIntervals_);
            return -(integ - sum*2.0/Base::sampleWeight_);
        }
    };
 
    template <typename Numeric, class BwFcn, class WFcn>
    class LOrPE1DRLCVFunctorHelper : public LOrPE1DCVFunctorHelper<Numeric, BwFcn, WFcn>
    {
    public:
        typedef LOrPE1DCVFunctorHelper<Numeric, BwFcn, WFcn> Base;
 
        inline LOrPE1DRLCVFunctorHelper(
            const int i_symbetaPower,
            const double i_leftBoundary,
            const double i_rightBoundary,
            const BoundaryHandling& i_bh,
            const BwFcn& i_bandwidthFunctor,
            const WFcn& i_localizingWeight,
            const double i_localizingWeightXmin,
            const double i_localizingWeightXmax,
            const unsigned i_nIntegIntervals,
            const unsigned i_nIntegPoints,
            const bool i_normalizeLOrPEEstimate,
            const std::vector<Numeric>& i_coords,
            const double i_rlcvAlpha = 0.5)
            : Base(
                i_symbetaPower, i_leftBoundary, i_rightBoundary,
                i_bh, i_bandwidthFunctor, i_localizingWeight,
                i_localizingWeightXmin, i_localizingWeightXmax,
                i_nIntegIntervals, i_nIntegPoints,
                i_normalizeLOrPEEstimate, i_coords)
        {
            setRlcvAlpha(i_rlcvAlpha);
        }
 
        inline virtual ~LOrPE1DRLCVFunctorHelper() {}
 
        inline void setRlcvAlpha(const double alpha)
        {
            rlcvAlpha_ = alpha;
            minDensFactor_ = std::pow(Base::effSampleSize_, alpha);
        }
 
        inline virtual void setWeights(const double* weights, unsigned long nWeights)
        {
            Base::setWeights(weights, nWeights);
            minDensFactor_ = std::pow(Base::effSampleSize_, rlcvAlpha_);
        }
 
        inline double rlcvAlpha() const {return rlcvAlpha_;}
 
    protected:
        inline virtual double pointContribution(
            const double selfC, const double leaveOneOutEstimate,
            const double pointWeight, const double localWeight) const
        {
            const double minDens = selfC/minDensFactor_;
            return pointWeight*localWeight*std::log(std::max(leaveOneOutEstimate, minDens));
        }
 
        inline virtual double valueToMaximize(const double sum) const
        {
            return sum;
        }
 
    private:
        double rlcvAlpha_;
        double minDensFactor_;
    };
 
    /** A convenience function for creating LOrPE1DLSCVFunctorHelper objects */
    template <typename Numeric, class BwFcn, class WFcn>
    inline LOrPE1DLSCVFunctorHelper<Numeric, BwFcn, WFcn>
    LOrPE1DLSCVFunctor(const int i_symbetaPower,
                       const double i_leftBoundary,
                       const double i_rightBoundary,
                       const BoundaryHandling& i_bh,
                       const BwFcn& i_bandwidthFunctor,
                       const WFcn& i_localizingWeight,
                       const double i_localizingWeightXmin,
                       const double i_localizingWeightXmax,
                       const unsigned i_nIntegIntervals,
                       const unsigned i_nIntegPoints,
                       const bool i_normalizeLOrPEEstimate,
                       const std::vector<Numeric>& i_coords)
    {
        return LOrPE1DLSCVFunctorHelper<Numeric, BwFcn, WFcn>(
            i_symbetaPower, i_leftBoundary, i_rightBoundary,
            i_bh, i_bandwidthFunctor, i_localizingWeight,
            i_localizingWeightXmin, i_localizingWeightXmax,
            i_nIntegIntervals, i_nIntegPoints,
            i_normalizeLOrPEEstimate, i_coords);
    }
 
    /** 
    // A convenience function for creating LOrPE1DLSCVFunctorHelper objects
    // that perform global cross-validation (i.e., without CV localization
    // function).
    */
    template <typename Numeric, class BwFcn>
    inline LOrPE1DLSCVFunctorHelper<Numeric, BwFcn, double>
    LOrPE1DGlobLSCVFunctor(const int i_symbetaPower,
                           const double i_leftBoundary,
                           const double i_rightBoundary,
                           const BoundaryHandling& i_bh,
                           const BwFcn& i_bandwidthFunctor,
                           const unsigned i_nIntegIntervals,
                           const unsigned i_nIntegPoints,
                           const bool i_normalizeLOrPEEstimate,
                           const std::vector<Numeric>& i_coords)
    {
        return LOrPE1DLSCVFunctorHelper<Numeric, BwFcn, double>(
            i_symbetaPower, i_leftBoundary, i_rightBoundary,
            i_bh, i_bandwidthFunctor,
            1.0, -DBL_MAX, DBL_MAX,
            i_nIntegIntervals, i_nIntegPoints,
            i_normalizeLOrPEEstimate, i_coords);
    }
 
    /** 
    // A convenience function for creating LOrPE1DLSCVFunctorHelper objects
    // that perform global cross-validation (i.e., without CV localization
    // function) and use constant bandwidth.
    */
    template <typename Numeric>
    inline LOrPE1DLSCVFunctorHelper<Numeric, double, double>
    LOrPE1DSimpleLSCVFunctor(const int i_symbetaPower,
                             const double i_leftBoundary,
                             const double i_rightBoundary,
                             const BoundaryHandling& i_bh,
                             const unsigned i_nIntegIntervals,
                             const unsigned i_nIntegPoints,
                             const bool i_normalizeLOrPEEstimate,
                             const std::vector<Numeric>& i_coords)
    {
        return LOrPE1DLSCVFunctorHelper<Numeric, double, double>(
            i_symbetaPower, i_leftBoundary, i_rightBoundary,
            i_bh, 1.0, 1.0, -DBL_MAX, DBL_MAX,
            i_nIntegIntervals, i_nIntegPoints,
            i_normalizeLOrPEEstimate, i_coords);
    }
 
    /** A convenience function for creating LOrPE1DRLCVFunctorHelper objects */
    template <typename Numeric, class BwFcn, class WFcn>
    inline LOrPE1DRLCVFunctorHelper<Numeric, BwFcn, WFcn>
    LOrPE1DRLCVFunctor(const int i_symbetaPower,
                       const double i_leftBoundary,
                       const double i_rightBoundary,
                       const BoundaryHandling& i_bh,
                       const BwFcn& i_bandwidthFunctor,
                       const WFcn& i_localizingWeight,
                       const double i_localizingWeightXmin,
                       const double i_localizingWeightXmax,
                       const unsigned i_nIntegIntervals,
                       const unsigned i_nIntegPoints,
                       const bool i_normalizeLOrPEEstimate,
                       const std::vector<Numeric>& i_coords,
                       const double i_rlcvAlpha = 0.5)
    {
        return LOrPE1DRLCVFunctorHelper<Numeric, BwFcn, WFcn>(
            i_symbetaPower, i_leftBoundary, i_rightBoundary,
            i_bh, i_bandwidthFunctor, i_localizingWeight,
            i_localizingWeightXmin, i_localizingWeightXmax,
            i_nIntegIntervals, i_nIntegPoints,
            i_normalizeLOrPEEstimate, i_coords, i_rlcvAlpha);
    }
 
    /** 
    // A convenience function for creating LOrPE1DRLCVFunctorHelper objects
    // that perform global cross-validation (i.e., without CV localization
    // function).
    */
    template <typename Numeric, class BwFcn>
    inline LOrPE1DRLCVFunctorHelper<Numeric, BwFcn, double>
    LOrPE1DGlobRLCVFunctor(const int i_symbetaPower,
                           const double i_leftBoundary,
                           const double i_rightBoundary,
                           const BoundaryHandling& i_bh,
                           const BwFcn& i_bandwidthFunctor,
                           const unsigned i_nIntegIntervals,
                           const unsigned i_nIntegPoints,
                           const bool i_normalizeLOrPEEstimate,
                           const std::vector<Numeric>& i_coords,
                           const double i_rlcvAlpha = 0.5)
    {
        return LOrPE1DRLCVFunctorHelper<Numeric, BwFcn, double>(
            i_symbetaPower, i_leftBoundary, i_rightBoundary,
            i_bh, i_bandwidthFunctor,
            1.0, -DBL_MAX, DBL_MAX,
            i_nIntegIntervals, i_nIntegPoints,
            i_normalizeLOrPEEstimate, i_coords, i_rlcvAlpha);
    }
 
    /** 
    // A convenience function for creating LOrPE1DRLCVFunctorHelper objects
    // that perform global cross-validation (i.e., without CV localization
    // function) and use constant bandwidth.
    */
    template <typename Numeric>
    inline LOrPE1DRLCVFunctorHelper<Numeric, double, double>
    LOrPE1DSimpleRLCVFunctor(const int i_symbetaPower,
                             const double i_leftBoundary,
                             const double i_rightBoundary,
                             const BoundaryHandling& i_bh,
                             const unsigned i_nIntegIntervals,
                             const unsigned i_nIntegPoints,
                             const bool i_normalizeLOrPEEstimate,
                             const std::vector<Numeric>& i_coords,
                             const double i_rlcvAlpha = 0.5)
    {
        return LOrPE1DRLCVFunctorHelper<Numeric, double, double>(
            i_symbetaPower, i_leftBoundary, i_rightBoundary,
            i_bh, 1.0, 1.0, -DBL_MAX, DBL_MAX,
            i_nIntegIntervals, i_nIntegPoints,
            i_normalizeLOrPEEstimate, i_coords, i_rlcvAlpha);
    }
 
    // Convenience functions for creating either RLCV or LSCV
    // functors. The convention is that LSCV will be used in case
    // rlcvAlpha parameter is negative.
    template <typename Numeric, class BwFcn, class WFcn>
    inline std::unique_ptr<LOrPE1DCVFunctorHelper<Numeric, BwFcn, WFcn> >
    LOrPE1DCVFunctor(const int i_symbetaPower,
                     const double i_leftBoundary,
                     const double i_rightBoundary,
                     const BoundaryHandling& i_bh,
                     const BwFcn& i_bandwidthFunctor,
                     const WFcn& i_localizingWeight,
                     const double i_localizingWeightXmin,
                     const double i_localizingWeightXmax,
                     const unsigned i_nIntegIntervals,
                     const unsigned i_nIntegPoints,
                     const bool i_normalizeLOrPEEstimate,
                     const std::vector<Numeric>& i_coords,
                     const double i_rlcvAlpha)
    {
        typedef LOrPE1DCVFunctorHelper<Numeric, BwFcn, WFcn> Base;
        typedef LOrPE1DRLCVFunctorHelper<Numeric, BwFcn, WFcn> RLCV;
        typedef LOrPE1DLSCVFunctorHelper<Numeric, BwFcn, WFcn> LSCV;
 
        Base* ptr = 0;
        if (i_rlcvAlpha < 0.0)
            ptr = new LSCV(
                i_symbetaPower, i_leftBoundary, i_rightBoundary,
                i_bh, i_bandwidthFunctor, i_localizingWeight,
                i_localizingWeightXmin, i_localizingWeightXmax,
                i_nIntegIntervals, i_nIntegPoints,
                i_normalizeLOrPEEstimate, i_coords);
        else
            ptr = new RLCV(
                i_symbetaPower, i_leftBoundary, i_rightBoundary,
                i_bh, i_bandwidthFunctor, i_localizingWeight,
                i_localizingWeightXmin, i_localizingWeightXmax,
                i_nIntegIntervals, i_nIntegPoints,
                i_normalizeLOrPEEstimate, i_coords, i_rlcvAlpha);
        return std::unique_ptr<Base>(ptr);
    }
 
    template <typename Numeric, class BwFcn>
    inline std::unique_ptr<LOrPE1DCVFunctorHelper<Numeric, BwFcn, double> >
    LOrPE1DGlobCVFunctor(const int i_symbetaPower,
                         const double i_leftBoundary,
                         const double i_rightBoundary,
                         const BoundaryHandling& i_bh,
                         const BwFcn& i_bandwidthFunctor,
                         const unsigned i_nIntegIntervals,
                         const unsigned i_nIntegPoints,
                         const bool i_normalizeLOrPEEstimate,
                         const std::vector<Numeric>& i_coords,
                         const double i_rlcvAlpha)
    {
        typedef LOrPE1DCVFunctorHelper<Numeric, BwFcn, double> Base;
        typedef LOrPE1DRLCVFunctorHelper<Numeric, BwFcn, double> RLCV;
        typedef LOrPE1DLSCVFunctorHelper<Numeric, BwFcn, double> LSCV;
 
        Base* ptr = 0;
        if (i_rlcvAlpha < 0.0)
            ptr = new LSCV(
                i_symbetaPower, i_leftBoundary, i_rightBoundary,
                i_bh, i_bandwidthFunctor,
                1.0, -DBL_MAX, DBL_MAX,
                i_nIntegIntervals, i_nIntegPoints,
                i_normalizeLOrPEEstimate, i_coords);
        else
            ptr = new RLCV(
                i_symbetaPower, i_leftBoundary, i_rightBoundary,
                i_bh, i_bandwidthFunctor,
                1.0, -DBL_MAX, DBL_MAX,
                i_nIntegIntervals, i_nIntegPoints,
                i_normalizeLOrPEEstimate, i_coords, i_rlcvAlpha);
        return std::unique_ptr<Base>(ptr);
    }
 
    template <typename Numeric>
    inline std::unique_ptr<LOrPE1DCVFunctorHelper<Numeric, double, double> >
    LOrPE1DSimpleCVFunctor(const int i_symbetaPower,
                           const double i_leftBoundary,
                           const double i_rightBoundary,
                           const BoundaryHandling& i_bh,
                           const unsigned i_nIntegIntervals,
                           const unsigned i_nIntegPoints,
                           const bool i_normalizeLOrPEEstimate,
                           const std::vector<Numeric>& i_coords,
                           const double i_rlcvAlpha)
    {
        typedef LOrPE1DCVFunctorHelper<Numeric, double, double> Base;
        typedef LOrPE1DRLCVFunctorHelper<Numeric, double, double> RLCV;
        typedef LOrPE1DLSCVFunctorHelper<Numeric, double, double> LSCV;
 
        Base* ptr = 0;
        if (i_rlcvAlpha < 0.0)
            ptr = new LSCV(
                i_symbetaPower, i_leftBoundary, i_rightBoundary,
                i_bh, 1.0, 1.0, -DBL_MAX, DBL_MAX,
                i_nIntegIntervals, i_nIntegPoints,
                i_normalizeLOrPEEstimate, i_coords);
        else
            ptr = new RLCV(
                i_symbetaPower, i_leftBoundary, i_rightBoundary,
                i_bh, 1.0, 1.0, -DBL_MAX, DBL_MAX,
                i_nIntegIntervals, i_nIntegPoints,
                i_normalizeLOrPEEstimate, i_coords, i_rlcvAlpha);
        return std::unique_ptr<Base>(ptr);
    }
}
 
#endif // NPSTAT_LORPE1DCV_HH_