|
Go to the documentation of this file. 1 #ifndef NPSTAT_CONTORTHOPOLY1D_HH_
2 #define NPSTAT_CONTORTHOPOLY1D_HH_
21 #include "geners/CPP11_auto_ptr.hh"
29 #ifdef USE_LAPACK_QUAD
36 class StorablePolySeries1D;
37 class ContOrthoPoly1DDeg;
42 typedef npstat::PreciseType<double>::type precise_type;
46 typedef std::pair<double,double> MeasurePoint;
48 #ifdef USE_LAPACK_QUAD
63 template< typename Numeric1, typename Numeric2>
65 const std::vector<std::pair<Numeric1,Numeric2> >& measure,
67 bool shiftMeasureCoords = true);
70 template< typename Numeric>
72 const std::vector<Numeric>& coords,
74 bool shiftMeasureCoords = true);
77 template< typename Numeric1, typename Numeric2, class OrthoPoly>
79 const std::vector<std::pair<Numeric1,Numeric2> >& measure,
80 const OrthoPoly& ops, bool shiftMeasureCoords = false);
91 template< typename Numeric, class OrthoPoly>
93 const std::vector<Numeric>& coords,
94 const OrthoPoly& ops, bool shiftMeasureCoords = false);
106 inline unsigned long measureLength() const { return measure_.size();}
107 inline double weight( const unsigned long index) const
108 { return measure_.at(index).second;}
109 inline precise_type sumOfWeights() const { return wsum_;}
110 inline precise_type sumOfWeightSquares() const { return wsumsq_;}
111 inline bool areAllWeightsEqual() const { return allWeightsEqual_;}
112 inline double minCoordinate() const { return minX_;}
113 inline double maxCoordinate() const { return maxX_;}
114 inline precise_type measureShift() const { return shiftX_;}
115 inline precise_type meanCoordinate() const { return meanX_;}
116 inline precise_type location() const { return meanX_;}
117 inline precise_type scale() const { return 1.0;}
125 inline MeasurePoint measure( const unsigned long index) const
126 { return measure_.at(index);}
127 inline double coordinate( const unsigned long index) const
128 { return measure_.at(index).first;}
135 double poly( unsigned deg, double x) const;
139 unsigned deg1, unsigned deg2, double x) const;
146 template < typename Numeric>
147 inline void allPolys( const double x, const unsigned maxdeg,
148 Numeric* polyValues) const
149 {getAllPolys(x - shiftX_, maxdeg, polyValues);}
155 inline void allpoly( const long double x, long double* values,
156 const unsigned maxdeg) const
157 {getAllPolys(x - shiftX_, maxdeg, values);}
164 template < typename Numeric>
166 double* averages, unsigned maxdeg) const;
176 template < class Pair>
178 double* averages, unsigned maxdeg) const;
186 template < typename Numeric>
188 unsigned long lenCoords,
189 unsigned maxdeg) const;
192 template < typename Numeric>
193 double series( const Numeric* coeffs, unsigned deg, double x) const;
202 template < class Functor, typename Numeric>
204 Numeric* coeffs, unsigned maxdeg) const;
221 template < class Functor>
223 double *coeffs, unsigned maxdeg,
224 double xmin, double xmax,
225 double integralValue = 1.0) const;
244 template < class Functor>
246 const double *coeffs, unsigned maxdeg,
247 double xmin, double xmax,
248 unsigned integrationRulePoints=0U) const;
254 template < class Functor>
256 unsigned maxdeg) const;
280 unsigned deg3, unsigned deg4) const;
297 double xmin, double xmax) const;
304 unsigned deg3, double xmin,
312 const std::vector<std::pair<unsigned,unsigned> >& pairs,
313 unsigned maxdeg, double xmin, double xmax) const;
324 template < class Functor>
326 unsigned deg1, unsigned deg2,
327 double xmin, double xmax,
328 unsigned integrationRulePoints) const;
335 template < class Functor>
338 unsigned maxdeg, double xmin, double xmax,
339 double* averages) const;
345 template < class Functor>
349 double xmin, double xmax) const;
362 unsigned deg2, unsigned power2) const;
369 bool degreesSortedIncreasingOrder = false) const;
378 unsigned deg3, unsigned deg4) const;
380 unsigned deg3, unsigned deg4,
381 unsigned deg5, unsigned deg6) const;
383 unsigned deg3, unsigned deg4,
384 unsigned deg5, unsigned deg6,
385 unsigned deg7, unsigned deg8) const;
389 double cov4( unsigned a, unsigned b, unsigned c, unsigned d) const;
392 double cov6( unsigned a, unsigned b, unsigned c, unsigned d,
393 unsigned e, unsigned f) const;
396 double cov8( unsigned a, unsigned b, unsigned c, unsigned d,
397 unsigned e, unsigned f, unsigned g, unsigned h) const;
400 double covCov4( unsigned a, unsigned b, unsigned c, unsigned d,
401 unsigned e, unsigned f, unsigned g, unsigned h) const;
407 double slowCov8( unsigned a, unsigned b, unsigned c, unsigned d,
408 unsigned e, unsigned f, unsigned g, unsigned h) const;
422 unsigned m, unsigned n, bool highOrder) const;
428 const std::vector<std::pair<unsigned,unsigned> >& pairs,
429 bool highOrder) const;
433 double xmin, double xmax,
434 const double *coeffs=0, unsigned maxdeg=0) const;
437 inline std::pair<precise_type,precise_type>
440 const Recurrence& rc(rCoeffs_.at(k));
441 return std::pair<precise_type,precise_type>(rc.alpha, rc.sqrbeta);
448 inline std::pair<precise_type,precise_type>
451 const Recurrence& rc(rCoeffs_.at(k));
452 return std::pair<precise_type,precise_type>(rc.alpha, rc.beta);
456 static const precise_type precise_zero;
460 inline Recurrence( const precise_type a,
461 const precise_type b)
465 sqrbeta = std::sqrt(beta);
470 precise_type sqrbeta;
474 friend class PolyFcn;
476 class PolyFcn : public Functor1<precise_type, precise_type>
480 const unsigned d1, const unsigned power)
481 : poly(p), deg1(d1), polypow(power) {}
483 inline virtual ~PolyFcn() {}
485 inline precise_type operator()( const precise_type& x) const
487 precise_type p = 1.0;
490 p = poly.normpoly(deg1, x);
499 const precise_type myPow = polypow;
500 p = std::pow(p, myPow);
514 friend class TripleProdFcn;
516 class TripleProdFcn : public Functor1<precise_type, precise_type>
520 const unsigned d2, const unsigned d3)
526 degmax = std::max(d1, std::max(d2, d3));
529 inline virtual ~TripleProdFcn() {}
531 inline precise_type operator()( const precise_type& x) const
532 { return poly.normpolyprod(degs, 3, degmax, x);}
540 template< class Funct> class EWPolyProductFcn;
541 template< class Funct> friend class EWPolyProductFcn;
543 template< class Funct>
544 class EWPolyProductFcn : public Functor1<precise_type, precise_type>
548 const unsigned d1, const unsigned d2)
553 degmax = std::max(d1, d2);
556 inline virtual ~EWPolyProductFcn() {}
558 inline precise_type operator()( const precise_type& x) const
559 { return poly.normpolyprod(degs,2,degmax,x- poly.shiftX_)*wf(x);}
571 inline MemoKey( const unsigned d0, const unsigned d1,
572 const unsigned d2, const unsigned d3)
579 std::sort(degs_, degs_+nDegs_);
582 inline MemoKey( const unsigned d0, const unsigned d1,
583 const unsigned d2, const unsigned d3,
584 const unsigned d4, const unsigned d5)
593 std::sort(degs_, degs_+nDegs_);
596 inline MemoKey( const unsigned d0, const unsigned d1,
597 const unsigned d2, const unsigned d3,
598 const unsigned d4, const unsigned d5,
599 const unsigned d6, const unsigned d7)
610 std::sort(degs_, degs_+nDegs_);
613 inline const unsigned* degrees() const { return degs_;}
614 inline unsigned nDegrees() const { return nDegs_;}
616 inline bool operator<( const MemoKey& r) const
618 if (nDegs_ < r.nDegs_)
620 if (r.nDegs_ < nDegs_)
622 for ( unsigned i=0; i<nDegs_; ++i)
624 if (degs_[i] < r.degs_[i])
626 if (r.degs_[i] < degs_[i])
632 inline bool operator==( const MemoKey& r) const
634 if (nDegs_ != r.nDegs_)
636 for ( unsigned i=0; i<nDegs_; ++i)
637 if (degs_[i] != r.degs_[i])
642 inline bool operator!=( const MemoKey& r) const
643 { return !(* this == r);}
650 template< typename Numeric1, typename Numeric2>
651 void copyMeasure( const std::vector<std::pair<Numeric1,Numeric2> >& inMeasure);
653 template< typename Numeric>
654 void copyInputCoords( const std::vector<Numeric>& inCoords);
656 void calcMeanXandWeightSums( bool shiftTheMeasure);
659 void calcRecurrenceStieltjes();
660 void calcRecurrenceLanczos();
662 template< class OrthoPoly>
663 long double sigma( const OrthoPoly& ops, int k, int l,
664 Matrix<long double>& sigMatrix,
665 Matrix<int>& flagMatrix) const;
667 template< class OrthoPoly>
668 void calculateMonicAverages( const OrthoPoly& ops, long double* averages,
671 template< class OrthoPoly>
672 void calcRecurrenceModifiedChebyshev( const OrthoPoly& ops);
674 precise_type monicpoly( unsigned degree, precise_type x) const;
676 std::pair<precise_type,precise_type> monicInnerProducts(
677 unsigned degree) const;
679 precise_type normpoly( unsigned degree, precise_type x) const;
681 std::pair<precise_type,precise_type> twonormpoly(
682 unsigned deg1, unsigned deg2, precise_type x) const;
684 precise_type normpolyprod( const unsigned* degrees, unsigned nDeg,
685 unsigned degmax, precise_type x) const;
687 template < typename Numeric>
688 precise_type normseries( const Numeric* coeffs, unsigned maxdeg,
689 precise_type x) const;
691 template < typename Numeric>
692 void getAllPolys( double x, unsigned maxdeg, Numeric *polyValues) const;
694 double getCachedAverage( const MemoKey& key) const;
696 std::vector<MeasurePoint> measure_;
697 std::vector<Recurrence> rCoeffs_;
698 mutable std::map<MemoKey,double> cachedAverages_;
700 precise_type wsumsq_;
702 precise_type shiftX_;
706 bool allWeightsEqual_;
710 inline StorablePolySeries1D* makeStorablePolySeries_2(
711 double xmin, double xmax, const std::vector<double>& coeffs) const
713 CPP11_auto_ptr<StorablePolySeries1D> ptr;
719 return ptr.release();
722 inline double sumOfWeights_2() const
723 { return sumOfWeights();}
725 inline double sumOfWeightSquares_2() const
726 { return sumOfWeightSquares();}
728 inline double meanCoordinate_2() const
729 { return meanCoordinate();}
731 inline double empiricalKroneckerDelta_2( const unsigned deg1,
732 const unsigned deg2) const
735 inline Matrix<double> jacobiMatrix_2( const unsigned n) const
738 inline double integratePoly_2( const unsigned degree, const unsigned power,
739 const double xmin, const double xmax) const
742 inline double integrateTripleProduct_2( const unsigned deg1, const unsigned deg2,
743 const unsigned deg3, const double xmin,
744 const double xmax) const
747 inline Matrix<double> tripleProductMatrix_2(
748 const std::vector<std::pair<unsigned,unsigned> >& pairs,
749 const unsigned maxdeg, const double xmin, const double xmax) const
752 inline double location2() const { return location();}
754 inline double scale2() const { return scale();}
756 inline std::pair<double,double> monicRecurrenceCoeffs2( const unsigned deg) const
759 return std::pair<double,double>(p.first, p.second);
762 inline std::pair<double,double> recurrenceCoeffs2( const unsigned deg) const
765 return std::pair<double,double>(p.first, p.second);
780 const unsigned degree,
781 const long double normfactor=1.0L)
782 : fcn_(fcn), norm_(normfactor), deg_(degree)
784 if (deg_ > fcn_. maxDegree()) throw std::invalid_argument(
785 "In npstat::ContOrthoPoly1DDeg::constructor: "
786 "degree argument is out of range");
791 inline virtual long double operator()( const long double& a) const
792 { return norm_*fcn_. poly(deg_, a);}
803 #include "npstat/nm/ContOrthoPoly1D.icc"
Base class for quadratures on the [-1, 1] interval.
bool operator==(const npstat::BoxND< Numeric > &l, const npstat::BoxND< Numeric > &r)
Gauss-Legendre quadratures in quad precision.
Gauss-Legendre quadratures in long double precision.
Enumeration of methods used to create orthogonal polynomials with discrete weights.
Compile-time deduction of an appropriate precise numeric type.
Interface definitions and concrete simple functors for a variety of functor-based calculations.
Definition: AbsIntervalQuadrature1D.hh:19
Definition: ContOrthoPoly1D.hh:777
Definition: ContOrthoPoly1D.hh:40
CPP11_auto_ptr< StorablePolySeries1D > makeStorablePolySeries(double xmin, double xmax, const double *coeffs=0, unsigned maxdeg=0) const
MeasurePoint measure(const unsigned long index) const Definition: ContOrthoPoly1D.hh:125
ContOrthoPoly1D(unsigned maxDegree, const std::vector< Numeric > &coords, OrthoPolyMethod m=OPOLY_STIELTJES, bool shiftMeasureCoords=true)
ContOrthoPoly1D(unsigned maxDegree, const std::vector< Numeric > &coords, const OrthoPoly &ops, bool shiftMeasureCoords=false)
double integrateEWPolyProduct(const Functor &weight, unsigned deg1, unsigned deg2, double xmin, double xmax, unsigned integrationRulePoints) const
void sampleAverages(const Numeric *coords, unsigned long lenCoords, double *averages, unsigned maxdeg) const
double empiricalKroneckerCovariance(unsigned deg1, unsigned deg2, unsigned deg3, unsigned deg4) const
double cov8(unsigned a, unsigned b, unsigned c, unsigned d, unsigned e, unsigned f, unsigned g, unsigned h) const
Matrix< double > epsCovarianceMatrix(const std::vector< std::pair< unsigned, unsigned > > &pairs, bool highOrder) const
double epsExpectation(unsigned m, unsigned n, bool highOrder) const
double integratedSquaredError(const Functor &fcn, const double *coeffs, unsigned maxdeg, double xmin, double xmax, unsigned integrationRulePoints=0U) const
void weightedPointsAverages(const Pair *points, unsigned long numPoints, double *averages, unsigned maxdeg) const
double slowCov8(unsigned a, unsigned b, unsigned c, unsigned d, unsigned e, unsigned f, unsigned g, unsigned h) const
ContOrthoPoly1D(unsigned maxDegree, const std::vector< std::pair< Numeric1, Numeric2 > > &measure, const OrthoPoly &ops, bool shiftMeasureCoords=false)
Matrix< double > sampleProductAverages(const Numeric *coords, unsigned long lenCoords, unsigned maxdeg) const
void extWeightAverages(const Functor &weight, const AbsIntervalQuadrature1D &quad, unsigned maxdeg, double xmin, double xmax, double *averages) const
double covCov4(unsigned a, unsigned b, unsigned c, unsigned d, unsigned e, unsigned f, unsigned g, unsigned h) const
double constrainedCoeffs(const Functor &fcn, double *coeffs, unsigned maxdeg, double xmin, double xmax, double integralValue=1.0) const
ContOrthoPoly1D(unsigned maxDegree, const std::vector< std::pair< Numeric1, Numeric2 > > &measure, OrthoPolyMethod m=OPOLY_STIELTJES, bool shiftMeasureCoords=true)
void allPolys(const double x, const unsigned maxdeg, Numeric *polyValues) const Definition: ContOrthoPoly1D.hh:147
precise_type integratePoly(unsigned degree, unsigned power, double xmin, double xmax) const
double powerAverage(unsigned deg, unsigned power) const
double series(const Numeric *coeffs, unsigned deg, double x) const
std::pair< precise_type, precise_type > monicRecurrenceCoeffs(const unsigned k) const Definition: ContOrthoPoly1D.hh:449
void weightCoeffs(double *coeffs, unsigned maxdeg) const
void allpoly(const long double x, long double *values, const unsigned maxdeg) const Definition: ContOrthoPoly1D.hh:155
double cov4(unsigned a, unsigned b, unsigned c, unsigned d) const
Matrix< double > extWeightProductAverages(const Functor &weight, const AbsIntervalQuadrature1D &quad, unsigned maxdeg, double xmin, double xmax) const
std::pair< double, double > polyPair(unsigned deg1, unsigned deg2, double x) const
precise_type empiricalKroneckerDelta(unsigned deg1, unsigned deg2) const
double cachedJointAverage(unsigned deg1, unsigned deg2, unsigned deg3, unsigned deg4) const
Matrix< precise_type > tripleProductMatrix(const std::vector< std::pair< unsigned, unsigned > > &pairs, unsigned maxdeg, double xmin, double xmax) const
double poly(unsigned deg, double x) const
ContOrthoPoly1D * clone() const Definition: ContOrthoPoly1D.hh:101
std::pair< precise_type, precise_type > recurrenceCoeffs(const unsigned k) const Definition: ContOrthoPoly1D.hh:438
double weightedSquaredError(const Functor &fcn, const double *coeffs, unsigned maxdeg) const
double epsCovariance(unsigned i, unsigned j, unsigned m, unsigned n, bool highOrder) const
double cov6(unsigned a, unsigned b, unsigned c, unsigned d, unsigned e, unsigned f) const
unsigned maxDegree() const Definition: ContOrthoPoly1D.hh:105
precise_type integrateTripleProduct(unsigned deg1, unsigned deg2, unsigned deg3, double xmin, double xmax) const
void calculateCoeffs(const Functor &fcn, Numeric *coeffs, unsigned maxdeg) const
Matrix< precise_type > jacobiMatrix(unsigned n) const
double jointAverage(const unsigned *degrees, unsigned nDegrees, bool degreesSortedIncreasingOrder=false) const
double effectiveSampleSize() const
double jointPowerAverage(unsigned deg1, unsigned power1, unsigned deg2, unsigned power2) const
Matrix< precise_type > empiricalKroneckerMatrix(unsigned maxdeg) const
Definition: GaussLegendreQuadratureQ.hh:26
Definition: GaussLegendreQuadrature.hh:27
Definition: AbsArrayProjector.hh:14
OrthoPolyMethod Definition: OrthoPolyMethod.hh:20
Definition: SimpleFunctors.hh:58
|