ScalableClassicalOrthoPoly1D.hh

Go to the documentation of this file.

#ifndef NPSTAT_SCALABLECLASSICALORTHOPOLY1D_HH_
#define NPSTAT_SCALABLECLASSICALORTHOPOLY1D_HH_
 
/*!
// \file ScalableClassicalOrthoPoly1D.hh
//
// \brief Class for shifting and scaling classical continuous orthonormal
//        polynomials
//
// Author: I. Volobouev
//
// October 2020
*/
 
#include "npstat/nm/AbsClassicalOrthoPoly1D.hh"
 
namespace npstat {
    class ScalableOrthoPoly1DDeg;
    class ScalableOrthoPoly1DWeight;
 
    class ScalableClassicalOrthoPoly1D
    {
    public:
        typedef ScalableOrthoPoly1DDeg degree_functor;
        typedef ScalableOrthoPoly1DWeight weight_functor;
 
        ScalableClassicalOrthoPoly1D(const AbsClassicalOrthoPoly1D& poly,
                                     long double location, long double scale);
 
        ScalableClassicalOrthoPoly1D(const AbsClassicalOrthoPoly1D& poly,
                                     double location, double scale);
 
        ScalableClassicalOrthoPoly1D(const ScalableClassicalOrthoPoly1D&);
        ScalableClassicalOrthoPoly1D& operator=(const ScalableClassicalOrthoPoly1D&);
 
        inline ScalableClassicalOrthoPoly1D* clone() const
            {return new ScalableClassicalOrthoPoly1D(*this);}
 
        inline ~ScalableClassicalOrthoPoly1D() {delete poly_;}
 
        inline long double location() const {return location_;}
        inline long double scale() const {return scale_;}
 
        // Note that the Jacobi matrix principal minor is returned
        // for the unscaled OPS
        inline Matrix<long double> jacobiMatrix(const unsigned n) const
            {return poly_->jacobiMatrix(n);}
 
        inline void setLocation(const long double v) {location_ = v;}
        inline void setScale(const long double v)
        {
            if (v <= 0.0) throw std::invalid_argument(
                "In npstat::ScalableClassicalOrthoPoly1D::setScale: "
                "scale parameter must be positive");
            scale_ = v;
        }
 
        inline long double weight(const long double x) const
            {return poly_->weight((x - location_)/scale_)/scale_;}
 
        inline double xmin() const {return scale_*poly_->xmin() + location_;}
 
        inline double xmax() const {return scale_*poly_->xmax() + location_;}
 
        inline unsigned maxDegree() const {return poly_->maxDegree();}
 
        inline long double poly(const unsigned deg, const long double x) const
            {return poly_->poly(deg, (x - location_)/scale_);}
 
        inline long double monic(const unsigned deg, const long double x) const
            {return poly_->monic(deg, (x - location_)/scale_);}
 
        /** 
        // Values of two orthonormal polynomials.
        // Faster than calling "poly" two times.
        */
        inline std::pair<long double,long double> twopoly(
            const unsigned deg1, const unsigned deg2, const long double x) const
            {return poly_->twopoly(deg1, deg2, (x - location_)/scale_);}
 
        /**
        // Values of all orthonormal polynomials up to some degree.
        // Faster than calling "poly" multiple times. The size of
        // the "values" array should be at least maxdeg + 1.
        */
        inline void allpoly(const long double x,
                            long double* values, const unsigned maxdeg) const
            {poly_->allpoly((x - location_)/scale_, values, maxdeg);}
 
        inline void allmonic(const long double x,
                            long double* values, const unsigned maxdeg) const
            {poly_->allmonic((x - location_)/scale_, values, maxdeg);}
 
        inline double series(const double* coeffs, const unsigned maxdeg,
                             const double x) const
            {return poly_->series(coeffs, maxdeg, (x - location_)/scale_);}
 
        inline double integrateSeries(const double* coeffs, const unsigned maxdeg) const
            {return scale_*poly_->integrateSeries(coeffs, maxdeg);}
 
        /**
        // Build the coefficients of the orthogonal polynomial series
        // for the given function. The length of the array "coeffs"
        // should be at least maxdeg + 1. Note that the coefficients
        // are returned in the order of increasing degree (same order
        // is used by the "series" function).
        */
        template <class Functor, class Quadrature>
        void calculateCoeffs(const Functor& fcn, const Quadrature& quad,
                             double* coeffs, unsigned maxdeg) const;
 
        /**
        // Build the coefficients of the orthogonal polynomial series
        // for the given sample of points (empirical density function).
        // The length of the array "coeffs" should be at least maxdeg + 1.
        // Note that the coefficients are returned in the order of
        // increasing degree (same order is used by the "series" function).
        */
        template <class Numeric>
        inline void sampleCoeffs(const Numeric* coords, const unsigned long lenCoords,
                                 double* coeffs, const unsigned maxdeg) const
        {
            const Numeric mu = static_cast<Numeric>(location_);
            const Numeric s = static_cast<Numeric>(scale_);
            poly_->sampleCoeffs(coords, lenCoords, mu, s, coeffs, maxdeg);
        }
 
        template <class Numeric>
        inline void sampleCoeffVars(const Numeric* coords, const unsigned long lenCoords,
                                    const double* coeffs, const unsigned maxdeg,
                                    double* variances) const
        {
            const Numeric mu = static_cast<Numeric>(location_);
            const Numeric s = static_cast<Numeric>(scale_);
            poly_->sampleCoeffVars(coords, lenCoords, mu, s, coeffs, maxdeg, variances);
        }
 
        template <class Numeric>
        inline npstat::Matrix<double>
        sampleCoeffCovariance(const Numeric* coords, const unsigned long lenCoords,
                              const double* coeffs, const unsigned maxdeg) const
        {
            const Numeric mu = static_cast<Numeric>(location_);
            const Numeric s = static_cast<Numeric>(scale_);
            return poly_->sampleCoeffCovariance(coords, lenCoords, mu, s, coeffs, maxdeg);
        }
 
        /**
        // Build the coefficients of the orthogonal polynomial series for
        // the given sample of weighted points (empirical density function).
        // The first element of the pair will be treated as the coordinate
        // and the second element will be treated as weight. Weights must
        // not be negative. The length of the array "coeffs" should be
        // at least maxdeg + 1. Note that the coefficients are returned
        // in the order of increasing degree (same order is used by the
        // "series" function).
        */
        template <class Pair>
        inline void weightedSampleCoeffs(const Pair* points, const unsigned long numPoints,
                                         double* coeffs, const unsigned maxdeg) const
        {
            typedef typename Pair::first_type Numeric;
            const Numeric mu = static_cast<Numeric>(location_);
            const Numeric s = static_cast<Numeric>(scale_);
            return poly_->weightedSampleCoeffs(points, numPoints, mu, s, coeffs, maxdeg);
        }
 
        template <class Pair>
        inline void weightedSampleCoeffVars(const Pair* points, const unsigned long nPoints,
                                            const double* coeffs, const unsigned maxdeg,
                                            double* variances) const
        {
            typedef typename Pair::first_type Numeric;
            const Numeric mu = static_cast<Numeric>(location_);
            const Numeric s = static_cast<Numeric>(scale_);
            return poly_->weightedSampleCoeffVars(points,nPoints,mu,s,coeffs,maxdeg,variances);
        }
 
        template <class Pair>
        inline npstat::Matrix<double>
        weightedSampleCoeffCovariance(const Pair* points, const unsigned long nPoints,
                                      const double* coeffs, const unsigned maxdeg) const
        {
            typedef typename Pair::first_type Numeric;
            const Numeric mu = static_cast<Numeric>(location_);
            const Numeric s = static_cast<Numeric>(scale_);
            return poly_->weightedSampleCoeffCovariance(points, nPoints, mu, s, coeffs, maxdeg);
        }
 
        /**
        // Unweighted averages of the polynomial values for the given sample.
        // The length of array "averages" should be at least maxdeg + 1.
        */
        template <class Numeric, typename Real>
        void sampleAverages(const Numeric* coords, unsigned long lenCoords,
                            Real* averages, unsigned maxdeg) const;
 
        /**
        // Averages of the polynomial values for the given sample of
        // weighted points (not weighting by the polynomial weight function).
        // The first element of the pair will be treated as the coordinate
        // and the second element will be treated as weight. Weights must not
        // be negative. The length of array "averages" should be at least
        // maxdeg + 1.
        */
        template <class Pair, typename Real>
        void weightedPointsAverages(const Pair* points, unsigned long numPoints,
                                    Real* averages, unsigned maxdeg) const;
 
        /**
        // Averages of the polynomial values weighted by an external
        // weight function. The length of array "averages" should be
        // at least maxdeg + 1.
        */
        template <class Functor, typename Real>
        void extWeightAverages(const Functor& extWeight,
                               const AbsIntervalQuadrature1D& quad,
                               Real* averages, unsigned maxdeg) const;
 
        /**
        // Pairwise products of the polynomial values weighted by
        // an external weight function. The returned matrix will be
        // symmetric and will have the dimensions (maxdeg + 1) x (maxdeg + 1).
        */
        template <class Functor>
        Matrix<double> extWeightProductAverages(const Functor& extWeight,
            const AbsIntervalQuadrature1D& quad, unsigned maxdeg) const;
 
    private:
        AbsClassicalOrthoPoly1D* poly_;
        long double location_;
        long double scale_;
 
#ifdef SWIG
    public:
        inline void setLocation2(const double v)
            {setLocation(v);}
 
        inline void setScale2(const double v)
            {setScale(v);}
 
        inline double location2() const
            {return location();}
 
        inline double scale2() const
            {return scale();}
 
        inline double weight2(const double x) const
            {return weight(x);}
 
        inline double poly2(const unsigned deg, const double x) const
            {return poly(deg, x);}
 
        inline double monic2(const unsigned deg, const double x) const
            {return monic(deg, x);}
 
        inline std::pair<double,double> twopoly2(const unsigned deg1,
                                                 const unsigned deg2,
                                                 const double x) const
        {
            const std::pair<long double,long double>& p = twopoly(deg1, deg2, x);
            return std::pair<double,double>(p.first, p.second);
        }
 
        inline std::vector<double> allpoly2(const unsigned maxdeg, const double x) const
        {
            std::vector<long double> p(maxdeg+1);
            allpoly(x, &p[0], p.size());
            return std::vector<double>(p.begin(), p.end());
        }
 
        inline std::vector<double> allmonic2(const unsigned maxdeg, const double x) const
        {
            std::vector<long double> p(maxdeg+1);
            allmonic(x, &p[0], p.size());
            return std::vector<double>(p.begin(), p.end());
        }
#endif // SWIG
    };
 
    /**
    // A functor for the weight function of the given ortho poly system.
    // The poly system is not copied, only a reference is used. It is
    // a responsibility of the user to make sure that the lifetime of
    // the poly system object exceeds the lifetime of the functor.
    */
    class ScalableOrthoPoly1DWeight : public Functor1<long double, long double>
    {
    public:
        inline explicit ScalableOrthoPoly1DWeight(
            const ScalableClassicalOrthoPoly1D& fcn,
            const long double normfactor)
            : fcn_(fcn), norm_(normfactor) {}
 
        // Separate constructor here because swig gets confused
        // by long doubles with defaults
        inline explicit ScalableOrthoPoly1DWeight(
            const ScalableClassicalOrthoPoly1D& fcn)
            : fcn_(fcn), norm_(1.0L) {}
 
        inline virtual ~ScalableOrthoPoly1DWeight() {}
 
        inline virtual long double operator()(const long double& a) const
            {return norm_*fcn_.weight(a);}
 
    private:
        ScalableOrthoPoly1DWeight();
 
        const ScalableClassicalOrthoPoly1D& fcn_;
        long double norm_;
    };
 
    /**
    // A functor for a particular degree polynomial of the given ortho
    // poly system. The poly system is not copied, only a reference is used.
    // It is a responsibility of the user to make sure that the lifetime of
    // the poly system object exceeds the lifetime of the functor.
    */
    class ScalableOrthoPoly1DDeg : public Functor1<long double, long double>
    {
    public:
        inline ScalableOrthoPoly1DDeg(const ScalableClassicalOrthoPoly1D& fcn,
                                      const unsigned degree,
                                      const long double normfactor=1.0L)
            : fcn_(fcn), norm_(normfactor), deg_(degree)
        {
            if (deg_ > fcn_.maxDegree()) throw std::invalid_argument(
                "In npstat::ScalableOrthoPoly1DDeg::constructor: "
                "degree argument is out of range");
        }
 
        inline virtual ~ScalableOrthoPoly1DDeg() {}
 
        inline virtual long double operator()(const long double& a) const
            {return norm_*fcn_.poly(deg_, a);}
 
    private:
        ScalableOrthoPoly1DDeg();
 
        const ScalableClassicalOrthoPoly1D& fcn_;
        long double norm_;
        unsigned deg_;
    };
}
 
#include "npstat/nm/ScalableClassicalOrthoPoly1D.icc"
 
#endif // NPSTAT_SCALABLECLASSICALORTHOPOLY1D_HH_