Index: trunk/npstat/nm/ScalableClassicalOrthoPoly1D.hh =================================================================== --- trunk/npstat/nm/ScalableClassicalOrthoPoly1D.hh (revision 722) +++ trunk/npstat/nm/ScalableClassicalOrthoPoly1D.hh (revision 723) @@ -1,175 +1,233 @@ #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 ScalableClassicalOrthoPoly1D { public: ScalableClassicalOrthoPoly1D(const AbsClassicalOrthoPoly1D& poly, long double location, long double scale); ScalableClassicalOrthoPoly1D(const ScalableClassicalOrthoPoly1D&); ScalableClassicalOrthoPoly1D& operator=(const ScalableClassicalOrthoPoly1D&); inline ~ScalableClassicalOrthoPoly1D() {delete poly_;} inline long double location() const {return location_;} inline long double scale() const {return scale_;} 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_);} /** // Values of two orthonormal polynomials. // Faster than calling "poly" two times. */ inline std::pair 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 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 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 inline void sampleCoeffs(const Numeric* coords, const unsigned long lenCoords, double* coeffs, const unsigned maxdeg) const { const Numeric mu = static_cast(location_); const Numeric s = static_cast(scale_); poly_->sampleCoeffs(coords, lenCoords, mu, s, coeffs, maxdeg); } template inline void sampleCoeffVars(const Numeric* coords, const unsigned long lenCoords, const double* coeffs, const unsigned maxdeg, double* variances) const { const Numeric mu = static_cast(location_); const Numeric s = static_cast(scale_); poly_->sampleCoeffVars(coords, lenCoords, mu, s, coeffs, maxdeg, variances); } template inline npstat::Matrix sampleCoeffCovariance(const Numeric* coords, const unsigned long lenCoords, const double* coeffs, const unsigned maxdeg) const { const Numeric mu = static_cast(location_); const Numeric s = static_cast(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 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(location_); const Numeric s = static_cast(scale_); return poly_->weightedSampleCoeffs(points, numPoints, mu, s, coeffs, maxdeg); } template 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(location_); const Numeric s = static_cast(scale_); return poly_->weightedSampleCoeffVars(points,nPoints,mu,s,coeffs,maxdeg,variances); } template inline npstat::Matrix 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(location_); const Numeric s = static_cast(scale_); return poly_->weightedSampleCoeffCovariance(points, nPoints, mu, s, coeffs, maxdeg); } private: AbsClassicalOrthoPoly1D* poly_; long double location_; long double scale_; }; + + /** + // 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 + { + public: + inline explicit ScalableOrthoPoly1DWeight( + const ScalableClassicalOrthoPoly1D& fcn, + const long double normfactor=1.0L) + : fcn_(fcn), norm_(normfactor) {} + + 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 + { + 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_