ContOrthoPoly1DQ.hh

Go to the documentation of this file.

#ifndef NPSTAT_CONTORTHOPOLY1DQ_HH_
#define NPSTAT_CONTORTHOPOLY1DQ_HH_
 
/*!
// \file ContOrthoPoly1DQ.hh
//
// \brief Continuous orthogonal polynomial series in one dimension
//        generated for a discrete measure, using quadruple precision
//
// Author: I. Volobouev
//
// July 2023
*/
 
#include <vector>
#include <utility>
 
#include "npstat/nm/OrthoPolyMethod.hh"
#include "npstat/nm/RecurrenceQ.hh"
#include "npstat/nm/Matrix.hh"
 
namespace npstat {
    class ContOrthoPoly1DQ
    {
    public:
        typedef std::pair<quad_float,quad_float> MeasurePoint;
 
        /**
        // Main constructor, with obvious arguments. The first element
        // of the measure pair is the coordinate and the second element
        // is the weight (all weights must be non-negative). Internally,
        // the measure will be sorted in the order of increasing weight
        // and the measure coordinates will normally be shifted so that
        // their weighted mean is at 0.
        */
        template<typename Numeric1, typename Numeric2>
        ContOrthoPoly1DQ(unsigned maxDegree,
                         const std::vector<std::pair<Numeric1,Numeric2> >& measure,
                         OrthoPolyMethod m = OPOLY_STIELTJES,
                         bool shiftMeasureCoords = true);
 
        /** Constructor which assumes that all weights are 1.0 */
        template<typename Numeric>
        ContOrthoPoly1DQ(unsigned maxDegree,
                         const std::vector<Numeric>& coords,
                         OrthoPolyMethod m = OPOLY_STIELTJES,
                         bool shiftMeasureCoords = true);
 
        /** Constructor that uses the modified Chebyshev algorithm */
        template<typename Numeric1, typename Numeric2, class OrthoPoly>
        ContOrthoPoly1DQ(unsigned maxDegree,
                         const std::vector<std::pair<Numeric1,Numeric2> >& measure,
                         const OrthoPoly& ops, bool shiftMeasureCoords = false);
 
        /**
        // Constructor that uses the modified Chebyshev algorithm
        // and assumes that all weights are 1.0
        */
        template<typename Numeric, class OrthoPoly>
        ContOrthoPoly1DQ(unsigned maxDegree,
                         const std::vector<Numeric>& coords,
                         const OrthoPoly& ops, bool shiftMeasureCoords = false);
 
        /**
        // Function added for interface compatibility with
        // AbsClassicalOrthoPoly1D
        */
        inline ContOrthoPoly1DQ* clone() const {return new ContOrthoPoly1DQ(*this);}
 
        //@{
        /** A simple inspector of object properties */
        inline unsigned maxDegree() const {return maxdeg_;}
        inline bool areAllWeightsEqual() const {return allWeightsEqual_;}
        inline quad_float measureShift() const {return shiftX_;}
        inline quad_float meanCoordinate() const {return meanX_;}
        inline quad_float location() const {return meanX_;}
        inline quad_float scale() const {return 1.0;}
        //@}
 
        /** Kish's effective sample size for the measure */
        quad_float effectiveSampleSize() const;
 
        /** Return the value of one of the normalized polynomials */
        quad_float poly(unsigned deg, quad_float x) const;
 
        /** Return the values of two of the normalized polynomials */
        std::pair<quad_float,quad_float> polyPair(
            unsigned deg1, unsigned deg2, quad_float x) const;
 
        /**
        // Return the values of all orthonormal polynomials up to some
        // maximum degree. Size of the "polyValues" array should be
        // at least maxdeg + 1.
        */
        template <typename Numeric>
        inline void allPolys(const quad_float x, const unsigned maxdeg,
                             Numeric* polyValues) const
            {getAllPolys(x - shiftX_, maxdeg, polyValues);}
 
        /**
        // Function added for interface compatibility with
        // AbsClassicalOrthoPoly1D
        */
        inline void allpoly(const long double x, long double* values,
                            const unsigned maxdeg) const
            {getAllPolys(x - shiftX_, maxdeg, values);}
 
        /**
        // This method is useful for testing the numerical precision
        // of the orthonormalization procedure. It returns the scalar
        // products between various polynomials.
        */
        quad_float empiricalKroneckerDelta(unsigned deg1, unsigned deg2) const;
 
        /**
        // A faster way to generate Kronecker deltas if the whole matrix
        // of them is needed. The argument must not exceed the return
        // value of the "maxDegree" function. The resulting matrix will
        // be dimensioned (maxdeg + 1) x (maxdeg + 1).
        */
        Matrix<quad_float> empiricalKroneckerMatrix(unsigned maxdeg) const;
 
    private:
        template<typename Numeric1, typename Numeric2>
        void copyMeasure(const std::vector<std::pair<Numeric1,Numeric2> >& inMeasure);
 
        template<typename Numeric>
        void copyInputCoords(const std::vector<Numeric>& inCoords);
 
        void calcMeanXandWeightSums(bool shiftTheMeasure);
 
        template<class OrthoPoly>
        quad_float sigma(const OrthoPoly& ops, int k, int l,
                         Matrix<quad_float>& sigMatrix,
                         Matrix<int>& flagMatrix) const;
 
        template<class OrthoPoly>
        void calculateMonicAverages(const OrthoPoly& ops, quad_float* averages,
                                    unsigned maxdeg);
 
        template<class OrthoPoly>
        void calcRecurrenceModifiedChebyshev(const OrthoPoly& ops);
 
        quad_float normpoly(unsigned degree, quad_float x) const;
 
        std::pair<quad_float,quad_float> twonormpoly(
            unsigned deg1, unsigned deg2, quad_float x) const;
 
        template <typename Numeric>
        void getAllPolys(quad_float x, unsigned maxdeg, Numeric *polyValues) const;
 
        std::vector<MeasurePoint> measure_;
        std::vector<Private::RecurrenceQ> rCoeffs_;
        quad_float wsum_;
        quad_float wsumsq_;
        quad_float meanX_;
        quad_float shiftX_;
        quad_float minX_;
        quad_float maxX_;
        unsigned maxdeg_;
        bool allWeightsEqual_;
    };
}
 
#include "npstat/nm/ContOrthoPoly1DQ.icc"
 
#endif // NPSTAT_CONTORTHOPOLY1DQ_HH_