Index: trunk/npstat/nm/ContOrthoPoly1D.icc =================================================================== --- trunk/npstat/nm/ContOrthoPoly1D.icc (revision 551) +++ trunk/npstat/nm/ContOrthoPoly1D.icc (revision 552) @@ -1,106 +1,121 @@ #include #include "npstat/nm/GaussLegendreQuadrature.hh" namespace npstat { template void ContOrthoPoly1D::calculateCoeffs(const Functor& fcn, double *coeffs, const unsigned maxdeg) const { if (maxdeg > maxdeg_) throw std::invalid_argument( "In npstat::ContOrthoPoly1D::calculateCoeffs: " "maximum degree out of range"); assert(coeffs); std::vector scalarProducts(maxdeg+1U, 0.0L); const unsigned long mSize = measure_.size(); for (unsigned long i = 0; i < mSize; ++i) { const double x = measure_[i].first; const double w = measure_[i].second; const long double f = fcn(x + meanX_); long double polyk = 1.0L, polykm1 = 0.0L; for (unsigned k=0; k double ContOrthoPoly1D::weightedSquaredError(const Functor& fcn, const double *coeffs, const unsigned maxdeg) const { if (maxdeg > maxdeg_) throw std::invalid_argument( "In npstat::ContOrthoPoly1D::weightedSquaredError: " "maximum degree out of range"); assert(coeffs); long double sum = 0.0L; const unsigned long mSize = measure_.size(); for (unsigned long i = 0; i < mSize; ++i) { const double x = measure_[i].first; const long double delta = normseries(coeffs, maxdeg, x) - fcn(x + meanX_); sum += delta*delta*measure_[i].second; } return sum/wsum_; } template double ContOrthoPoly1D::integratedSquaredError( const Functor& infcn, const double *coeffs, const unsigned maxdeg, const double xmin, const double xmax, unsigned integrationRulePoints) const { if (maxdeg > maxdeg_) throw std::invalid_argument( "In npstat::ContOrthoPoly1D::integratedSquaredError: " "maximum degree out of range"); assert(coeffs); if (!integrationRulePoints) { integrationRulePoints = GaussLegendreQuadrature::minimalExactRule(2*maxdeg); if (!integrationRulePoints) throw std::invalid_argument( "In npstat::ContOrthoPoly1D::integratedSquaredError: " "order of requested integration rule is too high"); } DeltaSquaredSerFcn fcn(*this, infcn, coeffs, maxdeg); GaussLegendreQuadrature glq(integrationRulePoints); return glq.integrate(fcn, xmin, xmax); } + template + double ContOrthoPoly1D::integrateEWPolyProduct( + const Functor& weight, + const unsigned deg1, const unsigned deg2, + const double xmin, const double xmax, + const unsigned integrationRulePoints) const + { + if (deg1 > maxdeg_ || deg2 > maxdeg_) throw std::invalid_argument( + "In npstat::ContOrthoPoly1D::integrateEWPolyProduct: " + "degree argument out of range"); + GaussLegendreQuadrature glq(integrationRulePoints); + EWPolyProductFcn fcn(*this, weight, deg1, deg2); + return glq.integrate(fcn, xmin, xmax); + } + template void ContOrthoPoly1D::getAllPolys(const double x, const unsigned maxdeg, Numeric *polyValues) const { if (maxdeg > maxdeg_) throw std::invalid_argument( "In npstat::ContOrthoPoly1D::getAllPolys: " "maximum degree out of range"); assert(polyValues); long double polyk = 1.0L, polykm1 = 0.0L; polyValues[0] = polyk; for (unsigned k=0; k #include #include #include #include #include "npstat/nm/OrthoPolyMethod.hh" #include "npstat/nm/Matrix.hh" #include "npstat/nm/SimpleFunctors.hh" namespace npstat { class ContOrthoPoly1D { public: // The first element of the pair is the coordinate and the // second measure is the weight (all weights must be non-negative) typedef std::pair MeasurePoint; /** // Main constructor, with obvious arguments. Internally, the // measure will be sorted in the order of increasing weight // and the measure coordinates will be shifted so that their // weighted mean is at 0. */ ContOrthoPoly1D(unsigned maxDegree, const std::vector& measure, OrthoPolyMethod m = OPOLY_STIELTJES); /** Constructor which assumes that all weights are 1.0 */ ContOrthoPoly1D(unsigned maxDegree, const std::vector& coords, OrthoPolyMethod m = OPOLY_STIELTJES); //@{ /** A simple inspector of object properties */ inline unsigned maxDegree() const {return maxdeg_;} inline unsigned long measureLength() const {return measure_.size();} inline MeasurePoint measure(const unsigned index) const {return measure_.at(index);} inline double coordinate(const unsigned index) const {return measure_.at(index).first;} inline double weight(const unsigned index) const {return measure_.at(index).second;} inline double sumOfWeights() const {return wsum_;} inline double sumOfWeightSquares() const {return wsumsq_;} inline bool areAllWeightsEqual() const {return allWeightsEqual_;} inline double minCoordinate() const {return minX_;} inline double maxCoordinate() const {return maxX_;} inline double meanCoordinate() const {return meanX_;} //@} /** Kish's effective sample size for the measure */ double effectiveSampleSize() const; /** Return the value of one of the normalized polynomials */ double poly(unsigned deg, double x) const; /** Return the values of two of the normalized polynomials */ std::pair polyPair( unsigned deg1, unsigned deg2, double 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. */ void allPolys(double x, unsigned maxdeg, double *polyValues) const; void allPolys(double x, unsigned maxdeg, long double *polyValues) const; //@} /** Return the value of the orthonormal polynomial series */ double series(const double *coeffs, unsigned maxdeg, double x) const; /** // 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, double *coeffs, unsigned maxdeg) const; /** // Integrated squared error between the given function and the // polynomial series. Parameter "integrationRulePoints" specifies // the parameter "npoints" for the GaussLegendreQuadrature class. // If "integrationRulePoints" is 0, the rule will be picked // automatically in such a way that it integrates a polynomial // of degree maxdeg*2 exactly. */ template double integratedSquaredError(const Functor& fcn, const double *coeffs, unsigned maxdeg, double xmin, double xmax, unsigned integrationRulePoints=0U) const; /** // Squared error between the given function and the polynomial // series, weighted according to the measure */ template double weightedSquaredError(const Functor& fcn, const double *coeffs, unsigned maxdeg) const; /** // This method is useful for testing the numerical precision // of the orthonormalization procedure. It returns the scalar // products between various polynomials. */ double empiricalKroneckerDelta(unsigned deg1, unsigned deg2) const; /** // If the Kronecker deltas were calculated from a sample, they // would be random and would have a covariance matrix. This // is an estimate of the terms in that covariance matrix. The // covariance is between delta_{deg1,deg2} and delta_{deg3,deg4}. */ double empiricalKroneckerCovariance(unsigned deg1, unsigned deg2, unsigned deg3, unsigned deg4) const; /** // Examine the recurrence coefficients. The first element of // the returned pair is alpha, and the second element is beta. */ std::pair recurrenceCoeffs(unsigned deg) const; /** // Generate principal minor of order n of the Jacobi matrix. // n must not exceed "maxDegree" */ Matrix jacobiMatrix(unsigned n) const; /** // Roots of the polynomial with the given degree. Naturally, // the degree argument must not exceed "maxDegree". */ void calculateRoots(double *roots, unsigned degree) const; /** // Integral of an orthonormal polynomials to some power without // the weight (so this is not an inner product with the poly of // degree 0). */ double integratePoly(unsigned degree, unsigned power, double xmin, double xmax) const; /** // Integral of the product of three orthonormal polynomials // without the weight (so this is not an inner product) */ double integrateTripleProduct(unsigned deg1, unsigned deg2, unsigned deg3, double xmin, double xmax) const; /** + // Integral of the product of two orthonormal polynomials + // with some external weight function (e.g., oracle density). + // "EW" in the method name stands for "externally weighted". + // + // "integrationRulePoints" argument will be passed to the + // GaussLegendreQuadrature constructor and must be supported + // by that class. + */ + template + double integrateEWPolyProduct(const Functor& weight, + unsigned deg1, unsigned deg2, + double xmin, double xmax, + unsigned integrationRulePoints) const; + + /** // A measure-weighted average of orthonormal poly values // to some power */ double powerAverage(unsigned deg, unsigned power) const; /** // A measure-weighted average of the product of two orthonormal // poly values to some powers */ double jointPowerAverage(unsigned deg1, unsigned power1, unsigned deg2, unsigned power2) const; /** // A measure-weighted average of the product of several // orthonormal poly values */ double jointAverage(const unsigned* degrees, unsigned nDegrees, bool degreesSortedIncreasingOrder = false) const; //@{ /** // A measure-weighted average that will be remembered internally // so that another call to this function with the same arguments // will return the answer quickly */ double cachedJointAverage(unsigned deg1, unsigned deg2, unsigned deg3, unsigned deg4) const; double cachedJointAverage(unsigned deg1, unsigned deg2, unsigned deg3, unsigned deg4, unsigned deg5, unsigned deg6) const; double cachedJointAverage(unsigned deg1, unsigned deg2, unsigned deg3, unsigned deg4, unsigned deg5, unsigned deg6, unsigned deg7, unsigned deg8) const; //@} /** Covariance between v_{ab} and v_{cd} */ double cov4(unsigned a, unsigned b, unsigned c, unsigned d) const; /** Covariance between v_{ab}*v_{cd} and v_{ef} */ double cov6(unsigned a, unsigned b, unsigned c, unsigned d, unsigned e, unsigned f) const; /** Covariance between v_{ab}*v_{cd} and v_{ef}*v_{gh} */ double cov8(unsigned a, unsigned b, unsigned c, unsigned d, unsigned e, unsigned f, unsigned g, unsigned h) const; /** // Alternative implementation of "cov8". It is slower than "cov8" // but easier to program and verify. Useful mainly for testing "cov8". */ double slowCov8(unsigned a, unsigned b, unsigned c, unsigned d, unsigned e, unsigned f, unsigned g, unsigned h) const; double epsExpectation(unsigned m, unsigned n, bool highOrder) const; double epsCovariance(unsigned m1, unsigned n1, unsigned m2, unsigned n2, bool highOrder) const; private: struct Recurrence { inline Recurrence(const long double a, const long double b) : alpha(a), beta(b) { assert(beta > 0.0L); sqrbeta = sqrtl(beta); } long double alpha; long double beta; long double sqrbeta; }; class PolyFcn; friend class PolyFcn; class PolyFcn : public Functor1 { public: inline PolyFcn(const ContOrthoPoly1D& p, const unsigned d1, const unsigned power) : poly(p), deg1(d1), polypow(power) {} inline long double operator()(const long double& x) const { long double p = 1.0L; if (polypow) { p = poly.normpoly(deg1, x); switch (polypow) { case 1U: break; case 2U: p *= p; break; default: p = powl(p, polypow); break; } } return p; } private: const ContOrthoPoly1D& poly; unsigned deg1; unsigned polypow; }; class TripleProdFcn; friend class TripleProdFcn; class TripleProdFcn : public Functor1 { public: inline TripleProdFcn(const ContOrthoPoly1D& p, const unsigned d1, const unsigned d2, const unsigned d3) : poly(p) { degs[0] = d1; degs[1] = d2; degs[2] = d3; degmax = std::max(d1, std::max(d2, d3)); } inline long double operator()(const long double& x) const {return poly.normpolyprod(degs, 3, degmax, x);} private: const ContOrthoPoly1D& poly; unsigned degs[3]; unsigned degmax; }; + template class EWPolyProductFcn; + template friend class EWPolyProductFcn; + + template + class EWPolyProductFcn : public Functor1 + { + public: + inline EWPolyProductFcn(const ContOrthoPoly1D& p, const Funct& w, + const unsigned d1, const unsigned d2) + : poly(p), wfun(w) + { + degs[0] = d1; + degs[1] = d2; + degmax = std::max(d1, d2); + } + + inline long double operator()(const long double& x) const + {return poly.normpolyprod(degs, 2, degmax, x)*wfun(x);} + + private: + const ContOrthoPoly1D& poly; + const Funct& wfun; + unsigned degs[2]; + unsigned degmax; + }; + class MemoKey { public: inline MemoKey(const unsigned d0, const unsigned d1, const unsigned d2, const unsigned d3) : nDegs_(4U) { degs_[0] = d0; degs_[1] = d1; degs_[2] = d2; degs_[3] = d3; std::sort(degs_, degs_+nDegs_); } inline MemoKey(const unsigned d0, const unsigned d1, const unsigned d2, const unsigned d3, const unsigned d4, const unsigned d5) : nDegs_(6U) { degs_[0] = d0; degs_[1] = d1; degs_[2] = d2; degs_[3] = d3; degs_[4] = d4; degs_[5] = d5; std::sort(degs_, degs_+nDegs_); } inline MemoKey(const unsigned d0, const unsigned d1, const unsigned d2, const unsigned d3, const unsigned d4, const unsigned d5, const unsigned d6, const unsigned d7) : nDegs_(8U) { degs_[0] = d0; degs_[1] = d1; degs_[2] = d2; degs_[3] = d3; degs_[4] = d4; degs_[5] = d5; degs_[6] = d6; degs_[7] = d7; std::sort(degs_, degs_+nDegs_); } inline const unsigned* degrees() const {return degs_;} inline unsigned nDegrees() const {return nDegs_;} inline bool operator<(const MemoKey& r) const { if (nDegs_ < r.nDegs_) return true; if (r.nDegs_ < nDegs_) return false; for (unsigned i=0; i monicInnerProducts( unsigned degree) const; long double normpoly(unsigned degree, long double x) const; std::pair twonormpoly( unsigned deg1, unsigned deg2, long double x) const; long double normpolyprod(const unsigned* degrees, unsigned nDeg, unsigned degmax, long double x) const; long double normseries(const double *coeffs, unsigned maxdeg, long double x) const; template void getAllPolys(double x, unsigned maxdeg, Numeric *polyValues) const; double getCachedAverage(const MemoKey& key) const; std::vector measure_; std::vector rCoeffs_; mutable std::map cachedAverages_; long double wsum_; long double wsumsq_; double meanX_; double minX_; double maxX_; unsigned maxdeg_; bool allWeightsEqual_; }; } #include "npstat/nm/ContOrthoPoly1D.icc" #endif // NPSTAT_CONTORTHOPOLY1D_HH_ Index: trunk/npstat/stat/orthoPoly1DVProducts.hh =================================================================== --- trunk/npstat/stat/orthoPoly1DVProducts.hh (revision 551) +++ trunk/npstat/stat/orthoPoly1DVProducts.hh (revision 552) @@ -1,106 +1,125 @@ #ifndef NPSTAT_ORTHOPOLY1DVPRODUCTS_HH_ #define NPSTAT_ORTHOPOLY1DVPRODUCTS_HH_ /*! // \file orthoPoly1DVProducts.hh // // \brief Utility functions for calculating statistical properties of // 1-d orthogonal polynomials // // Author: I. Volobouev // // September 2017 */ #include "npstat/nm/ArrayND.hh" #include "npstat/nm/AbsClassicalOrthoPoly1D.hh" namespace npstat { typedef std::pair UUPair; + // See the "Empirical Chaos Polynomials" write-up for the + // definition of v_{ab} and various related formulae for + // expectations and covariances template double expectedVProduct(const AbsClassicalOrthoPoly1D& poly, const Quadrature& quad, unsigned long nPoints, UUPair ab, UUPair cd); template double expectedVProduct(const AbsClassicalOrthoPoly1D& poly, const Quadrature& quad, unsigned long nPoints, UUPair ab, UUPair cd, UUPair ef); template double expectedVProduct(const AbsClassicalOrthoPoly1D& poly, const Quadrature& quad, unsigned long nPoints, UUPair ab, UUPair cd, UUPair ef, UUPair gh); + // Covariance between v_{ab} and v_{cd} template double expectedVCovariance(const AbsClassicalOrthoPoly1D& poly, const Quadrature& quad, unsigned long nPoints, UUPair ab, UUPair cd); + // Covariance between v_{ab}*v_{cd} and v_{ef} template double expectedVCovariance(const AbsClassicalOrthoPoly1D& poly, const Quadrature& quad, unsigned long nPoints, UUPair ab, UUPair cd, UUPair ef); + // Covariance between v_{ab}*v_{cd} and v_{ef}*v_{gh} template double expectedVCovariance(const AbsClassicalOrthoPoly1D& poly, const Quadrature& quad, unsigned long nPoints, UUPair ab, UUPair cd, UUPair ef, UUPair gh); + // v_{ab} calculated for a given sample of unweighted points. + // For samples generated using the polynomial weight function + // as the density, the expectation value of v_{ab} is 0. template double sampleVProduct(const AbsClassicalOrthoPoly1D& poly, const Numeric* coords, unsigned long nPoints, UUPair ab); + // v_{ab}*v_{cd} calculated for a given sample of unweighted points template double sampleVProduct(const AbsClassicalOrthoPoly1D& poly, const Numeric* coords, unsigned long nPoints, UUPair ab, UUPair cd); + // v_{ab}*v_{cd}*v_{ef} calculated for a given sample template double sampleVProduct(const AbsClassicalOrthoPoly1D& poly, const Numeric* coords, unsigned long nPoints, UUPair ab, UUPair cd, UUPair ef); + // v_{ab}*v_{cd}*v_{ef}*v_{gh} calculated for a given sample template double sampleVProduct(const AbsClassicalOrthoPoly1D& poly, const Numeric* coords, unsigned long nPoints, UUPair ab, UUPair cd, UUPair ef, UUPair gh); // Oracle counterpart of the "epsExpectation" method of the // ContOrthoPoly1D class. If "highOrder" parameter is "true", // the calculations will be performed to O(N^{-3/2}). If "highOrder" // is "false", the calculations will be to O(N^{-1}). template double oracleEpsExpectation(const AbsClassicalOrthoPoly1D& poly, const Quadrature& quad, unsigned long nPoints, unsigned m, unsigned n, bool highOrder); // Oracle counterpart of the "epsCovariance" method of the // ContOrthoPoly1D class template double oracleEpsCovariance(const AbsClassicalOrthoPoly1D& poly, const Quadrature& quad, unsigned long nPoints, unsigned m1, unsigned n1, unsigned m2, unsigned n2, bool highOrder); // Fill a 4-d array with oracle covariances. The code knows that the // covariances are symmetric with respect to permuting all indices, and // does not recalculate them when permuting indices is sufficient. The // array span will be maxdeg+1 in each dimension. template ArrayND oracleEpsCovarianceArray( const AbsClassicalOrthoPoly1D& poly, const Quadrature& quad, unsigned long nPoints, unsigned maxdeg, bool highOrder); - // Calculate the "eps" value for a given sample of (random) values + // Calculate the "eps_{mn}" approximation for a given sample + // of unweighted points (up to O(N^{-1}) or O(N^{-3/2})). In + // this approximation, all v_{ij} values are calculated using + // sample averages of oracle polynomials. Note that the result + // is not the exact "eps_{mn}" for the sample but just this + // particular special approximation. It is useful for checking + // the effect of replacing the empirical polys with the oracle + // ones in the "eps_{mn}" approximation. template double sampleEpsValue(const AbsClassicalOrthoPoly1D& poly, const Numeric* coords, unsigned long nPoints, unsigned m, unsigned n, bool highOrder); } #include "npstat/stat/orthoPoly1DVProducts.icc" #endif // NPSTAT_ORTHOPOLY1DVPRODUCTS_HH_