Index: trunk/npstat/stat/orthoPoly1DVProducts.icc =================================================================== --- trunk/npstat/stat/orthoPoly1DVProducts.icc (revision 567) +++ trunk/npstat/stat/orthoPoly1DVProducts.icc (revision 568) @@ -1,467 +1,522 @@ #include #include #include #include #include namespace npstat { namespace Private { inline void validateNPoints(const char* where, const unsigned long nPoints) { if (!nPoints) { std::string message("In npstat::"); message += where; message += ": sample size must be positive"; throw std::invalid_argument(message); } } inline int kdelta(const UUPair ab) { return ab.first == ab.second ? 1 : 0; } inline bool isZeroPair(const UUPair pair) { return !(pair.first || pair.second); } } template double expectedVProduct(const AbsClassicalOrthoPoly1D& poly, const Quadrature& quad, const unsigned long nPoints, const UUPair ab, const UUPair cd) { Private::validateNPoints("expectedVProduct", nPoints); if (Private::isZeroPair(ab) || Private::isZeroPair(cd)) return 0.0; else return (poly.jointAverage(quad, ab.first, ab.second, cd.first, cd.second) - Private::kdelta(ab)*Private::kdelta(cd))/nPoints; } template inline double expectedVCovariance(const AbsClassicalOrthoPoly1D& poly, const Quadrature& quad, const unsigned long nPoints, const UUPair ab, const UUPair cd) { Private::validateNPoints("expectedVCovariance", nPoints); return expectedVProduct(poly, quad, nPoints, ab, cd); } template double expectedVProduct(const AbsClassicalOrthoPoly1D& poly, const Quadrature& quad, const unsigned long nPoints, const UUPair ab, const UUPair cd, const UUPair ef) { Private::validateNPoints("expectedVProduct", nPoints); if (Private::isZeroPair(ab) || Private::isZeroPair(cd) || Private::isZeroPair(ef)) return 0.0; const double pcdef = poly.jointAverage(quad, cd.first, cd.second, ef.first, ef.second); const double pabef = poly.jointAverage(quad, ab.first, ab.second, ef.first, ef.second); const double pabcd = poly.jointAverage(quad, ab.first, ab.second, cd.first, cd.second); unsigned degs[6]; degs[0] = ab.first; degs[1] = ab.second; degs[2] = cd.first; degs[3] = cd.second; degs[4] = ef.first; degs[5] = ef.second; const double pabcdef = poly.jointAverage(quad, degs, 6); const int deltaprod = Private::kdelta(ab)*Private::kdelta(cd)* Private::kdelta(ef); return (pabcdef - Private::kdelta(ab)*pcdef - Private::kdelta(cd)*pabef - Private::kdelta(ef)*pabcd + 2.0*deltaprod)/nPoints/nPoints; } template inline double expectedVCovariance(const AbsClassicalOrthoPoly1D& poly, const Quadrature& quad, const unsigned long nPoints, const UUPair ab, const UUPair cd, const UUPair ef) { Private::validateNPoints("expectedVCovariance", nPoints); return expectedVProduct(poly, quad, nPoints, ab, cd, ef); } template double expectedVProduct(const AbsClassicalOrthoPoly1D& poly, const Quadrature& quad, const unsigned long nPoints, const UUPair ab, const UUPair cd, const UUPair ef, const UUPair gh) { Private::validateNPoints("expectedVProduct", nPoints); if (Private::isZeroPair(ab) || Private::isZeroPair(cd) || Private::isZeroPair(ef) || Private::isZeroPair(gh)) return 0.0; const double pabcd = poly.jointAverage(quad, ab.first, ab.second, cd.first, cd.second); const double pefgh = poly.jointAverage(quad, ef.first, ef.second, gh.first, gh.second); const double pabef = poly.jointAverage(quad, ab.first, ab.second, ef.first, ef.second); const double pcdgh = poly.jointAverage(quad, cd.first, cd.second, gh.first, gh.second); const double pabgh = poly.jointAverage(quad, ab.first, ab.second, gh.first, gh.second); const double pcdef = poly.jointAverage(quad, cd.first, cd.second, ef.first, ef.second); unsigned degs[8]; degs[0] = ab.first; degs[1] = ab.second; degs[2] = cd.first; degs[3] = cd.second; degs[4] = ef.first; degs[5] = ef.second; degs[6] = gh.first; degs[7] = gh.second; const double pabcdefgh = poly.jointAverage(quad, degs, 8); const double pabcdef = poly.jointAverage(quad, degs, 6); degs[4] = gh.first; degs[5] = gh.second; const double pabcdgh = poly.jointAverage(quad, degs, 6); degs[2] = ef.first; degs[3] = ef.second; const double pabefgh = poly.jointAverage(quad, degs, 6); degs[0] = cd.first; degs[1] = cd.second; const double pcdefgh = poly.jointAverage(quad, degs, 6); const int deltaprod = Private::kdelta(ab)*Private::kdelta(cd)* Private::kdelta(ef)*Private::kdelta(gh); const double tmp = Private::kdelta(ef)*Private::kdelta(gh)*pabcd + Private::kdelta(cd)*Private::kdelta(gh)*pabef + Private::kdelta(cd)*Private::kdelta(ef)*pabgh + Private::kdelta(ab)*Private::kdelta(cd)*pefgh + Private::kdelta(ab)*Private::kdelta(ef)*pcdgh + Private::kdelta(ab)*Private::kdelta(gh)*pcdef; const double term1 = pabcd*pefgh + pabef*pcdgh + pabgh*pcdef + 3.0*deltaprod - tmp; const double term2 = pabcdefgh - 6.0*deltaprod - Private::kdelta(ab)*pcdefgh - Private::kdelta(cd)*pabefgh - Private::kdelta(ef)*pabcdgh - Private::kdelta(gh)*pabcdef - pabcd*pefgh - pabef*pcdgh - pabgh*pcdef + 2.0*tmp; return (term1 + term2/nPoints)/nPoints/nPoints; } template inline double expectedVCovariance(const AbsClassicalOrthoPoly1D& poly, const Quadrature& quad, const unsigned long nPoints, const UUPair ab, const UUPair cd, const UUPair ef, const UUPair gh) { Private::validateNPoints("expectedVCovariance", nPoints); return expectedVProduct(poly, quad, nPoints, ab, cd, ef, gh) - expectedVProduct(poly, quad, nPoints, ab, cd)* expectedVProduct(poly, quad, nPoints, ef, gh); } - template + template inline double sampleVProduct(const AbsClassicalOrthoPoly1D& poly, const Numeric* coords, const unsigned long nPoints, const UUPair ab) { Private::validateNPoints("sampleVProduct", nPoints); assert(coords); if (Private::isZeroPair(ab)) return 0.0; else { unsigned degs[2]; degs[0] = ab.first; degs[1] = ab.second; return poly.jointSampleAverage(coords, nPoints, degs, 2) - Private::kdelta(ab); } } - template + template inline double sampleVProduct(const AbsClassicalOrthoPoly1D& poly, const Numeric* coords, const unsigned long nPoints, const UUPair ab, const UUPair cd) { Private::validateNPoints("sampleVProduct", nPoints); assert(coords); if (Private::isZeroPair(ab) || Private::isZeroPair(cd)) return 0.0; else { unsigned degs1[2]; degs1[0] = ab.first; degs1[1] = ab.second; unsigned degs2[2]; degs2[0] = cd.first; degs2[1] = cd.second; const std::pair& dp = poly.twoJointSampleAverages( coords, nPoints, degs1, 2, degs2, 2); return (dp.first - Private::kdelta(ab))* (dp.second - Private::kdelta(cd)); } } - template + template + inline double sampleVCovariance(const AbsClassicalOrthoPoly1D& poly, + const Numeric* points, + const unsigned long nPoints, + const UUPair ab, const UUPair cd) + { + Private::validateNPoints("sampleVCovariance", nPoints); + return sampleVProduct(poly, points, nPoints, ab, cd); + } + + template inline double sampleVProduct(const AbsClassicalOrthoPoly1D& poly, const Numeric* coords, const unsigned long nPoints, const UUPair ab, const UUPair cd, const UUPair ef) { Private::validateNPoints("sampleVProduct", nPoints); assert(coords); if (Private::isZeroPair(ab) || Private::isZeroPair(cd) || Private::isZeroPair(ef)) return 0.0; else return sampleVProduct(poly, coords, nPoints, ab, cd)* sampleVProduct(poly, coords, nPoints, ef); } - template + template + inline double sampleVCovariance(const AbsClassicalOrthoPoly1D& poly, + const Numeric* points, + const unsigned long nPoints, + const UUPair ab, const UUPair cd, + const UUPair ef) + { + Private::validateNPoints("sampleVCovariance", nPoints); + return sampleVProduct(poly, points, nPoints, ab, cd, ef); + } + + template inline double sampleVProduct(const AbsClassicalOrthoPoly1D& poly, const Numeric* coords, const unsigned long nPoints, const UUPair ab, const UUPair cd, const UUPair ef, const UUPair gh) { Private::validateNPoints("sampleVProduct", nPoints); assert(coords); if (Private::isZeroPair(ab) || Private::isZeroPair(cd) || Private::isZeroPair(ef) || Private::isZeroPair(gh)) return 0.0; else return sampleVProduct(poly, coords, nPoints, ab, cd)* sampleVProduct(poly, coords, nPoints, ef, gh); } + template + inline double sampleVCovariance(const AbsClassicalOrthoPoly1D& poly, + const Numeric* points, + const unsigned long nPoints, + const UUPair ab, const UUPair cd, + const UUPair ef, const UUPair gh) + { + Private::validateNPoints("sampleVCovariance", nPoints); + return sampleVProduct(poly, points, nPoints, ab, cd, ef, gh) - + sampleVProduct(poly, points, nPoints, ab, cd)* + sampleVProduct(poly, points, nPoints, ef, gh); + } + template double oracleEpsExpectation(const AbsClassicalOrthoPoly1D& poly, const Quadrature& quad, const unsigned long nPoints, const unsigned m_in, const unsigned n_in, const bool highOrder) { if (highOrder) { Private::validateNPoints("oracleEpsExpectation", nPoints); 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*poly.jointAverage(quad, k, k, m, n); } if (m == n) sum -= 1.0; else sum -= (poly.jointAverage(quad, m, m, m, n) + poly.jointAverage(quad, m, n, n, n))/2.0; } return sum/nPoints; } else return 0.0; } - template + template double sampleEpsValue(const AbsClassicalOrthoPoly1D& poly, const Numeric* coords, unsigned long nPoints, const unsigned m_in, const unsigned n_in, const bool highOrder) { Private::validateNPoints("sampleEpsValue", nPoints); assert(coords); 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); const double vmn = sampleVProduct(poly, coords, nPoints, UUPair(m, n)); if (m == n) sum -= (vmn + vmn*vmn); else { const double vnn = sampleVProduct(poly, coords, nPoints, UUPair(n, n)); const double vmm = sampleVProduct(poly, coords, nPoints, UUPair(m, m)); sum -= (vmn + 0.5*vmn*(vnn + vmm)); } for (unsigned k=0; k<=n; ++k) { const double f = k > m ? 1.0 : 2.0; const double vkmkn = sampleVProduct(poly, coords, nPoints, UUPair(k, m), UUPair(k, n)); sum += f*vkmkn; } } return sum; } else return -sampleVProduct(poly, coords, nPoints, UUPair(m_in, n_in)); } +#define npstat_EpsCovarianceCalculation_mainbody(calculator, argument) do { \ + 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 */ \ + sum += calculator (poly, argument, nPoints, \ + UUPair(m2, n2), UUPair(m1, n1)); \ + sum += calculator (poly, argument, nPoints, \ + UUPair(m2, n2), UUPair(n2, n2), UUPair(m1, n1))/2.0; \ + sum += calculator (poly, argument, nPoints, \ + UUPair(m2, n2), UUPair(m2, m2), UUPair(m1, n1))/2.0; \ + for (unsigned k=0; k<=n2; ++k) \ + { \ + const double factor = k > m2 ? 1.0 : 2.0; \ + sum -= factor*calculator (poly, argument, nPoints, \ + UUPair(k, m2), UUPair(k, n2), UUPair(m1, n1));\ + } \ + /* Process the term -v_{m1,n1}/2 (v_{n1,n1} + v_{m1,m1}) */ \ + unsigned idx[2]; \ + idx[0] = n1; \ + idx[1] = m1; \ + for (unsigned ii=0; ii<2; ++ii) \ + { \ + const unsigned mOrN = idx[ii]; \ + sum += calculator (poly, argument, nPoints, \ + UUPair(m1, n1), UUPair(mOrN, mOrN), \ + UUPair(m2, n2))/2.0; \ + sum += calculator (poly, argument, nPoints, \ + UUPair(m1, n1), UUPair(mOrN, mOrN), \ + UUPair(m2, n2), UUPair(n2, n2))/4.0; \ + sum += calculator (poly, argument, nPoints, \ + UUPair(m1, n1), UUPair(mOrN, mOrN), \ + UUPair(m2, n2), UUPair(m2, m2))/4.0; \ + for (unsigned k=0; k<=n2; ++k) \ + { \ + const double factor = k > m2 ? 1.0 : 2.0; \ + sum -= factor/2.0*calculator (poly, argument, nPoints, \ + UUPair(m1, n1), UUPair(mOrN, mOrN), \ + UUPair(k, m2), UUPair(k, n2)); \ + } \ + } \ + /* Process the sum */ \ + for (unsigned k1=0; k1<=n1; ++k1) \ + { \ + const double f1 = k1 > m1 ? 1.0 : 2.0; \ + sum -= f1*calculator (poly, argument, nPoints, \ + UUPair(k1, m1), UUPair(k1, n1), UUPair(m2, n2)); \ + sum -= f1*calculator (poly, argument, nPoints, \ + UUPair(k1, m1), UUPair(k1, n1), \ + UUPair(m2, n2), UUPair(n2, n2))/2.0; \ + sum -= f1*calculator (poly, argument, nPoints, \ + UUPair(k1, m1), UUPair(k1, n1), \ + UUPair(m2, n2), UUPair(m2, m2))/2.0; \ + for (unsigned k=0; k<=n2; ++k) \ + { \ + const double factor = k > m2 ? 1.0 : 2.0; \ + sum += f1*factor*calculator (poly, argument, nPoints, \ + UUPair(k1, m1), UUPair(k1, n1), \ + UUPair(k, m2), UUPair(k, n2)); \ + } \ + } \ + return sum; \ + } \ + else \ + return calculator (poly, argument, nPoints, \ + UUPair(m2_in, n2_in), UUPair(m1_in, n1_in)); \ +} while(0); + template double oracleEpsCovariance(const AbsClassicalOrthoPoly1D& poly, const Quadrature& quad, const unsigned long nPoints, const unsigned m1_in, const unsigned n1_in, const unsigned m2_in, const unsigned n2_in, const bool highOrder) { Private::validateNPoints("oracleEpsCovariance", nPoints); - 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 - sum += expectedVCovariance(poly, quad, nPoints, - UUPair(m2, n2), UUPair(m1, n1)); - sum += expectedVCovariance(poly, quad, nPoints, - UUPair(m2, n2), UUPair(n2, n2), UUPair(m1, n1))/2.0; - sum += expectedVCovariance(poly, quad, nPoints, - UUPair(m2, n2), UUPair(m2, m2), UUPair(m1, n1))/2.0; - for (unsigned k=0; k<=n2; ++k) - { - const double factor = k > m2 ? 1.0 : 2.0; - sum -= factor*expectedVCovariance(poly, quad, nPoints, - UUPair(k, m2), UUPair(k, n2), UUPair(m1, n1)); - } - - // Process the term -v_{m1,n1}/2 (v_{n1,n1} + v_{m1,m1}) - unsigned idx[2]; - idx[0] = n1; - idx[1] = m1; - for (unsigned ii=0; ii<2; ++ii) - { - const unsigned mOrN = idx[ii]; - sum += expectedVCovariance(poly, quad, nPoints, - UUPair(m1, n1), UUPair(mOrN, mOrN), - UUPair(m2, n2))/2.0; - sum += expectedVCovariance(poly, quad, nPoints, - UUPair(m1, n1), UUPair(mOrN, mOrN), - UUPair(m2, n2), UUPair(n2, n2))/4.0; - sum += expectedVCovariance(poly, quad, nPoints, - UUPair(m1, n1), UUPair(mOrN, mOrN), - UUPair(m2, n2), UUPair(m2, m2))/4.0; - for (unsigned k=0; k<=n2; ++k) - { - const double factor = k > m2 ? 1.0 : 2.0; - sum -= factor/2.0*expectedVCovariance(poly, quad, nPoints, - UUPair(m1, n1), UUPair(mOrN, mOrN), - UUPair(k, m2), UUPair(k, n2)); - } - } - - // Process the sum - for (unsigned k1=0; k1<=n1; ++k1) - { - const double f1 = k1 > m1 ? 1.0 : 2.0; - sum -= f1*expectedVCovariance(poly, quad, nPoints, - UUPair(k1, m1), UUPair(k1, n1), UUPair(m2, n2)); - sum -= f1*expectedVCovariance(poly, quad, nPoints, - UUPair(k1, m1), UUPair(k1, n1), - UUPair(m2, n2), UUPair(n2, n2))/2.0; - sum -= f1*expectedVCovariance(poly, quad, nPoints, - UUPair(k1, m1), UUPair(k1, n1), - UUPair(m2, n2), UUPair(m2, m2))/2.0; - for (unsigned k=0; k<=n2; ++k) - { - const double factor = k > m2 ? 1.0 : 2.0; - sum += f1*factor*expectedVCovariance(poly, quad, nPoints, - UUPair(k1, m1), UUPair(k1, n1), - UUPair(k, m2), UUPair(k, n2)); - } - } + npstat_EpsCovarianceCalculation_mainbody(expectedVCovariance, quad); + } - return sum; - } - else - return expectedVCovariance(poly, quad, nPoints, - UUPair(m2_in, n2_in), UUPair(m1_in, n1_in)); + template + double sampleEpsCovariance(const AbsClassicalOrthoPoly1D& poly, + const Numeric* coords, + const unsigned long nPoints, + const unsigned m1_in, const unsigned n1_in, + const unsigned m2_in, const unsigned n2_in, + const bool highOrder) + { + Private::validateNPoints("sampleEpsCovariance", nPoints); + npstat_EpsCovarianceCalculation_mainbody(sampleVCovariance, coords); } +#define npstat_EpsCovarianceArray_mainbody(calculator, argument) do { \ + const unsigned mdp1 = maxdeg + 1; \ + ArrayND arr(mdp1, mdp1, mdp1, mdp1); \ + unsigned long nDefined = 0; \ + for (unsigned i=0; i<=maxdeg; ++i) \ + for (unsigned j=0; j<=i; ++j) \ + for (unsigned k=0; k<=maxdeg; ++k) \ + for (unsigned m=0; m<=k; ++m) \ + if (std::make_pair(j,i) <= std::make_pair(m,k)) \ + { \ + arr(j, i, m, k) = calculator ( \ + poly, argument, nPoints, j, i, m, k, highOrder); \ + ++nDefined; \ + } \ + for (unsigned i=0; i<=maxdeg; ++i) \ + for (unsigned j=0; j<=maxdeg; ++j) \ + for (unsigned k=0; k<=maxdeg; ++k) \ + for (unsigned m=0; m<=maxdeg; ++m) \ + if (!(j<=i && m<=k && std::make_pair(j,i) <= std::make_pair(m,k)))\ + { \ + unsigned i1 = i; \ + unsigned j1 = j; \ + unsigned k1 = k; \ + unsigned m1 = m; \ + if (!(j1 <= i1)) \ + std::swap(i1, j1); \ + if (!(m1 <= k1)) \ + std::swap(m1, k1); \ + if (!(std::make_pair(j1,i1) <= std::make_pair(m1,k1))) \ + { \ + std::swap(m1, j1); \ + std::swap(k1, i1); \ + } \ + arr(j, i, m, k) = arr(j1, i1, m1, k1); \ + ++nDefined; \ + } \ + assert(nDefined == arr.length()); \ + return arr; \ +} while (0); + template ArrayND oracleEpsCovarianceArray( const AbsClassicalOrthoPoly1D& poly, const Quadrature& quad, const unsigned long nPoints, const unsigned maxdeg, const bool highOrder) { - const unsigned mdp1 = maxdeg + 1; - ArrayND arr(mdp1, mdp1, mdp1, mdp1); - unsigned long nDefined = 0; - - for (unsigned i=0; i<=maxdeg; ++i) - for (unsigned j=0; j<=i; ++j) - for (unsigned k=0; k<=maxdeg; ++k) - for (unsigned m=0; m<=k; ++m) - if (std::make_pair(j,i) <= std::make_pair(m,k)) - { - arr(j, i, m, k) = oracleEpsCovariance( - poly, quad, nPoints, j, i, m, k, highOrder); - ++nDefined; - } - - for (unsigned i=0; i<=maxdeg; ++i) - for (unsigned j=0; j<=maxdeg; ++j) - for (unsigned k=0; k<=maxdeg; ++k) - for (unsigned m=0; m<=maxdeg; ++m) - if (!(j<=i && m<=k && std::make_pair(j,i) <= std::make_pair(m,k))) - { - unsigned i1 = i; - unsigned j1 = j; - unsigned k1 = k; - unsigned m1 = m; - if (!(j1 <= i1)) - std::swap(i1, j1); - if (!(m1 <= k1)) - std::swap(m1, k1); - if (!(std::make_pair(j1,i1) <= std::make_pair(m1,k1))) - { - std::swap(m1, j1); - std::swap(k1, i1); - } - arr(j, i, m, k) = arr(j1, i1, m1, k1); - ++nDefined; - } + Private::validateNPoints("oracleEpsCovarianceArray", nPoints); + npstat_EpsCovarianceArray_mainbody(oracleEpsCovariance, quad); + } - assert(nDefined == arr.length()); - return arr; + template + ArrayND sampleEpsCovarianceArray( + const AbsClassicalOrthoPoly1D& poly, const Numeric* coords, + const unsigned long nPoints, const unsigned maxdeg, const bool highOrder) + { + Private::validateNPoints("sampleEpsCovarianceArray", nPoints); + npstat_EpsCovarianceArray_mainbody(sampleEpsCovariance, coords); } } Index: trunk/npstat/stat/orthoPoly1DVProducts.hh =================================================================== --- trunk/npstat/stat/orthoPoly1DVProducts.hh (revision 567) +++ trunk/npstat/stat/orthoPoly1DVProducts.hh (revision 568) @@ -1,125 +1,163 @@ #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 + template double sampleVProduct(const AbsClassicalOrthoPoly1D& poly, - const Numeric* coords, unsigned long nPoints, + const Numeric* coords, unsigned long nCoords, UUPair ab); // v_{ab}*v_{cd} calculated for a given sample of unweighted points - template + template double sampleVProduct(const AbsClassicalOrthoPoly1D& poly, - const Numeric* coords, unsigned long nPoints, + const Numeric* coords, unsigned long nCoords, UUPair ab, UUPair cd); // v_{ab}*v_{cd}*v_{ef} calculated for a given sample - template + template double sampleVProduct(const AbsClassicalOrthoPoly1D& poly, - const Numeric* coords, unsigned long nPoints, + const Numeric* coords, unsigned long nCoords, UUPair ab, UUPair cd, UUPair ef); // v_{ab}*v_{cd}*v_{ef}*v_{gh} calculated for a given sample - template + template double sampleVProduct(const AbsClassicalOrthoPoly1D& poly, - const Numeric* coords, unsigned long nPoints, + const Numeric* coords, unsigned long nCoords, UUPair ab, UUPair cd, UUPair ef, UUPair gh); + // Covariance between v_{ab} and v_{cd} estimated using + // polynomial averages taken from the given sample + template + double sampleVCovariance(const AbsClassicalOrthoPoly1D& poly, + const Numeric* coords, unsigned long nCoords, + UUPair ab, UUPair cd); + + // Covariance between v_{ab}*v_{cd} and v_{ef} estimated using + // polynomial averages taken from the given sample + template + double sampleVCovariance(const AbsClassicalOrthoPoly1D& poly, + const Numeric* coords, unsigned long nCoords, + UUPair ab, UUPair cd, UUPair ef); + + // Covariance between v_{ab}*v_{cd} and v_{ef}*v_{gh} estimated + // using polynomial averages taken from the given sample + template + double sampleVCovariance(const AbsClassicalOrthoPoly1D& poly, + const Numeric* coords, unsigned long nCoords, + 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 + // covariances are symmetric with respect to permuting some 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_{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 + template double sampleEpsValue(const AbsClassicalOrthoPoly1D& poly, - const Numeric* coords, unsigned long nPoints, + const Numeric* coords, unsigned long nCoords, unsigned m, unsigned n, bool highOrder); + + // eps_{mn} covariance approximation using polynomial product + // averages from the given sample of points + template + double sampleEpsCovariance(const AbsClassicalOrthoPoly1D& poly, + const Numeric* coords, unsigned long nCoords, + unsigned m1, unsigned n1, + unsigned m2, unsigned n2, bool highOrder); + + // Fill a 4-d array with sample covariances. The code knows that the + // covariances are symmetric with respect to permuting some indices, and + // does not recalculate them when permuting indices is sufficient. The + // array span will be maxdeg+1 in each dimension. + template + ArrayND sampleEpsCovarianceArray( + const AbsClassicalOrthoPoly1D& poly, const Numeric* coords, + unsigned long nCoords, unsigned maxdeg, bool highOrder); } #include "npstat/stat/orthoPoly1DVProducts.icc" #endif // NPSTAT_ORTHOPOLY1DVPRODUCTS_HH_