Index: trunk/npstat/nm/ContOrthoPoly1D.hh =================================================================== --- trunk/npstat/nm/ContOrthoPoly1D.hh (revision 557) +++ trunk/npstat/nm/ContOrthoPoly1D.hh (revision 558) @@ -1,483 +1,491 @@ #ifndef NPSTAT_CONTORTHOPOLY1D_HH_ #define NPSTAT_CONTORTHOPOLY1D_HH_ /*! // \file ContOrthoPoly1D.hh // // \brief Continuous orthogonal polynomial series in one dimension // generated for a discrete measure // // Author: I. Volobouev // // May 2017 */ #include #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; /** + // Build the coefficients of the orthogonal polynomial series + // for the weight function itself (that is, fcn(x_i) = w_i, + // using the terminology of the "calculateCoeffs" method). + // This can sometimes be useful for weighted measures. + */ + void weightsSquaredCoeffs(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; /** Estimate of eps_{mn} expectation */ double epsExpectation(unsigned m, unsigned n, bool highOrder) const; /** // Estimate of covariance between eps_{ij} and eps_{mn}. // // The result should be identical with "empiricalKroneckerCovariance" // in case "highOrder" argument is "false". However, this method // caches results of various calculations of averages and should // perform better inside cycles over indices. */ double epsCovariance(unsigned i, unsigned j, unsigned m, unsigned n, 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), wf(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-poly.meanX_)*wf(x);} private: const ContOrthoPoly1D& poly; const Funct& wf; 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/nm/ContOrthoPoly1D.cc =================================================================== --- trunk/npstat/nm/ContOrthoPoly1D.cc (revision 557) +++ trunk/npstat/nm/ContOrthoPoly1D.cc (revision 558) @@ -1,880 +1,913 @@ #include #include #include #include "npstat/nm/ContOrthoPoly1D.hh" #include "npstat/nm/PairCompare.hh" #include "npstat/nm/allocators.hh" inline static int kdelta(const unsigned i, const unsigned j) { return i == j ? 1 : 0; } namespace npstat { ContOrthoPoly1D::ContOrthoPoly1D(const unsigned maxDegree, const std::vector& inCoords, const OrthoPolyMethod m) : wsum_(inCoords.size()), wsumsq_(inCoords.size()), minX_(DBL_MAX), maxX_(-DBL_MAX), maxdeg_(maxDegree), allWeightsEqual_(true) { if (inCoords.empty()) throw std::invalid_argument( "In npstat::ContOrthoPoly1D constructor: empty list of coordinates"); const unsigned long sz = inCoords.size(); if (sz <= maxDegree) throw std::invalid_argument( "In npstat::ContOrthoPoly1D constructor: insufficient number of points"); long double xsum = std::accumulate(inCoords.begin(), inCoords.end(), 0.0L); meanX_ = xsum/sz; // Improve the numerical precision of the mean x subtraction xsum = 0.0L; for (unsigned long i = 0; i < sz; ++i) { const double x = inCoords[i]; if (x < minX_) minX_ = x; if (x > maxX_) maxX_ = x; xsum += (x - meanX_); } const double dx = xsum/sz; // Fill the measure measure_.reserve(sz); for (unsigned long i = 0; i < sz; ++i) measure_.push_back(MeasurePoint(inCoords[i] - meanX_ - dx, 1.0)); meanX_ += dx; calcRecurrenceCoeffs(m); } ContOrthoPoly1D::ContOrthoPoly1D(const unsigned maxDegree, const std::vector& inMeasure, const OrthoPolyMethod m) : measure_(inMeasure), wsum_(0.0L), wsumsq_(0.0L), minX_(DBL_MAX), maxX_(-DBL_MAX), maxdeg_(maxDegree), allWeightsEqual_(true) { // Check argument validity if (measure_.empty()) throw std::invalid_argument( "In npstat::ContOrthoPoly1D constructor: empty measure"); // We expect all weights to be equal quite often. // Check if this is indeed the case. If not, sort the // weights in the increasing order, hopefully reducing // the round-off error. const unsigned long mSize = measure_.size(); const double w0 = measure_[0].second; for (unsigned long i = 1; i < mSize && allWeightsEqual_; ++i) allWeightsEqual_ = (w0 == measure_[i].second); if (!allWeightsEqual_) std::sort(measure_.begin(), measure_.end(), LessBySecond()); if (measure_[0].second < 0.0) throw std::invalid_argument( "In npstat::ContOrthoPoly1D constructor: negative measure entry found"); unsigned long nZeroWeights = 0; if (!allWeightsEqual_) for (unsigned long i = 0; i < mSize; ++i) { if (measure_[i].second == 0.0) ++nZeroWeights; else break; } if (mSize <= maxDegree + nZeroWeights) throw std::invalid_argument( "In npstat::ContOrthoPoly1D constructor: insufficient number of measure points"); // Sum up the weights long double xsum = 0.0L; if (allWeightsEqual_) { const long double lw0 = w0; wsum_ = mSize*lw0; wsumsq_ = mSize*lw0*lw0; for (unsigned long i = 0; i < mSize; ++i) { const double x = measure_[i].first; xsum += x; if (x < minX_) minX_ = x; if (x > maxX_) maxX_ = x; } xsum *= lw0; } else { for (unsigned long i = 0; i < mSize; ++i) { const double x = measure_[i].first; const double w = measure_[i].second; wsum_ += w; wsumsq_ += w*w; xsum += w*x; if (x < minX_) minX_ = x; if (x > maxX_) maxX_ = x; } } if (wsum_ == 0.0L) throw std::invalid_argument( "In npstat::ContOrthoPoly1D constructor: measure is not positive"); meanX_ = xsum/wsum_; // Shift the measure for (unsigned long i = 0; i < mSize; ++i) measure_[i].first -= meanX_; // Try to improve the mean xsum = 0.0L; for (unsigned long i = 0; i < mSize; ++i) xsum += measure_[i].first*measure_[i].second; const double dx = xsum/wsum_; for (unsigned long i = 0; i < mSize; ++i) measure_[i].first -= dx; meanX_ += dx; calcRecurrenceCoeffs(m); } void ContOrthoPoly1D::calcRecurrenceCoeffs(const OrthoPolyMethod m) { rCoeffs_.clear(); rCoeffs_.reserve(maxdeg_ + 1); switch (m) { case OPOLY_STIELTJES: calcRecurrenceStieltjes(); break; case OPOLY_LANCZOS: calcRecurrenceLanczos(); break; default: throw std::runtime_error( "In npstat::ContOrthoPoly1D::calcRecurrenceCoeffs: " "incomplete switch statement. This is a bug. Please report."); } assert(rCoeffs_.size() == maxdeg_ + 1); } double ContOrthoPoly1D::effectiveSampleSize() const { if (allWeightsEqual_) return measure_.size(); else return wsum_*wsum_/wsumsq_; } long double ContOrthoPoly1D::monicpoly(const unsigned degree, const long double x) const { long double polyk = 1.0L, polykm1 = 0.0L; for (unsigned k=0; k maxdeg_) throw std::invalid_argument( "In npstat::ContOrthoPoly1D::powerAverage: " "degree argument is out of range"); const unsigned long mSize = measure_.size(); long double sum = 0.0L; if (power == 1U) for (unsigned long i = 0; i < mSize; ++i) sum += measure_[i].second*normpoly(deg, measure_[i].first); else for (unsigned long i = 0; i < mSize; ++i) sum += measure_[i].second*powl(normpoly(deg, measure_[i].first), power); return sum/wsum_; } std::pair ContOrthoPoly1D::twonormpoly( const unsigned deg1, const unsigned deg2, const long double x) const { long double p1 = 0.0L, p2 = 0.0L, polyk = 1.0L, polykm1 = 0.0L; const unsigned degmax = std::max(deg1, deg2); for (unsigned k=0; k(p1, p2); } long double ContOrthoPoly1D::normpolyprod(const unsigned* degrees, const unsigned nDeg, const unsigned degmax, const long double x) const { long double prod = 1.0L, polyk = 1.0L, polykm1 = 0.0L; for (unsigned k=0; k ContOrthoPoly1D::monicInnerProducts(const unsigned degree) const { if (degree) { long double sum = 0.0L, xsum = 0.0L; const unsigned long mSize = measure_.size(); for (unsigned long i = 0; i < mSize; ++i) { const long double x = measure_[i].first; const long double p = monicpoly(degree, x); const long double pprod = p*p*measure_[i].second; sum += pprod; xsum += x*pprod; } return std::pair(xsum/wsum_, sum/wsum_); } else // Average x should be 0 for degree == 0 return std::pair(0.0L, 1.0L); } + void ContOrthoPoly1D::weightsSquaredCoeffs(double *coeffs, + const unsigned maxdeg) const + { + if (maxdeg > maxdeg_) throw std::invalid_argument( + "In npstat::ContOrthoPoly1D::weightsSquaredCoeffs: " + "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 = w; + + long double polyk = 1.0L, polykm1 = 0.0L; + for (unsigned k=0; k ip = monicInnerProducts(deg); rCoeffs_.push_back(Recurrence(ip.first/ip.second, ip.second/prevSprod)); prevSprod = ip.second; } } void ContOrthoPoly1D::calcRecurrenceLanczos() { typedef long double Real; const unsigned long mSize = measure_.size(); std::vector dmem(mSize*2UL); Real* p0 = &dmem[0]; for (unsigned long i=0; i maxdeg_) throw std::invalid_argument( "In npstat::ContOrthoPoly1D::poly: degree argument is out of range"); return normpoly(deg, x - meanX_); } std::pair ContOrthoPoly1D::polyPair( const unsigned deg1, const unsigned deg2, const double x) const { if (deg1 > maxdeg_ || deg2 > maxdeg_) throw std::invalid_argument( "In npstat::ContOrthoPoly1D::polyPair: degree argument is out of range"); const std::pair& ld = twonormpoly(deg1, deg2, x - meanX_); return std::pair(ld.first, ld.second); } double ContOrthoPoly1D::series(const double *coeffs, const unsigned deg, const double x) const { assert(coeffs); if (deg > maxdeg_) throw std::invalid_argument( "In npstat::ContOrthoPoly1D::series: degree argument is out of range"); return normseries(coeffs, deg, x - meanX_); } void ContOrthoPoly1D::allPolys(const double x, const unsigned deg, double *polyValues) const { getAllPolys(x - meanX_, deg, polyValues); } void ContOrthoPoly1D::allPolys(const double x, const unsigned deg, long double *polyValues) const { getAllPolys(x - meanX_, deg, polyValues); } double ContOrthoPoly1D::empiricalKroneckerDelta( const unsigned ideg1, const unsigned ideg2) const { if (ideg1 > maxdeg_ || ideg2 > maxdeg_) throw std::invalid_argument( "In npstat::ContOrthoPoly1D::empiricalKroneckerDelta: " "degree argument is out of range"); long double sum = 0.0L; const unsigned long mSize = measure_.size(); for (unsigned long i = 0; i < mSize; ++i) { const std::pair& p = twonormpoly(ideg1, ideg2, measure_[i].first); sum += measure_[i].second*p.first*p.second; } return sum/wsum_; } double ContOrthoPoly1D::jointPowerAverage( const unsigned ideg1, const unsigned power1, const unsigned ideg2, const unsigned power2) const { if (power1 == 1U && power2 == 1U) return empiricalKroneckerDelta(ideg1, ideg2); else { if (ideg1 > maxdeg_ || ideg2 > maxdeg_) throw std::invalid_argument( "In npstat::ContOrthoPoly1D::jointPowerAverage: " "degree argument is out of range"); long double sum = 0.0L; const unsigned long mSize = measure_.size(); for (unsigned long i = 0; i < mSize; ++i) { const std::pair& p = twonormpoly(ideg1, ideg2, measure_[i].first); sum += measure_[i].second*powl(p.first,power1)*powl(p.second,power2); } return sum/wsum_; } } double ContOrthoPoly1D::jointAverage( const unsigned* degrees, const unsigned nDegrees, const bool sorted) const { if (nDegrees) { assert(degrees); unsigned degmax; if (sorted) degmax = degrees[nDegrees-1U]; else degmax = *std::max_element(degrees, degrees+nDegrees); if (degmax > maxdeg_) throw std::invalid_argument( "In npstat::ContOrthoPoly1D::jointAverage: " "degree argument is out of range"); long double sum = 0.0L; const unsigned long mSize = measure_.size(); for (unsigned long i = 0; i < mSize; ++i) sum += measure_[i].second*normpolyprod( degrees, nDegrees, degmax, measure_[i].first); return sum/wsum_; } else return 1.0; } double ContOrthoPoly1D::empiricalKroneckerCovariance( const unsigned deg1, const unsigned deg2, const unsigned deg3, const unsigned deg4) const { double cov = 0.0; if (!((deg1 == 0 && deg2 == 0) || (deg3 == 0 && deg4 == 0))) { unsigned degs[4]; degs[0] = deg1; degs[1] = deg2; degs[2] = deg3; degs[3] = deg4; cov = (jointAverage(degs, 4) - kdelta(deg1, deg2)*kdelta(deg3, deg4))/ effectiveSampleSize(); if (std::min(deg1, deg2) == std::min(deg3, deg4) && std::max(deg1, deg2) == std::max(deg3, deg4)) if (cov < 0.0) cov = 0.0; } return cov; } std::pair ContOrthoPoly1D::recurrenceCoeffs(const unsigned deg) const { if (deg > maxdeg_) throw std::invalid_argument( "In npstat::ContOrthoPoly1D::recurrenceCoeffs: " "degree argument is out of range"); const Recurrence& rc(rCoeffs_[deg]); return std::pair(rc.alpha, rc.beta); } Matrix ContOrthoPoly1D::jacobiMatrix(const unsigned n) const { if (!n) throw std::invalid_argument( "In npstat::ContOrthoPoly1D::jacobiMatrix: " "can not build matrix of zero size"); if (n > maxdeg_) throw std::invalid_argument( "In npstat::ContOrthoPoly1D::jacobiMatrix: " "matrix size is out of range"); Matrix mat(n, n, 0); unsigned ip1 = 1; for (unsigned i=0; i& mat = jacobiMatrix(deg); if (deg == 1U) roots[0] = mat[0][0]; else mat.tdSymEigen(roots, deg); for (unsigned i=0; i maxdeg_) throw std::invalid_argument( "In npstat::ContOrthoPoly1D::integratePoly: " "degree argument is out of range"); const unsigned nGood = GaussLegendreQuadrature::minimalExactRule(deg1*power); if (!nGood) throw std::invalid_argument( "In npstat::ContOrthoPoly1D::integratePoly: " "product of poly degree and power is too large"); GaussLegendreQuadrature glq(nGood); PolyFcn fcn(*this, deg1, power); return glq.integrate(fcn, xmin - meanX_, xmax - meanX_); } double ContOrthoPoly1D::integrateTripleProduct( const unsigned deg1, const unsigned deg2, const unsigned deg3, const double xmin, const double xmax) const { const unsigned maxdex = std::max(std::max(deg1, deg2), deg3); if (maxdex > maxdeg_) throw std::invalid_argument( "In npstat::ContOrthoPoly1D::integrateTripleProduct: " "degree argument is out of range"); const unsigned sumdeg = deg1 + deg2 + deg3; const unsigned nGood = GaussLegendreQuadrature::minimalExactRule(sumdeg); if (!nGood) throw std::invalid_argument( "In npstat::ContOrthoPoly1D::integrateTripleProduct: " "sum of the polynomial degrees is too large"); GaussLegendreQuadrature glq(nGood); TripleProdFcn fcn(*this, deg1, deg2, deg3); return glq.integrate(fcn, xmin - meanX_, xmax - meanX_); } double ContOrthoPoly1D::cachedJointAverage(const unsigned deg1, const unsigned deg2, const unsigned deg3, const unsigned deg4) const { MemoKey key(deg1, deg2, deg3, deg4); return getCachedAverage(key); } double ContOrthoPoly1D::cachedJointAverage(const unsigned deg1, const unsigned deg2, const unsigned deg3, const unsigned deg4, const unsigned deg5, const unsigned deg6) const { MemoKey key(deg1, deg2, deg3, deg4, deg5, deg6); return getCachedAverage(key); } double ContOrthoPoly1D::cachedJointAverage(const unsigned deg1, const unsigned deg2, const unsigned deg3, const unsigned deg4, const unsigned deg5, const unsigned deg6, const unsigned deg7, const unsigned deg8) const { MemoKey key(deg1, deg2, deg3, deg4, deg5, deg6, deg7, deg8); return getCachedAverage(key); } double ContOrthoPoly1D::getCachedAverage(const MemoKey& key) const { double value; std::map::const_iterator it = cachedAverages_.find(key); if (it == cachedAverages_.end()) { value = jointAverage(key.degrees(), key.nDegrees(), true); cachedAverages_.insert(std::pair(key, value)); } else value = it->second; return value; } double ContOrthoPoly1D::cov4(const unsigned deg1, const unsigned deg2, const unsigned deg3, const unsigned deg4) const { const bool has0 = (deg1 == 0 && deg2 == 0) || (deg3 == 0 && deg4 == 0); double cov = 0.0; if (!has0) { const double N = effectiveSampleSize(); cov = (cachedJointAverage(deg1, deg2, deg3, deg4) - kdelta(deg1, deg2)*kdelta(deg3, deg4))/N; if (std::min(deg1, deg2) == std::min(deg3, deg4) && std::max(deg1, deg2) == std::max(deg3, deg4)) if (cov < 0.0) cov = 0.0; } return cov; } double ContOrthoPoly1D::cov6(const unsigned a, const unsigned b, const unsigned c, const unsigned d, const unsigned e, const unsigned f) const { const bool has0 = (a == 0 && b == 0) || (c == 0 && d == 0) || (e == 0 && f == 0); double cov = 0.0; if (!has0) { double sum = cachedJointAverage(a, b, c, d, e, f); double add = 2.0; if (kdelta(a, b)) sum -= cachedJointAverage(c, d, e, f); else add = 0.0; if (kdelta(c, d)) sum -= cachedJointAverage(a, b, e, f); else add = 0.0; if (kdelta(e, f)) sum -= cachedJointAverage(a, b, c, d); else add = 0.0; const double N = effectiveSampleSize(); cov = (sum + add)/N/N; } return cov; } double ContOrthoPoly1D::slowCov8(const unsigned a, const unsigned b, const unsigned c, const unsigned d, const unsigned e, const unsigned f, const unsigned g, const unsigned h) const { const bool has0 = (a == 0 && b == 0) || (c == 0 && d == 0) || (e == 0 && f == 0) || (g == 0 && h == 0); double cov = 0.0; if (!has0) { const double pabcd = cachedJointAverage(a, b, c, d); const double pefgh = cachedJointAverage(e, f, g, h); const double pabef = cachedJointAverage(a, b, e, f); const double pcdgh = cachedJointAverage(c, d, g, h); const double pabgh = cachedJointAverage(a, b, g, h); const double pcdef = cachedJointAverage(c, d, e, f); const double pabcdef = cachedJointAverage(a, b, c, d, e, f); const double pabcdgh = cachedJointAverage(a, b, c, d, g, h); const double pabefgh = cachedJointAverage(a, b, e, f, g, h); const double pcdefgh = cachedJointAverage(c, d, e, f, g, h); const double pabcdefgh = cachedJointAverage(a, b, c, d, e, f, g, h); const double deltaprod = kdelta(a, b)*kdelta(c, d)*kdelta(e, f)*kdelta(g, h); const double tmp = kdelta(e,f)*kdelta(g,h)*pabcd + kdelta(c,d)*kdelta(g,h)*pabef + kdelta(c,d)*kdelta(e,f)*pabgh + kdelta(a,b)*kdelta(c,d)*pefgh + kdelta(a,b)*kdelta(e,f)*pcdgh + kdelta(a,b)*kdelta(g,h)*pcdef; const double term1 = pabcd*pefgh + pabef*pcdgh + pabgh*pcdef + 3.0*deltaprod - tmp; const double term2 = pabcdefgh - 6.0*deltaprod - kdelta(a,b)*pcdefgh - kdelta(c,d)*pabefgh - kdelta(e,f)*pabcdgh - kdelta(g,h)*pabcdef - pabcd*pefgh - pabef*pcdgh - pabgh*pcdef + 2.0*tmp; const double nPoints = effectiveSampleSize(); const double prod8 = (term1 + term2/nPoints)/nPoints/nPoints; cov = prod8 - cov4(a, b, c, d)*cov4(e, f, g, h); } return cov; } double ContOrthoPoly1D::cov8(const unsigned a, const unsigned b, const unsigned c, const unsigned d, const unsigned e, const unsigned f, const unsigned g, const unsigned h) const { const bool has0 = (a == 0 && b == 0) || (c == 0 && d == 0) || (e == 0 && f == 0) || (g == 0 && h == 0); double cov = 0.0; if (!has0) { const bool includeCubicPart = true; // First, calculate the O(N^-2) part const double pabef = cachedJointAverage(a, b, e, f); const double pcdgh = cachedJointAverage(c, d, g, h); const double pabgh = cachedJointAverage(a, b, g, h); const double pcdef = cachedJointAverage(c, d, e, f); const double deltaprod = kdelta(a, b)*kdelta(c, d)*kdelta(e, f)*kdelta(g, h); const double tmp2 = kdelta(c, d)*kdelta(g, h)*pabef + kdelta(c, d)*kdelta(e, f)*pabgh + kdelta(a, b)*kdelta(e, f)*pcdgh + kdelta(a, b)*kdelta(g, h)*pcdef; const double term2 = pabef*pcdgh + pabgh*pcdef + 2.0*deltaprod - tmp2; double term3 = 0.0; if (includeCubicPart) { // Additional terms needed to calculate the O(N^-3) part const double pabcdefgh = cachedJointAverage(a, b, c, d, e, f, g, h); double sixsum = 0.0; if (kdelta(a, b)) sixsum += cachedJointAverage(c, d, e, f, g, h); if (kdelta(c, d)) sixsum += cachedJointAverage(a, b, e, f, g, h); if (kdelta(e, f)) sixsum += cachedJointAverage(a, b, c, d, g, h); if (kdelta(g, h)) sixsum += cachedJointAverage(a, b, c, d, e, f); const double pabcd = cachedJointAverage(a, b, c, d); const double pefgh = cachedJointAverage(e, f, g, h); const double tmp3 = tmp2 + kdelta(e, f)*kdelta(g, h)*pabcd + kdelta(a, b)*kdelta(c, d)*pefgh; term3 = pabcdefgh - 6.0*deltaprod - sixsum - (pabcd*pefgh + pabef*pcdgh + pabgh*pcdef) + 2.0*tmp3; } const double N = effectiveSampleSize(); cov = (term2 + term3/N)/N/N; // const bool isVariance = ?; // if (isVariance) // if (cov < 0.0) // cov = 0.0; } return cov; } double ContOrthoPoly1D::epsExpectation(const unsigned m_in, const unsigned n_in, const bool highOrder) const { if (highOrder) { long double sum = 0.0L; if (m_in || n_in) { const unsigned m = std::min(m_in, n_in); const unsigned n = std::max(m_in, n_in); for (unsigned k=0; k<=n; ++k) { const double f = k == n ? 1.0 : (k > m ? 1.0 : 2.0); sum += f*cachedJointAverage(k, k, m, n); } if (m == n) sum -= 1.0; else sum -= (cachedJointAverage(m, m, m, n) + cachedJointAverage(m, n, n, n))/2.0; } return sum/effectiveSampleSize(); } else return 0.0; } double ContOrthoPoly1D::epsCovariance(const unsigned m1_in, const unsigned n1_in, const unsigned m2_in, const unsigned n2_in, const bool highOrder) const { const bool has0 = (m1_in == 0 && n1_in == 0) || (m2_in == 0 && n2_in == 0); if (has0) return 0.0; if (highOrder) { const unsigned m1 = std::min(m1_in, n1_in); const unsigned n1 = std::max(m1_in, n1_in); const unsigned m2 = std::min(m2_in, n2_in); const unsigned n2 = std::max(m2_in, n2_in); long double sum = 0.0L; // Process the -v_{m1,n1} term (i.e., the linear one) of eps_{m1,n1} sum += cov4(m2, n2, m1, n1); sum += cov6(m2, n2, n2, n2, m1, n1)/2.0; sum += cov6(m2, n2, m2, m2, m1, n1)/2.0; for (unsigned k=0; k<=n2; ++k) { const double factor = k > m2 ? 1.0 : 2.0; sum -= factor*cov6(k, m2, k, n2, m1, n1); } // Process the term -v_{m1,n1}/2 (v_{n1,n1} + v_{m1,m1}) of eps_{m1,n1} unsigned idx[2]; idx[0] = n1; idx[1] = m1; for (unsigned ii=0; ii<2; ++ii) { const unsigned mOrN = idx[ii]; sum += cov6(m1, n1, mOrN, mOrN, m2, n2)/2.0; sum += cov8(m1, n1, mOrN, mOrN, m2, n2, n2, n2)/4.0; sum += cov8(m1, n1, mOrN, mOrN, m2, n2, m2, m2)/4.0; for (unsigned k=0; k<=n2; ++k) { const double factor = k > m2 ? 1.0 : 2.0; sum -= factor/2.0*cov8(m1, n1, mOrN, mOrN, k, m2, k, n2); } } // Process the sum in eps_{m1,n1} for (unsigned k1=0; k1<=n1; ++k1) { const double f1 = k1 > m1 ? 1.0 : 2.0; sum -= f1*cov6(k1, m1, k1, n1, m2, n2); sum -= f1*cov8(k1, m1, k1, n1, m2, n2, n2, n2)/2.0; sum -= f1*cov8(k1, m1, k1, n1, m2, n2, m2, m2)/2.0; for (unsigned k=0; k<=n2; ++k) { const double factor = k > m2 ? 1.0 : 2.0; sum += f1*factor*cov8(k1, m1, k1, n1, k, m2, k, n2); } } return sum; } else return cov4(m2_in, n2_in, m1_in, n1_in); } }