ClassicalOrthoPoly1DFromWeight.hh

Go to the documentation of this file.

#ifndef NPSTAT_CLASSICALORTHOPOLY1DFROMWEIGHT_HH_
#define NPSTAT_CLASSICALORTHOPOLY1DFROMWEIGHT_HH_
 
/*!
// \file ClassicalOrthoPoly1DFromWeight.hh
//
// \brief Orthogonal polynomials on an interval with an arbitrary
//        user-provided weight function
//
// Author: I. Volobouev
//
// November 2020
*/
 
#include <vector>
 
#include "npstat/nm/AbsClassicalOrthoPoly1D.hh"
#include "npstat/nm/Recurrence.hh"
#include "npstat/nm/OrthoPolyMethod.hh"
 
namespace npstat {
    template <class Functor1D>
    class ClassicalOrthoPoly1DFromWeight : public AbsClassicalOrthoPoly1D
    {
    public:
        /**
        // "nIntegPoints" parameter is the number of integration points
        // for the Fejer or rectangular rule quadrature. Note that Fejer
        // quadrature with n points can integrate exactly polynomials of
        // degree n-1 and that this quadrature will be used for calculating
        // various inner products in which the density is used as the weight.
        // Therefore, "nIntegPoints" should be at least 2*maxDegree + 1 + C,
        // where "C" is the polynomial degree that can be used to approximate
        // the density reasonably well.
        //
        // Rectangular rule quadrature can be used for densities whose
        // multiple derivatives become very close to 0 near the support
        // boundaries. In particular, rectangular quadrature is expected
        // to be useful for KDE densities estimated with symmetric beta
        // kernels without boundary correction and with support on
        // an appropriately extended interval.
        */
        inline ClassicalOrthoPoly1DFromWeight(
            const Functor1D& fcn,
            const unsigned maxDegree, const unsigned nIntegPoints,
            const long double i_xmin, const long double i_xmax,
            const OrthoPolyMethod m = OPOLY_STIELTJES,
            const bool useFejerQuadrature = true)
            : fcn_(fcn), xmin_(i_xmin), xmax_(i_xmax), maxdeg_(maxDegree)
        {
            initialize(nIntegPoints, m, useFejerQuadrature);
        }
 
        inline virtual ~ClassicalOrthoPoly1DFromWeight() {}
 
        inline virtual ClassicalOrthoPoly1DFromWeight* clone() const
            {return new ClassicalOrthoPoly1DFromWeight(*this);}
 
        inline virtual long double weight(const long double x) const
            {return fcn_(x)/wsum_;}
        inline virtual double xmin() const {return xmin_;}
        inline virtual double xmax() const {return xmax_;}
        inline virtual unsigned maxDegree() const {return maxdeg_;}
 
        inline virtual std::pair<long double,long double>
        monicRecurrenceCoeffs(const unsigned deg) const
        {
            const Private::Recurrence& r = rCoeffs_.at(deg);
            return std::pair<long double,long double>(r.alpha, r.beta);
        }
 
        inline virtual std::pair<long double,long double>
        recurrenceCoeffs(const unsigned deg) const
        {
            const Private::Recurrence& r = rCoeffs_.at(deg);
            return std::pair<long double,long double>(r.alpha, r.sqrbeta);
        }
 
    private:
        void initialize(unsigned nIntegPoints, OrthoPolyMethod m,
                        bool useFejerQuadrature);
 
        template<class Quad1D>
        void calcCoeffs(const Quad1D& quad, OrthoPolyMethod m);
 
        Functor1D fcn_;
        std::vector<Private::Recurrence> rCoeffs_;
        long double wsum_;
        long double xmin_;
        long double xmax_;
        unsigned maxdeg_;
 
#ifdef SWIG
    public:
        inline ClassicalOrthoPoly1DFromWeight(
            const Functor1D& fcn,
            const unsigned maxDegree, const unsigned nIntegPoints,
            const double i_xmin, const double i_xmax,
            const OrthoPolyMethod m = OPOLY_STIELTJES,
            const bool useFejerQuadrature = true)
            : fcn_(fcn), xmin_(i_xmin), xmax_(i_xmax), maxdeg_(maxDegree)
        {
            initialize(nIntegPoints, m, useFejerQuadrature);
        }
#endif
    };
}
 
#include "npstat/nm/ClassicalOrthoPoly1DFromWeight.icc"
 
#endif // NPSTAT_CLASSICALORTHOPOLY1DFROMWEIGHT_HH_