Index: trunk/npstat/nm/ContOrthoPoly1D.cc =================================================================== --- trunk/npstat/nm/ContOrthoPoly1D.cc (revision 709) +++ trunk/npstat/nm/ContOrthoPoly1D.cc (revision 710) @@ -1,960 +1,1033 @@ +#include #include #include "npstat/nm/ContOrthoPoly1D.hh" #include "npstat/nm/PairCompare.hh" #include "npstat/nm/allocators.hh" #include "npstat/nm/StorablePolySeries1D.hh" +#include "npstat/nm/Triple.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& 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; 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::weightCoeffs(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 { // Process various simple special cases first if (!power1) return powerAverage(ideg2, power2); if (!power2) return powerAverage(ideg1, power1); if (power1 == 1U && power2 == 1U) return kdelta(ideg1, ideg2); // General calculation 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); // See if we can avoid a direct calculation if (nDegrees == 1U) return kdelta(degrees[0], 0U); if (nDegrees == 2U) return kdelta(degrees[0], degrees[1]); // Check if we can remove leading zeros { unsigned nonZeroDegs = nDegrees; while (nonZeroDegs && !degrees[0]) { ++degrees; --nonZeroDegs; } if (nonZeroDegs < nDegrees) return jointAverage(degrees, nonZeroDegs, sorted); } // Check if we can remove any other zeros if (!sorted) { bool allNonZero = true; for (unsigned ideg=0; ideg 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_); } + Matrix ContOrthoPoly1D::tripleProductMatrix( + const std::vector >& pairs, + const unsigned maxdeg, const double xmin, const double xmax) const + { + typedef npstat::Triple Key3; + typedef std::map TripleProductMap; + + const unsigned nPairs = pairs.size(); + if (!nPairs) throw std::invalid_argument( + "In npstat::ContOrthoPoly1D::tripleProductMatrix: " + "empty vector of degree specifications"); + const unsigned nPolys = maxdeg + 1U; + + npstat::Matrix design0(nPairs, nPolys); + TripleProductMap tmap; + + // Avoid calculating the same integral more than once + for (unsigned ipair=0; ipair where = + tmap.insert(std::make_pair(key, value)); + assert(where.second); + it = where.first; + } + design0[ipair][ipoly] = it->second; + } + } + + return design0; + } + 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::covCov4(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 pabcdefgh = cachedJointAverage(a, b, c, d, e, f, g, h); const double pabcd = cachedJointAverage(a, b, c, d); const double pefgh = cachedJointAverage(e, f, g, h); const double N = effectiveSampleSize(); cov = (pabcdefgh - pabcd*pefgh)/N/N/N; } 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); } + Matrix ContOrthoPoly1D::epsCovarianceMatrix( + const std::vector >& pairs, + const bool highOrder) const + { + if (pairs.empty()) throw std::invalid_argument( + "In npstat::ContOrthoPoly1D::epsCovarianceMatrix: empty list of pairs"); + + const unsigned nPairs = pairs.size(); + Matrix cov(nPairs, nPairs); + for (unsigned ipair1=0; ipair1 ContOrthoPoly1D::makeStorablePolySeries( const double i_xmin, const double i_xmax, const double *coeffs, const unsigned maxdeg) const { const unsigned sz = rCoeffs_.size(); std::vector > rc(sz); for (unsigned i=0; i( new StorablePolySeries1D(rc, i_xmin, i_xmax, meanX_, coeffs, maxdeg)); } } Index: trunk/npstat/nm/Matrix.hh =================================================================== --- trunk/npstat/nm/Matrix.hh (revision 709) +++ trunk/npstat/nm/Matrix.hh (revision 710) @@ -1,826 +1,827 @@ #ifndef NPSTAT_MATRIX_HH_ #define NPSTAT_MATRIX_HH_ /*! // \file Matrix.hh // // \brief Template matrix class // // Author: I. Volobouev // // November 2008 */ #include #include #include #include #include "geners/ClassId.hh" #include "npstat/nm/EigenMethod.hh" #include "npstat/nm/SvdMethod.hh" #include "npstat/nm/GeneralizedComplex.hh" namespace npstat { /** // A simple helper class for matrix manipulations. Depending on how much // space is provided with the "Len" parameter, the data will be placed // either on the stack or on the heap. Do not set "Len" to 0, the code // assumes that the stack space is available for at least one element. // // C storage convention is used for internal data. In the element access // operations, array bounds are not checked. // // A number of operations (solving linear systems, calculating eigenvalues // and eigenvectors, singular value decomposition, etc) are performed // by calling appropriate routines from LAPACK. This usually limits // the Numeric template parameter types for which these operations are // available to float and double. If the operation is unavailable for // the template parameter used, std::invalid_argument exception is raised. // // Note that for simple matrix operations this class is likely to be // slower than matrix classes based on expression templates (such as // those in boost uBLAS or in Blitz++). If speed is really important // for your calculations, consider using a dedicated matrix library. */ template class Matrix { template friend class Matrix; public: typedef Numeric value_type; /** // Default constructor creates an unitialized matrix // which can be assigned from other matrices */ Matrix(); /** // This constructor creates an unitialized matrix // which can be assigned element-by-element or from // another matrix with the same dimensions */ Matrix(unsigned nrows, unsigned ncols); /** // This constructor initializes the matrix as follows: // // initCode = 0: all elements are initialized to 0. // // initCode = 1: matrix must be square; diagonal elements are // initialized to 1, all other elements to 0. */ Matrix(unsigned nrows, unsigned ncols, int initCode); /** // This constructor initializes the matrix from the // given 1-d array using C storage conventions */ Matrix(unsigned nrows, unsigned ncols, const Numeric* data); /** Copy constructor */ Matrix(const Matrix&); /** Move constructor */ Matrix(Matrix&&); /** Converting copy constructor */ template Matrix(const Matrix&); /** // Constructor from a subrange of another matrix. The minimum row // and column numbers are included, maximum numbers are excluded. */ template Matrix(const Matrix&, unsigned rowMin, unsigned rowMax, unsigned columnMin, unsigned columnMax); ~Matrix(); /** // Assignment operator. The matrix on the left must either be // in an uninitialized state or have compatible dimensions with // the matrix on the right. */ Matrix& operator=(const Matrix&); /** Move assignment operator */ Matrix& operator=(Matrix&&); /** Converting assignment operator */ template Matrix& operator=(const Matrix&); /** Set from triplets (for compatibility with Eigen) */ template Matrix& setFromTriplets(Iterator first, Iterator last); /** Set all data from an external array */ template Matrix& setData(const Num2* data, unsigned dataLength); /** Set a row from an external array */ template Matrix& setRow(unsigned row, const Num2* data, unsigned dataLength); /** Set a column from an external array */ template Matrix& setColumn(unsigned col, const Num2* data, unsigned dataLength); /** // Tag this matrix as diagonal. This may improve the speed // of certain operations. Of course, this should be done only // if you know for sure that this matrix is, indeed, diagonal. // Works for square matrices only. // // This method can also be used to tag matrices as non-diagonal. */ Matrix& tagAsDiagonal(bool b=true); /** // Fill the main diagonal and set all other elements to zero. // This operation tags the matrix as diagonal. */ template Matrix& diagFill(const Num2* data, unsigned dataLen); //@{ /** A simple inspection of matrix properties */ inline unsigned nRows() const {return nrows_;} inline unsigned nColumns() const {return ncols_;} inline unsigned long length() const {return static_cast(nrows_)*ncols_;} inline Numeric* data() const {return data_;} inline bool isSquare() const {return data_ && ncols_ == nrows_;} bool isSymmetric() const; bool isAntiSymmetric() const; bool isDiagonal() const; + bool isMainDiagonal() const; // Works for non-square matrices bool isZero() const; //@} /** // This method resets the object to an unintialized state // Can be applied in order to force the assignment operators to work. */ Matrix& uninitialize(); /** // Check if this matrix has the same number of rows and columns // as the other one. "false" will be returned in case any one of // the two matrices (or both) is not initialized. */ bool isCompatible(const Matrix& other) const; /** // This method changes the object dimensions. // All data is lost in the process. */ Matrix& resize(unsigned nrows, unsigned ncols); /** This method sets all elements to 0 */ Matrix& zeroOut(); /** // This method sets all diagonal elements to 0. // It can only be used with square matrices. */ Matrix& clearMainDiagonal(); /** // This method sets all non-diagonal elements to 0. // This operation is only valid for square matrices. // It tags the matrix as diagonal. */ Matrix& makeDiagonal(); /** This method sets all elements to the given value */ Matrix& constFill(Numeric c); /** // Method for transposing this matrix. Uses std::swap for square // matrices but might allocate memory from the heap if the matrix // is not square. */ Matrix& Tthis(); //@{ /** Compare two matrices for equality */ bool operator==(const Matrix&) const; bool operator!=(const Matrix&) const; //@} /** // Non-const access to the data (works like 2-d array in C). // No bounds checking. */ Numeric* operator[](unsigned); /** // Const access to the data (works like 2-d array in C). // No bounds checking. */ const Numeric* operator[](unsigned) const; /** Data modification method with bounds checking */ Matrix& set(unsigned row, unsigned column, Numeric value); /** Access by value without bounds checking */ Numeric operator()(unsigned row, unsigned column) const; /** Access by value with bounds checking */ Numeric at(unsigned row, unsigned column) const; /** Sum of the values in the given row */ Numeric rowSum(unsigned row) const; /** Sum of the values in the given column */ Numeric columnSum(unsigned column) const; /** Remove row and/or column. Indices out of range are ignored */ Matrix removeRowAndColumn(unsigned row, unsigned column) const; /** Number of non-zero elements */ unsigned nonZeros() const; /** // Replace every n x m square block in the matrix by its sum. // The number of rows must be divisible by n. The number of // columns must be divisible by m. */ void coarseSum(unsigned n, unsigned m, Matrix* result) const; /** // Replace every n x m square block in the matrix by its average // The number of rows must be divisible by n. The number of // columns must be divisible by m. */ void coarseAverage(unsigned n, unsigned m, Matrix* result) const; /** Unary plus */ Matrix operator+() const; /** Unary minus */ Matrix operator-() const; //@{ /** Binary algebraic operation with matrix or scalar */ Matrix operator*(const Matrix& r) const; Matrix operator*(Numeric r) const; Matrix operator/(Numeric r) const; Matrix operator+(const Matrix& r) const; Matrix operator-(const Matrix& r) const; //@} /** Hadamard product (a.k.a. Schur product) */ Matrix hadamardProduct(const Matrix& r) const; /** Hadamard ratio. All elements of the denominator must be non-zero */ Matrix hadamardRatio(const Matrix& denominator) const; //@{ /** // Binary algebraic operation which writes its result into // an existing matrix potentially avoiding memory allocation */ void times(const Matrix& r, Matrix* result) const; void times(Numeric r, Matrix* result) const; void over(Numeric r, Matrix* result) const; void plus(const Matrix& r, Matrix* result) const; void minus(const Matrix& r, Matrix* result) const; //@} //@{ /** In-place algebraic operations with matrices and scalars */ Matrix& operator*=(Numeric r); Matrix& operator/=(Numeric r); Matrix& operator+=(const Matrix& r); Matrix& operator-=(const Matrix& r); //@} //@{ /** Multiplication by a vector represented by a pointer and array size */ template Matrix timesVector(const Num2* data, unsigned dataLen) const; template void timesVector(const Num2* data, unsigned dataLen, Num3* result, unsigned resultLen) const; //@} /** // Multiplication by a vector represented by a pointer and array size. // This function is useful when only one element of the result is needed // (or elements are needed one at a time). */ template Numeric timesVector(unsigned rowNumber, const Num2* data, unsigned dataLen) const; //@{ /** Multiplication by a row on the left */ template Matrix rowMultiply(const Num2* data, unsigned dataLen) const; template void rowMultiply(const Num2* data, unsigned dataLen, Num3* result, unsigned resultLen) const; //@} /** // Multiplication by a row on the left. This function is useful when // only one element of the result is needed (or elements are needed // one at a time). */ template Numeric rowMultiply(unsigned columnNumber, const Num2* data, unsigned dataLen) const; /** Bilinear form with a vector (v^T M v) */ template Numeric bilinear(const Num2* data, unsigned dataLen) const; /** Bilinear form with a matrix (r^T M r) */ Matrix bilinear(const Matrix& r) const; /** Bilinear form with a transposed matrix (r M r^T) */ Matrix bilinearT(const Matrix& r) const; /** // Solution of a linear system of equations M*x = rhs. // The matrix must be square. "false" is returned in case // the matrix is degenerate. */ bool solveLinearSystem(const Numeric* rhs, unsigned lenRhs, Numeric* solution) const; /** // Solution of linear systems of equations M*X = RHS. // The matrix M (this one) must be square. "false" is // returned in case this matrix is degenerate. */ bool solveLinearSystems(const Matrix& RHS, Matrix* X) const; /** // Solution of an underdetermined linear system of equations M*x = rhs // in the minimum norm sense (M is this matrix). The norm is defined by // sqrt(x^T*V^(-1)*x). The result is given by A*rhs, where A is the // right inverse of M (i.e., M*A = I) minimizing the norm. V plays // the role of the inverse Gram matrix. It must be symmetric and // positive-definite. // // This method can also be used for solving constrained least squares // problems. For example, minimization of (y - y0)^T*V^(-1)*(y - y0) // under the condition M*y = rhs can be reduced to solving an // underdetermined linear system using x = y - y0 and // M*x = rhs' = rhs - M*y0. The application code can call this function // assuming that V is the covariance matrix of y and using rhs' as the // "rhs" argument. Subsequently, y can be calculated as x + y0. // The covariance matrix of y can be subsequently calculated by // V_y = (I - P)*V*(I - P)^T, where P = A*M. Operator (I - P) actually // projects the x (i.e., y - y0) space onto the space of solutions of // the equation M*x = 0. // // "false" is returned in case of failure. */ bool underdeterminedLinearSystem(const Numeric* rhs, unsigned lenRhs, const Matrix& V, Numeric* solution, unsigned lenSolution, Numeric* resultNormSquared = 0, Matrix* A = 0) const; /** // Solution of a linear overdetermined system of equations M*x = rhs // in the least squares sense. The "lenRhs" parameter should be equal // to the number of matrix rows, and it should exceed "lenSolution" // ("lenSolution" should be equal to the number of columns). "false" // is returned in case of failure. */ bool linearLeastSquares(const Numeric* rhs, unsigned lenRhs, Numeric* solution, unsigned lenSolution) const; /** // Solution of a linear system of equations M*x = rhs1 in the least // squares sense, subject to the constraint B*x = rhs2. Number of // equations (i.e., dimensionality of rhs1) plus the number of // constraints (i.e., dimensionality of rhs2) should exceed the // dimensionality of x. "false" is returned in case of failure. If // the chi-square of the result is desired, make the "resultChiSquare" // argument point to some valid location. Same goes for the result // covariance matrix and all subsequent arguments which can be used // to extract various intermediate results: // // unconstrainedSolution -- Array to store the unconstrained solution. // Must have at least "lenSol" elements. // // unconstrainedCovmat -- Covariance matrix of the unconstrained // solution. // // projectionMatrix -- Projection matrix used to obtain a part // of the constrained solution from the // unconstrained one. // // A -- Matrix multiplied by rhs2 to get the second // part of the constrained solution. // // Note that requesting the result covariance matrix or any other // subsequent output switches the code to a different and slower // "textbook" algorithm instead of an optimized algorithm from Lapack. */ bool constrainedLeastSquares(const Numeric* rhs1, unsigned lenRhs1, const Matrix& B, const Numeric* rhs2, unsigned lenRhs2, Numeric* solution, unsigned lenSol, Numeric* resultChiSquare = 0, Matrix* resultCovarianceMatrix = 0, Numeric* unconstrainedSolution = 0, Matrix* unconstrainedCovmat = 0, Matrix* projectionMatrix = 0, Matrix* A = 0) const; /** // Weighted least squares problem. The inverse covariance matrix // must be symmetric and positive definite. If the chi-square of the // result is desired, make the "resultChiSquare" argument point to // some valid location. Same goes for the result covariance matrix. // "false" is returned in case of failure. */ bool weightedLeastSquares(const Numeric* rhs, unsigned lenRhs, const Matrix& inverseCovarianceMatrix, Numeric* solution, unsigned lenSolution, Numeric* resultChiSquare = 0, Matrix* resultCovarianceMatrix = 0) const; /** Return transposed matrix */ Matrix T() const; /** Return the product of this matrix transposed with this, M^T*M */ Matrix TtimesThis() const; /** Return the product of this matrix with its transpose, M*M^T */ Matrix timesT() const; /** Return the product of this matrix with the transposed argument */ Matrix timesT(const Matrix& r) const; /** Return the product of this matrix transposed with the argument */ Matrix Ttimes(const Matrix& r) const; /** // "directSum" method constructs a block-diagonal matrix with // this matrix and the argument matrix inside the diagonal blocks. // Block which contains this matrix is the top left one. */ Matrix directSum(const Matrix& added) const; /** Construct a symmetric matrix from this one (symmetrize) */ Matrix symmetrize() const; /** Construct an antisymmetric matrix from this one (antisymmetrize) */ Matrix antiSymmetrize() const; /** // Matrix outer product. The number of rows of the result equals // the product of the numbers of rows of this matrix and the argument // matrix. The number of columns of the result is the product of the // numbers of columns. */ Matrix outer(const Matrix& r) const; /** Minimum value for any row/column */ Numeric minValue() const; /** Maximum value for any row/column */ Numeric maxValue() const; /** Maximum absolute value for any row/column */ Numeric maxAbsValue() const; /** Row and column of the smallest matrix element */ std::pair argmin() const; /** Row and column of the largest matrix element */ std::pair argmax() const; /** // Row and column of the matrix element with the largest // absolute value */ std::pair argmaxAbs() const; /** Frobenius norm of this matrix */ double frobeniusNorm() const; //@{ /** Calculate the trace (matrix must be square) */ Numeric tr() const; inline Numeric sp() const {return tr();} //@} //@{ /** // Calculate the trace of a product of this matrix with another. // This works faster if only the trace is of interest but not the // product itself. */ Numeric productTr(const Matrix& r) const; inline Numeric productSp(const Matrix& r) const {return productTr(r);} //@} /** Calculate the determinant (matrix must be square) */ Numeric det() const; /** // Inverse of a real symmetric positive-definite matrix. Use this // function (it has better numerical precision than symInv()) if // you know that your matrix is symmetric and positive definite. // For example, a non-singular covariance matrix calculated for // some dataset has these properties. */ Matrix symPDInv() const; /** // Invert real symmetric positive-definite matrix using // eigendecomposition. Use this function if you know that your // matrix is supposed to be symmetric and positive-definite but // might not be due to numerical round-offs. The arguments of // this method are as folows: // // tol -- Eigenvalues whose ratios to the largest eigenvalue // are below this parameter will be either truncated // or extended (per value of the next argument). // This parameter must lie in (0, 1). // // extend -- If "true", the inverse of small eigenvalues will // be replaced by the inverse of the largest eigenvalue // divided by the tolerance parameter. If "false", // inverse of such eigenvalues will be set to 0. // // m -- LAPACK method to use for calculating the // eigendecomposition. */ Matrix symPDEigenInv(double tol, bool extend=true, EigenMethod m=EIGEN_SIMPLE) const; /** // Inverse of a real symmetric matrix. If the matrix is not symmetric, // its symmetrized version will be used internally. Use this function // (it has better numerical precision than inv()) if you know that // your matrix is symmetric. */ Matrix symInv() const; /** Inverse of a square tridiagonal matrix */ Matrix triDiagInv() const; /** // Inverse of a general matrix (which must be square and non-singular) */ Matrix inv() const; /** // Left pseudoinverse of a real, full rank matrix. This // pseudoinverse can be calculated when the number of columns // does not exceed the number of rows. */ Matrix leftInv() const; /** // Right pseudoinverse of a real, full rank matrix. This // pseudoinverse can be calculated when the number of rows // does not exceed the number of columns. */ Matrix rightInv() const; /** // Eigenvalues and eigenvectors of a real tridiagonal symmetric // matrix. The "eigenvectors" pointer can be NULL if only eigenvalues // are needed. If it is not NULL, the _rows_ of that matrix will be // set to the computed eigenvectors. // // The "m" parameter specifies the LAPACK method to use for // calculating the eigendecomposition. */ void tdSymEigen(Numeric* eigenvalues, unsigned lenEigenvalues, Matrix* eigenvectors=0, EigenMethod m=EIGEN_RRR) const; /** // Eigenvalues and eigenvectors of a real symmetric matrix. // The "eigenvectors" pointer can be NULL if only eigenvalues // are needed. If it is not NULL, the _rows_ of that matrix // will be set to the computed eigenvectors. // // The "m" parameter specifies the LAPACK method to use for // calculating the eigendecomposition. */ void symEigen(Numeric* eigenvalues, unsigned lenEigenvalues, Matrix* eigenvectors=0, EigenMethod m=EIGEN_SIMPLE) const; /** // Eigenvalues and eigenvectors of a real general matrix. // Either "rightEigenvectors" or "leftEigenvectors" or both // can be NULL if they are not needed. // // If an eigenvalue is real, the corresponding row of the // output eigenvector matrix is set to the computed eigenvector. // If an eigenvalue is complex, all such eigenvalues come in // complex conjugate pairs. Corresponding eigenvectors make up // such a pair as well. The matrix row which corresponds to the // first eigenvalue in the pair is set to the real part of the pair // and the matrix row which corresponds to the second eigenvalue // in the pair is set to the imaginary part of the pair. For the // right eigenvectors, the first complex eigenvector in the pair // needs to be constructed by adding the imaginary part. The second // eigenvector has the imaginary part subtracted. For the left // eigenvectors the order is reversed: the first eigenvector needs // to be constructed by subtracting the imaginary part, and the // second by adding it. */ void genEigen(typename GeneralizedComplex::type *eigenvalues, unsigned lenEigenvalues, Matrix* rightEigenvectors = 0, Matrix* leftEigenvectors = 0) const; /** // Singular value decomposition of a matrix. The length of the // singular values buffer must be at least min(nrows, ncolumns). // // On exit, matrices U and V will contain singular vectors of // the current matrix, row-by-row. The decomposition is then // // *this = U^T * diag(singularValues, nrows, ncolumns) * V */ void svd(Numeric* singularValues, unsigned lenSingularValues, Matrix* U, Matrix* V, SvdMethod m = SVD_D_AND_C) const; /** // Calculate entropy-based effective rank. It is assumed that // the matrix is symmetric positive semidefinite. // // The "tol" parameter specifies the eigenvalue cutoff: // all eigenvalues equal to tol*maxEigenvalue or smaller // will be ignored. This parameter must be non-negative. // // The "m" parameter specifies the LAPACK method to use for // calculating eigenvalues. // // If "eigenSum" pointer is not NULL, *eigenSum will be filled // on exit with the sum of eigenvalues above the cutoff divided // by the largest eigenvalue. */ double symPSDefEffectiveRank(double tol = 0.0, EigenMethod m = EIGEN_D_AND_C, double* eigenSum = 0) const; /** // Calculate a function of this matrix. It is assumed that the // matrix is symmetric and real. The calculation is performed by // calculating the eignenbasis, evaluating the function of the // eigenvalues, and then transforming back to the original basis. */ template Matrix symFcn(const Functor& fcn) const; /** Power of a matrix. Matrix must be square. */ Matrix pow(unsigned degree) const; /** // The following function derives a correlation matrix // out of a covariance matrix. In the returned matrix // all diagonal elements will be set to 1. */ Matrix covarToCorr() const; /** // The following function derives a covariance matrix // out of a correlation matrix and of the covariances. */ Matrix corrToCovar(const double* variances, unsigned lenVariances) const; //@{ /** Method related to "geners" I/O */ inline gs::ClassId classId() const {return gs::ClassId(*this);} bool write(std::ostream& of) const; //@} static const char* classname(); static inline unsigned version() {return 1;} static void restore(const gs::ClassId& id, std::istream& in, Matrix* m); private: void initialize(unsigned nrows, unsigned ncols); void calcDiagonal(); void makeCoarse(unsigned n, unsigned m, Matrix* result) const; void invertDiagonal(Matrix* result) const; Numeric local_[Len]; Numeric* data_; unsigned nrows_; unsigned ncols_; unsigned len_; bool isDiagonal_; bool diagonalityKnown_; #ifdef SWIG public: inline std::vector mainDiagonal() const { if (!nrows_) throw std::runtime_error( "In npstat::Matrix::mainDiagonal: uninitialized matrix"); const unsigned len = std::min(nrows_, ncols_); std::vector d(len); for (unsigned i=0; isetData(data, dataLen);} inline void setRow2(unsigned row, const Numeric* data, unsigned dataLen) {this->setRow(row, data, dataLen);} inline void setColumn2(unsigned col, const Numeric* data, unsigned dataLen) {this->setColumn(col, data, dataLen);} inline Matrix timesVector2(const Numeric* data, unsigned dataLen) const {return timesVector(data, dataLen);} inline Matrix rowMultiply2(const Numeric* data, unsigned dataLen) const {return rowMultiply(data, dataLen);} inline Numeric bilinear2(const Numeric* data, unsigned dataLen) const {return bilinear(data, dataLen);} inline std::vector solveLinearSystem2( const Numeric* data, unsigned dataLen) const { std::vector result(dataLen); if (!solveLinearSystem(data, dataLen, dataLen ? &result[0] : 0)) result.clear(); return result; } inline std::pair symPSDefEffectiveRank2( double tol = 0.0, EigenMethod m = EIGEN_D_AND_C) const { double eigenSum = 0; double erunk = this->symPSDefEffectiveRank(tol, m, &eigenSum); return std::pair(erunk, eigenSum); } inline std::vector symEigenValues() const { std::vector result(nrows_); symEigen(&result[0], nrows_, 0); return result; } inline Matrix rmul2(Numeric r) const {return *this * r;} // Perhaps, some day SWIG will be able to wrap the following. Currently, // its type deduction for GeneralizedComplex::type fails. // // inline std::vector::type> // genEigenValues() const // { // std::vector::type> r(nrows_); // genEigen(&r[0], nrows_); // return r; // } #endif // SWIG }; /** Utility for making square diagonal matrices from the given array */ template Matrix diag(const Numeric* data, unsigned dataLen); /** // Utility for making rectangular diagonal matrices from the given array. // The length of the array must be at least min(nrows, ncols). */ template Matrix diag(const Numeric* data, unsigned nrows, unsigned ncols); } #include template std::ostream& operator<<(std::ostream& os, const npstat::Matrix& m); template inline npstat::Matrix operator*(const Numeric& l, const npstat::Matrix& r) { return r*l; } #include "npstat/nm/Matrix.icc" #endif // NPSTAT_MATRIX_HH_ Index: trunk/npstat/nm/Matrix.icc =================================================================== --- trunk/npstat/nm/Matrix.icc (revision 709) +++ trunk/npstat/nm/Matrix.icc (revision 710) @@ -1,3061 +1,3076 @@ #include #include #include #include #include #include #include "geners/GenericIO.hh" #include "geners/IOException.hh" #include "npstat/nm/allocators.hh" #include "npstat/nm/lapack_interface.hh" #include "npstat/nm/PreciseType.hh" #include "npstat/nm/ComplexComparesAbs.hh" namespace npstat { template inline Matrix::Matrix() : data_(0), nrows_(0), ncols_(0), len_(0), isDiagonal_(false), diagonalityKnown_(false) { } template inline Matrix::Matrix(const Matrix& r) : data_(0), nrows_(r.nrows_), ncols_(r.ncols_), len_(r.len_), isDiagonal_(r.isDiagonal_), diagonalityKnown_(r.diagonalityKnown_) { if (len_) { if (len_ <= Len) data_ = local_; else data_ = new Numeric[len_]; for (unsigned i=0; i inline Matrix::Matrix(Matrix&& r) : data_(0), nrows_(r.nrows_), ncols_(r.ncols_), len_(r.len_), isDiagonal_(r.isDiagonal_), diagonalityKnown_(r.diagonalityKnown_) { if (len_) { if (len_ <= Len) { data_ = local_; for (unsigned i=0; i template inline Matrix::Matrix(const Matrix& r) : data_(0), nrows_(r.nrows_), ncols_(r.ncols_), len_(r.len_), isDiagonal_(false), diagonalityKnown_(false) { if (len_) { if (len_ <= Len) data_ = local_; else data_ = new Numeric[len_]; copyBuffer(data_, r.data_, len_); } else { nrows_ = 0U; ncols_ = 0U; } } template template inline Matrix& Matrix::setData( const Num2* data, const unsigned dataLength) { if (!nrows_) throw std::runtime_error( "In npstat::Matrix::setData: uninitialized matrix"); if (dataLength != len_) throw std::invalid_argument( "In npstat::Matrix::setData: incompatible input length"); assert(data); copyBuffer(data_, data, len_); diagonalityKnown_ = false; isDiagonal_ = false; return *this; } template template inline Matrix& Matrix::setRow( const unsigned row, const Num2* data, const unsigned dataLength) { if (!nrows_) throw std::runtime_error( "In npstat::Matrix::setRow: uninitialized matrix"); if (row >= nrows_) throw std::invalid_argument( "In npstat::Matrix::setRow: row number out of range"); if (dataLength != ncols_) throw std::invalid_argument( "In npstat::Matrix::setRow: incompatible input length"); assert(data); copyBuffer(data_+row*ncols_, data, len_); diagonalityKnown_ = false; isDiagonal_ = false; return *this; } template template inline Matrix& Matrix::setColumn( const unsigned col, const Num2* data, const unsigned dataLength) { if (!nrows_) throw std::runtime_error( "In npstat::Matrix::setColumn: uninitialized matrix"); if (col >= ncols_) throw std::invalid_argument( "In npstat::Matrix::setColumn: column number out of range"); if (dataLength != nrows_) throw std::invalid_argument( "In npstat::Matrix::setColumn: incompatible input length"); assert(data); Numeric* d = data_ + col; for (unsigned i=0; i template inline Matrix& Matrix::setFromTriplets( Iterator first, Iterator const last) { zeroOut(); for (; first != last; ++first) { const unsigned r = static_cast(first->row()); if (r >= nrows_) throw std::invalid_argument( "In npstat::Matrix::setFromTriplets: row number out of range"); const unsigned c = static_cast(first->col()); if (c >= ncols_) throw std::invalid_argument( "In npstat::Matrix::setFromTriplets: column number out of range"); data_[r*ncols_ + c] += first->value(); } return *this; } template template inline Matrix::Matrix( const Matrix& r, const unsigned rowMin, const unsigned rowMax, const unsigned columnMin, const unsigned columnMax) : data_(0), nrows_(rowMax - rowMin), ncols_(columnMax - columnMin), len_(nrows_*ncols_), isDiagonal_(false), diagonalityKnown_(false) { if (rowMin > rowMax || columnMin > columnMax) throw std::invalid_argument( "In npstat::Matrix subrange constructor: " "invalid subrange specification"); if (len_) { if (rowMax > r.nrows_) throw std::invalid_argument( "In npstat::Matrix subrange constructor: " "maximum row number out of range"); if (columnMax > r.ncols_) throw std::invalid_argument( "In npstat::Matrix subrange constructor: " "maximum column number out of range"); if (len_ <= Len) data_ = local_; else data_ = new Numeric[len_]; for (unsigned row=0; row inline void Matrix::initialize(const unsigned nrows, const unsigned ncols) { const unsigned long sz = static_cast(nrows)*ncols; if (sz > static_cast(UINT_MAX)) throw std::invalid_argument( "In npstat::Matrix::validateRowsCols: requested matrix size is too large"); if (!sz) throw std::invalid_argument( "In npstat::Matrix::validateRowsCols: both number of " "rows and number of columns must be positive"); len_ = sz; if (len_ <= Len) data_ = local_; else data_ = new Numeric[len_]; } template inline Matrix::Matrix(const unsigned nrows, const unsigned ncols) : data_(0), nrows_(nrows), ncols_(ncols), len_(0), isDiagonal_(false), diagonalityKnown_(false) { initialize(nrows, ncols); } template inline Matrix::Matrix(const unsigned nrows, const unsigned ncols, const Numeric* data) : data_(0), nrows_(nrows), ncols_(ncols), len_(0), isDiagonal_(false), diagonalityKnown_(false) { initialize(nrows, ncols); assert(data); for (unsigned i=0; i inline Matrix& Matrix::zeroOut() { const Numeric zero = Numeric(); for (unsigned i=0; i inline Matrix& Matrix::constFill(const Numeric c) { for (unsigned i=0; i inline Matrix::Matrix(const unsigned nrows, const unsigned ncols, const int initCode) : data_(0), nrows_(nrows), ncols_(ncols), len_(0), isDiagonal_(false), diagonalityKnown_(false) { initialize(nrows, ncols); zeroOut(); switch (initCode) { case 0: break; case 1: { if (nrows_ != ncols_) { if (len_ > Len) delete [] data_; throw std::invalid_argument( "In npstat::Matrix constructor: " "unit matrix must be square"); } const Numeric one(static_cast(1)); for (unsigned i=0; i Len) delete [] data_; throw std::invalid_argument( "In npstat::Matrix constructor: " "invalid initialization code"); } } } template inline Matrix::~Matrix() { if (data_ != local_) delete [] data_; } template inline Matrix& Matrix::uninitialize() { if (data_ != local_) delete [] data_; data_ = 0; nrows_ = 0; ncols_ = 0; len_ = 0; diagonalityKnown_ = false; isDiagonal_ = false; return *this; } template inline void Matrix::invertDiagonal(Matrix* result) const { assert(result != this); Numeric* rdata = result->data_; const unsigned len = nrows_*ncols_; const Numeric zero = Numeric(); const Numeric one = static_cast(1); for (unsigned i=0; idiagonalityKnown_ = true; result->isDiagonal_ = true; } template inline Matrix Matrix::symPDInv() const { if (!nrows_) throw std::runtime_error( "In npstat::Matrix::symPDInv: uninitialized matrix"); if (nrows_ != ncols_) throw std::domain_error( "In npstat::Matrix::symPDInv: matrix is not square"); Matrix result(nrows_, nrows_); if (isDiagonal()) invertDiagonal(&result); else invert_posdef_sym_matrix(data_, nrows_, result.data_); return result; } template inline Matrix Matrix::symInv() const { if (!nrows_) throw std::runtime_error( "In npstat::Matrix::symInv: uninitialized matrix"); if (nrows_ != ncols_) throw std::domain_error( "In npstat::Matrix::symInv: matrix is not square"); Matrix result(nrows_, nrows_); if (isDiagonal()) invertDiagonal(&result); else invert_sym_matrix(data_, nrows_, result.data_); return result; } template inline Matrix Matrix::inv() const { if (!nrows_) throw std::runtime_error( "In npstat::Matrix::inv: uninitialized matrix"); if (nrows_ != ncols_) throw std::domain_error( "In npstat::Matrix::inv: matrix is not square"); Matrix result(nrows_, nrows_); if (isDiagonal()) invertDiagonal(&result); else invert_general_matrix(data_, nrows_, result.data_); return result; } template inline Matrix Matrix::triDiagInv() const { if (!nrows_) throw std::runtime_error( "In npstat::Matrix::triDiagInv: uninitialized matrix"); if (nrows_ != ncols_) throw std::domain_error( "In npstat::Matrix::triDiagInv: matrix is not square"); Matrix result(nrows_, nrows_); if (isDiagonal()) invertDiagonal(&result); else invert_td_matrix(data_, nrows_, result.data_); return result; } template inline Matrix Matrix::leftInv() const { if (!nrows_) throw std::runtime_error( "In npstat::Matrix::leftInv: uninitialized matrix"); if (nrows_ < ncols_) throw std::domain_error( "In npstat::Matrix::leftInv: more columns than rows"); return TtimesThis().symInv().timesT(*this); } template inline Matrix Matrix::rightInv() const { if (!nrows_) throw std::runtime_error( "In npstat::Matrix::rightInv: uninitialized matrix"); if (nrows_ > ncols_) throw std::domain_error( "In npstat::Matrix::rightInv: more rows than columns"); return Ttimes(timesT().symInv()); } template inline void Matrix::genEigen( typename GeneralizedComplex::type *eigvals, const unsigned lenEigvals, Matrix* rightEigenvectors, Matrix* leftEigenvectors) const { if (!nrows_) throw std::runtime_error( "In npstat::Matrix::genEigen: uninitialized matrix"); if (nrows_ != ncols_) throw std::domain_error( "In npstat::Matrix::genEigen: matrix is not square"); if (lenEigvals < ncols_) throw std::invalid_argument( "In npstat::Matrix::genEigen: " "insufficient length of the eigenvalue buffer"); assert(eigvals); Numeric* right = 0; if (rightEigenvectors) { if (rightEigenvectors == this) throw std::invalid_argument( "In npstat::Matrix::genEigen: " "can not replace matrix with its own right eigenvectors"); if (rightEigenvectors->data_ == 0) *rightEigenvectors = *this; if (!(nrows_ == rightEigenvectors->nrows_ && nrows_ == rightEigenvectors->ncols_)) throw std::invalid_argument( "In npstat::Matrix::genEigen: " "incompatible matrix for storing right eigenvectors"); right = rightEigenvectors->data_; rightEigenvectors->isDiagonal_ = false; rightEigenvectors->diagonalityKnown_ = false; } Numeric* left = 0; if (leftEigenvectors) { if (leftEigenvectors == this) throw std::invalid_argument( "In npstat::Matrix::genEigen: " "can not replace matrix with its own left eigenvectors"); if (leftEigenvectors->data_ == 0) *leftEigenvectors = *this; if (!(nrows_ == leftEigenvectors->nrows_ && nrows_ == leftEigenvectors->ncols_)) throw std::invalid_argument( "In npstat::Matrix::genEigen: " "incompatible matrix for storing left eigenvectors"); left = leftEigenvectors->data_; leftEigenvectors->isDiagonal_ = false; leftEigenvectors->diagonalityKnown_ = false; } if (left || right) gen_matrix_eigensystem(data_, ncols_, eigvals, right, left); else gen_matrix_eigenvalues(data_, ncols_, eigvals); } template inline void Matrix::svd(Numeric* singularValues, const unsigned lenSingularValues, Matrix* U, Matrix* V, const SvdMethod m) const { if (!nrows_) throw std::runtime_error( "In npstat::Matrix::svd: uninitialized matrix"); if (lenSingularValues < std::min(nrows_, ncols_)) throw std::invalid_argument( "In npstat::Matrix::svd: " "insufficient length of the singular values buffer"); assert(singularValues); assert(U); assert(V); if (U == this || V == this) throw std::invalid_argument( "In npstat::Matrix::svd: " "can not replace matrix with its own decomposition"); U->resize(nrows_, nrows_); V->resize(ncols_, ncols_); bool status = false; switch (m) { case SVD_SIMPLE: status = gen_matrix_svd(data_, nrows_, ncols_, U->data_, singularValues, V->data_); break; case SVD_D_AND_C: status = gen_matrix_svd_dc(data_, nrows_, ncols_, U->data_, singularValues, V->data_); break; default: assert(!"In npstat::Matrix::svd: incomplete switch " "statement. This is a bug. Please report."); } if (!status) throw std::runtime_error( "In npstat::Matrix::svd: algorithm did not converge"); } template inline void Matrix::symEigen(Numeric* eigenvalues, const unsigned lenEigenvalues, Matrix* eigenvectors, const EigenMethod m) const { if (!nrows_) throw std::runtime_error( "In npstat::Matrix::symEigen: uninitialized matrix"); if (nrows_ != ncols_) throw std::domain_error( "In npstat::Matrix::symEigen: matrix is not square"); if (lenEigenvalues < ncols_) throw std::invalid_argument( "In npstat::Matrix::symEigen: " "insufficient length of the eigenvalue buffer"); assert(eigenvalues); if (eigenvectors) { if (eigenvectors == this) throw std::invalid_argument( "In npstat::Matrix::symEigen: " "can not replace matrix with its own eigenvectors"); if (eigenvectors->data_ == 0) *eigenvectors = *this; if (!(nrows_ == eigenvectors->nrows_ && nrows_ == eigenvectors->ncols_)) throw std::invalid_argument( "In npstat::Matrix::symEigen: " "incompatible matrix for storing eigenvectors"); eigenvectors->isDiagonal_ = false; eigenvectors->diagonalityKnown_ = false; switch (m) { case EIGEN_SIMPLE: sym_matrix_eigensystem(data_, ncols_, eigenvalues, eigenvectors->data_); break; case EIGEN_D_AND_C: sym_matrix_eigensystem_dc(data_, ncols_, eigenvalues, eigenvectors->data_); break; case EIGEN_RRR: sym_matrix_eigensystem_rrr(data_, ncols_, eigenvalues, eigenvectors->data_); break; default: assert(!"In npstat::Matrix::symEigen: incomplete switch " "statement 1. This is a bug. Please report."); } } else switch (m) { case EIGEN_SIMPLE: sym_matrix_eigenvalues(data_, ncols_, eigenvalues); break; case EIGEN_D_AND_C: sym_matrix_eigenvalues_dc(data_, ncols_, eigenvalues); break; case EIGEN_RRR: sym_matrix_eigenvalues_rrr(data_, ncols_, eigenvalues); break; default: assert(!"In npstat::Matrix::symEigen: incomplete switch " "statement 2. This is a bug. Please report."); } } template inline void Matrix::tdSymEigen(Numeric* eigenvalues, const unsigned lenEigenvalues, Matrix* eigenvectors, const EigenMethod m) const { if (!nrows_) throw std::runtime_error( "In npstat::Matrix::tdSymEigen: uninitialized matrix"); if (nrows_ != ncols_) throw std::domain_error( "In npstat::Matrix::tdSymEigen: matrix is not square"); if (lenEigenvalues < ncols_) throw std::invalid_argument( "In npstat::Matrix::tdSymEigen: " "insufficient length of the eigenvalue buffer"); assert(eigenvalues); if (eigenvectors) { if (eigenvectors == this) throw std::invalid_argument( "In npstat::Matrix::tdSymEigen: " "can not replace matrix with its own eigenvectors"); if (eigenvectors->data_ == 0) *eigenvectors = *this; if (!(nrows_ == eigenvectors->nrows_ && nrows_ == eigenvectors->ncols_)) throw std::invalid_argument( "In npstat::Matrix::tdSymEigen: " "incompatible matrix for storing eigenvectors"); eigenvectors->isDiagonal_ = false; eigenvectors->diagonalityKnown_ = false; switch (m) { case EIGEN_SIMPLE: td_sym_matrix_eigensystem(data_, ncols_, eigenvalues, eigenvectors->data_); break; case EIGEN_D_AND_C: td_sym_matrix_eigensystem_dc(data_, ncols_, eigenvalues, eigenvectors->data_); break; case EIGEN_RRR: td_sym_matrix_eigensystem_rrr(data_, ncols_, eigenvalues, eigenvectors->data_); break; default: assert(!"In npstat::Matrix::tdSymEigen: incomplete switch " "statement 1. This is a bug. Please report."); } } else switch (m) { case EIGEN_SIMPLE: td_sym_matrix_eigenvalues(data_, ncols_, eigenvalues); break; case EIGEN_D_AND_C: td_sym_matrix_eigenvalues_dc(data_, ncols_, eigenvalues); break; case EIGEN_RRR: td_sym_matrix_eigenvalues_rrr(data_, ncols_, eigenvalues); break; default: assert(!"In npstat::Matrix::tdSymEigen: incomplete switch " "statement 2. This is a bug. Please report."); } } template inline Matrix Matrix::symPDEigenInv( const double tol, const bool extend, const EigenMethod eigenm) const { if (!nrows_) throw std::runtime_error( "In npstat::Matrix::symPDEigenInv: uninitialized matrix"); if (nrows_ != ncols_) throw std::domain_error( "In npstat::Matrix::symPDEigenInv: matrix is not square"); if (!(tol > 0.0 && tol < 1.0)) throw std::domain_error( "In npstat::Matrix::symPDEigenInv: invalid tolerance parameter"); Matrix evec; std::vector evBuf(nrows_); Numeric* eigenvalues = &evBuf[0]; this->symEigen(eigenvalues, nrows_, &evec, eigenm); const unsigned nm1 = nrows_ - 1U; const double emax = eigenvalues[nm1]; if (emax <= 0.0) throw std::domain_error( "In npstat::Matrix::symPDEigenInv: " "largest eigenvalue is not positive"); const Numeric ecut = emax*tol; const Numeric replacementValue = 1.0/(emax*tol); const Numeric one = 1.0; eigenvalues[nm1] = one/eigenvalues[nm1]; for (unsigned i=0; i ecut) eigenvalues[i] = one/eigenvalues[i]; else if (extend) eigenvalues[i] = replacementValue; else eigenvalues[i] = 0; } Matrix dmat(nrows_, nrows_); dmat.diagFill(eigenvalues, nrows_); return dmat.bilinear(evec); } template inline double Matrix::symPSDefEffectiveRank( const double tol, const EigenMethod m, double* eigenSum) const { if (!nrows_) throw std::runtime_error( "In npstat::Matrix::symPSDefEffectiveRank: uninitialized matrix"); if (nrows_ != ncols_) throw std::domain_error( "In npstat::Matrix::symPSDefEffectiveRank: matrix is not square"); if (tol < 0.0) throw std::domain_error( "In npstat::Matrix::symPSDefEffectiveRank: " "tolerance can not be negative"); std::vector eigenvaluesVec(nrows_); this->symEigen(&eigenvaluesVec[0], nrows_, 0, m); const Numeric* eigenvalues = &eigenvaluesVec[0]; if (eigenvalues[nrows_ - 1] <= static_cast(0)) { if (eigenSum) *eigenSum = 0.0; return 0.0; } const Numeric cutoff = tol*eigenvalues[nrows_ - 1]; unsigned ifirst = nrows_ - 1U; for (; ifirst; --ifirst) if (eigenvalues[ifirst - 1] <= cutoff) break; long double sum = 0.0L, entropy = 0.0L; Numeric largest = Numeric(); for (unsigned i=ifirst; i largest) largest = eigenvalues[i]; entropy -= eigenvalues[i]*std::log(eigenvalues[i]); } if (eigenSum) *eigenSum = sum/static_cast(largest); return sum*std::exp(entropy/sum); } template template inline Matrix Matrix::symFcn(const Functor& fcn) const { if (!nrows_) throw std::runtime_error( "In npstat::Matrix::symFcn: uninitialized matrix"); if (nrows_ != ncols_) throw std::domain_error( "In npstat::Matrix::symFcn: matrix is not square"); Matrix evectors(nrows_, nrows_); Matrix evalues(nrows_, nrows_); symEigen(evalues.data_, nrows_, &evectors); const Numeric zero = Numeric(); for (unsigned i=1; i inline Matrix Matrix::pow(const unsigned degree) const { if (!nrows_) throw std::runtime_error( "In npstat::Matrix::pow: uninitialized matrix"); if (nrows_ != ncols_) throw std::domain_error( "In npstat::Matrix::pow: matrix is not square"); if (degree == 0U) return Matrix(nrows_, ncols_, 1); if (degree == 1U) return *this; Matrix result(*this); for (unsigned i=1; i inline bool Matrix::isSymmetric() const { if (nrows_ == 0) return false; if (nrows_ != ncols_) return false; for (unsigned i=1; i inline void Matrix::calcDiagonal() { isDiagonal_ = false; if (len_) { const Numeric null = Numeric(); isDiagonal_ = (nrows_ == ncols_); for (unsigned i=0; i inline bool Matrix::isDiagonal() const { if (!diagonalityKnown_) (const_cast(this))->calcDiagonal(); return isDiagonal_; } template + inline bool Matrix::isMainDiagonal() const + { + if (!nrows_) + return false; + if (nrows_ == ncols_) + return isDiagonal(); + const Numeric null = Numeric(); + bool diagonal = true; + for (unsigned i=0; i inline bool Matrix::isZero() const { if (!nrows_) throw std::runtime_error( "In npstat::Matrix::isZero: uninitialized matrix"); const Numeric zero = Numeric(); for (unsigned i=0; i inline bool Matrix::isAntiSymmetric() const { if (nrows_ == 0) return false; if (nrows_ != ncols_) return false; for (unsigned i=1; i inline Matrix Matrix::directSum(const Matrix& r) const { if (!nrows_) throw std::runtime_error( "In npstat::Matrix::directSum: uninitialized matrix"); if (!r.nrows_) throw std::runtime_error( "In npstat::Matrix::directSum: uninitialized argument"); const unsigned newCols = ncols_+r.ncols_; Matrix result(nrows_+r.nrows_, newCols, 0); for (unsigned row=0; row inline Matrix Matrix::symmetrize() const { if (!nrows_) throw std::runtime_error( "In npstat::Matrix::symmetrize: uninitialized matrix"); if (nrows_ != ncols_) throw std::domain_error( "In npstat::Matrix::symmetrize: matrix is not square"); Matrix result(ncols_, nrows_); const Numeric two(static_cast(2)); for (unsigned i=0; i inline Matrix Matrix::antiSymmetrize() const { if (!nrows_) throw std::runtime_error( "In npstat::Matrix::antiSymmetrize: uninitialized matrix"); if (nrows_ != ncols_) throw std::domain_error( "In npstat::Matrix::antiSymmetrize: matrix is not square"); Matrix result(ncols_, nrows_); const Numeric two(static_cast(2)); const Numeric zero = Numeric(); for (unsigned i=0; i inline Matrix& Matrix::resize( const unsigned nrows, const unsigned ncols) { if (nrows != nrows_ || ncols != ncols_) { uninitialize(); len_ = nrows*ncols; if (len_) { nrows_ = nrows; ncols_ = ncols; if (len_ <= Len) data_ = local_; else data_ = new Numeric[len_]; } } else { isDiagonal_ = false; diagonalityKnown_ = false; } return *this; } template inline Matrix& Matrix::tagAsDiagonal( const bool b) { if (!data_) throw std::runtime_error( "In npstat::Matrix::tagAsDiagonal: uninitialized matrix"); if (nrows_ != ncols_) throw std::domain_error( "In npstat::Matrix::tagAsDiagonal: matrix is not square"); diagonalityKnown_ = true; isDiagonal_ = b; return *this; } template inline Matrix& Matrix::clearMainDiagonal() { const Numeric null = Numeric(); if (nrows_ != ncols_) throw std::domain_error( "In npstat::Matrix::clearMainDiagonal: matrix is not square"); for (unsigned i=0; i inline Matrix& Matrix::Tthis() { if (!nrows_) throw std::runtime_error( "In npstat::Matrix::Tthis: uninitialized matrix"); if (!(diagonalityKnown_ && isDiagonal_)) { if (nrows_ == ncols_) { for (unsigned i=0; i inline Matrix& Matrix::makeDiagonal() { const Numeric null = Numeric(); tagAsDiagonal(); for (unsigned i=0; i template inline Matrix& Matrix::diagFill(const Num2* data, const unsigned dataLen) { if (dataLen) assert(data); if (data_) { if (nrows_ != ncols_) throw std::domain_error( "In npstat::Matrix::diagFill: matrix is not square"); if (nrows_ != dataLen) throw std::invalid_argument( "In npstat::Matrix::diagFill: " "incompatible size of the input array"); } else if (dataLen) { nrows_ = dataLen; ncols_ = dataLen; len_ = dataLen*dataLen; if (len_ <= Len) data_ = local_; else data_ = new Numeric[len_]; } if (dataLen) { if (!(diagonalityKnown_ && isDiagonal_)) zeroOut(); for (unsigned i=0; i(data[i]); diagonalityKnown_ = true; isDiagonal_ = true; } return *this; } template inline bool Matrix::isCompatible(const Matrix& r) const { if (!nrows_) return false; else return nrows_ == r.nrows_ && ncols_ == r.ncols_; } template inline Matrix& Matrix::operator=(const Matrix& r) { if (this == &r) return *this; if (data_) { if (!(nrows_ == r.nrows_ && ncols_ == r.ncols_)) throw std::domain_error( "In npstat::Matrix::operator=: " "incompatible argument dimensions"); } else { nrows_ = r.nrows_; ncols_ = r.ncols_; len_ = r.len_; if (len_ <= Len) data_ = local_; else data_ = new Numeric[len_]; } for (unsigned i=0; i inline Matrix& Matrix::operator=(Matrix&& r) { if (this == &r) return *this; if (data_) { if (!(nrows_ == r.nrows_ && ncols_ == r.ncols_)) throw std::domain_error( "In npstat::Matrix::operator=: " "incompatible argument dimensions"); if (data_ != local_) { delete [] data_; data_ = 0; } } nrows_ = r.nrows_; ncols_ = r.ncols_; len_ = r.len_; if (len_ <= Len) { data_ = local_; for (unsigned i=0; i template inline Matrix& Matrix::operator=(const Matrix& r) { if (reinterpret_cast(this) == reinterpret_cast(&r)) return *this; if (data_) { if (!(nrows_ == r.nrows_ && ncols_ == r.ncols_)) throw std::domain_error( "In npstat::Matrix::operator=: " "incompatible argument dimensions"); } else { nrows_ = r.nrows_; ncols_ = r.ncols_; len_ = r.len_; if (len_ <= Len) data_ = local_; else data_ = new Numeric[len_]; } copyBuffer(data_, r.data_, len_); if (r.diagonalityKnown_ && r.isDiagonal_) { diagonalityKnown_ = true; isDiagonal_ = true; } else { diagonalityKnown_ = false; isDiagonal_ = false; } return *this; } template inline bool Matrix::operator==(const Matrix& r) const { if (!(nrows_ && r.nrows_)) throw std::runtime_error( "In npstat::Matrix::operator==: uninitialized matrix"); if (!(nrows_ == r.nrows_ && ncols_ == r.ncols_)) throw std::domain_error( "In npstat::Matrix::operator==: " "incompatible argument dimensions"); for (unsigned i=0; i inline bool Matrix::operator!=(const Matrix& r) const { return !(*this == r); } template inline const Numeric* Matrix::operator[](const unsigned i) const { return data_ + i*ncols_; } template inline Numeric* Matrix::operator[](const unsigned i) { diagonalityKnown_ = false; return data_ + i*ncols_; } template inline Numeric Matrix::operator()( const unsigned row, const unsigned column) const { return data_[row*ncols_ + column]; } template inline Numeric Matrix::at( const unsigned row, const unsigned column) const { if (row >= nrows_) throw std::out_of_range("npstat::Matrix::at: " "row index out of range"); if (column >= ncols_) throw std::out_of_range("npstat::Matrix::at: " "column index out of range"); return data_[row*ncols_ + column]; } template inline Numeric Matrix::rowSum(const unsigned row) const { typedef typename PreciseType::type Precise; if (row >= nrows_) throw std::out_of_range("npstat::Matrix::rowSum: " "row index out of range"); Precise sum = Precise(); const Numeric* rData = data_ + row*ncols_; for (unsigned col=0; col(sum); } template inline Numeric Matrix::columnSum(const unsigned column) const { typedef typename PreciseType::type Precise; if (column >= ncols_) throw std::out_of_range("npstat::Matrix::columnSum: " "column index out of range"); Precise sum = Precise(); const Numeric* columnData = data_ + column; for (unsigned row=0; row(sum); } template inline Matrix& Matrix::set( const unsigned row, const unsigned column, const Numeric value) { if (row >= nrows_) throw std::out_of_range("npstat::Matrix::set: " "row index out of range"); if (column >= ncols_) throw std::out_of_range("npstat::Matrix::set: " "column index out of range"); data_[row*ncols_ + column] = value; diagonalityKnown_ = false; return *this; } template inline Matrix Matrix::operator+(const Matrix& r) const { if (!(nrows_ && r.nrows_)) throw std::runtime_error( "In npstat::Matrix::operator+: uninitialized matrix"); if (!(nrows_ == r.nrows_ && ncols_ == r.ncols_)) throw std::domain_error( "In npstat::Matrix::operator+: " "incompatible argument dimensions"); Matrix result(nrows_, ncols_); for (unsigned i=0; i inline Matrix Matrix::hadamardProduct(const Matrix& r) const { if (!(nrows_ && r.nrows_)) throw std::runtime_error( "In npstat::Matrix::hadamardProduct: uninitialized matrix"); if (!(nrows_ == r.nrows_ && ncols_ == r.ncols_)) throw std::domain_error( "In npstat::Matrix::hadamardProduct: " "incompatible argument dimensions"); Matrix result(nrows_, ncols_); for (unsigned i=0; i inline Matrix Matrix::hadamardRatio(const Matrix& r) const { if (!(nrows_ && r.nrows_)) throw std::runtime_error( "In npstat::Matrix::hadamardRatio: uninitialized matrix"); if (!(nrows_ == r.nrows_ && ncols_ == r.ncols_)) throw std::domain_error( "In npstat::Matrix::hadamardRatio: " "incompatible argument dimensions"); Matrix result(nrows_, ncols_); for (unsigned i=0; i(0)) throw std::domain_error( "In npstat::Matrix::hadamardRatio: division by zero"); result.data_[i] = data_[i]/r.data_[i]; } return result; } template inline void Matrix::plus(const Matrix& r, Matrix* result) const { if (!(nrows_ && r.nrows_)) throw std::runtime_error( "In npstat::Matrix::plus: uninitialized matrix"); if (!(nrows_ == r.nrows_ && ncols_ == r.ncols_)) throw std::domain_error( "In npstat::Matrix::plus: " "incompatible argument dimensions"); assert(result); if (result != this) result->resize(nrows_, ncols_); Numeric* const rdata = result->data_; for (unsigned i=0; i inline Matrix Matrix::operator-(const Matrix& r) const { if (!(nrows_ && r.nrows_)) throw std::runtime_error( "In npstat::Matrix::operator-: uninitialized matrix"); if (!(nrows_ == r.nrows_ && ncols_ == r.ncols_)) throw std::domain_error( "In npstat::Matrix::operator-: " "incompatible argument dimensions"); Matrix result(nrows_, ncols_); for (unsigned i=0; i inline void Matrix::minus(const Matrix& r, Matrix* result) const { if (!(nrows_ && r.nrows_)) throw std::runtime_error( "In npstat::Matrix::minus: uninitialized matrix"); if (!(nrows_ == r.nrows_ && ncols_ == r.ncols_)) throw std::domain_error( "In npstat::Matrix::minus: " "incompatible argument dimensions"); assert(result); result->resize(nrows_, ncols_); Numeric* const rdata = result->data_; for (unsigned i=0; i inline Matrix Matrix::T() const { if (!nrows_) throw std::runtime_error( "In npstat::Matrix::T: uninitialized matrix"); Matrix result(ncols_, nrows_); if (diagonalityKnown_ && isDiagonal_) copyBuffer(result.data_, data_, nrows_*ncols_); else transposeBuffer(result.data_, data_, nrows_, ncols_); result.isDiagonal_ = isDiagonal_; result.diagonalityKnown_ = diagonalityKnown_; return result; } template inline Matrix Matrix::TtimesThis() const { if (!nrows_) throw std::runtime_error( "In npstat::Matrix::TtimesThis: uninitialized matrix"); Matrix result(ncols_, ncols_); if (isDiagonal()) { result.zeroOut(); for (unsigned row=0; row inline Matrix Matrix::bilinearT(const Matrix& r) const { if (!(nrows_ && r.nrows_)) throw std::runtime_error( "In npstat::Matrix::bilinearT: uninitialized matrix"); if (nrows_ != ncols_) throw std::domain_error( "In npstat::Matrix::bilinearT: matrix is not square"); if (ncols_ != r.ncols_) throw std::domain_error( "In npstat::Matrix::bilinearT: incompatible matrix dimensions"); Matrix result(r.nrows_, r.nrows_); const bool thisDiag = isDiagonal(); const bool rDiag = r.isDiagonal(); if (thisDiag && rDiag) { result.zeroOut(); for (unsigned row=0; row Len) buf = new Numeric[nrows_]; for (unsigned col=0; col inline Matrix Matrix::bilinear(const Matrix& r) const { if (!(nrows_ && r.nrows_)) throw std::runtime_error( "In npstat::Matrix::bilinear: uninitialized matrix"); if (nrows_ != ncols_) throw std::domain_error( "In npstat::Matrix::bilinear: matrix is not square"); if (ncols_ != r.nrows_) throw std::domain_error( "In npstat::Matrix::bilinear: incompatible matrix dimensions"); Matrix result(r.ncols_, r.ncols_); const bool thisDiag = isDiagonal(); const bool rDiag = r.isDiagonal(); if (thisDiag && rDiag) { result.zeroOut(); for (unsigned row=0; row Len) buf = new Numeric[nrows_]; for (unsigned col=0; col inline Matrix Matrix::Ttimes(const Matrix& r) const { if (!(nrows_ && r.nrows_)) throw std::runtime_error( "In npstat::Matrix::Ttimes: uninitialized matrix"); if (nrows_ != r.nrows_) throw std::domain_error( "In npstat::Matrix::Ttimes: incompatible matrix dimensions"); Matrix result(ncols_, r.ncols_); const bool thisDiag = isDiagonal(); const bool rDiag = r.isDiagonal(); if (thisDiag && rDiag) { result.zeroOut(); for (unsigned row=0; row inline Matrix Matrix::timesT(const Matrix& r) const { if (!(nrows_ && r.nrows_)) throw std::runtime_error( "In npstat::Matrix::timesT: uninitialized matrix"); if (ncols_ != r.ncols_) throw std::domain_error( "In npstat::Matrix::timesT: incompatible matrix dimensions"); Matrix result(nrows_, r.nrows_); const bool thisDiag = isDiagonal(); const bool rDiag = r.isDiagonal(); if (thisDiag && rDiag) { result.zeroOut(); for (unsigned row=0; row inline Matrix Matrix::timesT() const { if (!nrows_) throw std::runtime_error( "In npstat::Matrix::timesT: uninitialized matrix"); Matrix result(nrows_, nrows_); if (isDiagonal()) { result.zeroOut(); for (unsigned row=0; row inline double Matrix::frobeniusNorm() const { if (!nrows_) throw std::runtime_error( "In npstat::Matrix::frobeniusNorm: uninitialized matrix"); long double sum = 0.0L; for (unsigned i=0; i inline Matrix Matrix::covarToCorr() const { if (!nrows_) throw std::runtime_error( "In npstat::Matrix::covarToCorr: uninitialized matrix"); if (nrows_ != ncols_) throw std::domain_error( "In npstat::Matrix::covarToCorr: matrix is not square"); Matrix result(ncols_, nrows_); Numeric* rdata = result.data_; for (unsigned row=0; row(1); return result; } template inline Matrix Matrix::corrToCovar( const double* variances, const unsigned lenVariances) const { if (!nrows_) throw std::runtime_error( "In npstat::Matrix::corrToCovar: uninitialized matrix"); if (nrows_ != ncols_) throw std::domain_error( "In npstat::Matrix::corrToCovar: matrix is not square"); if (lenVariances != nrows_) throw std::invalid_argument( "In npstat::Matrix::corrToCovar: " "incompatible length of covariance array"); assert(variances); for (unsigned row=0; row inline void Matrix::makeCoarse( const unsigned n, const unsigned m, Matrix* result) const { assert(result); if (result == this) throw std::invalid_argument( "In npstat::Matrix::makeCoarse: " "can not replace matrix with its own coarse version"); const unsigned newRows = nrows_ / n; const unsigned newCols = ncols_ / m; result->resize(newRows, newCols); if (m == 1U && n == 1U) { copyBuffer(result->data_, data_, len_); result->isDiagonal_ = isDiagonal_; result->diagonalityKnown_ = diagonalityKnown_; } else if (m == 1U) { // Sum the rows for (unsigned row=0; rowdata_ + row*ncols_; for (unsigned col=0; coldata_ + row*newCols; const Numeric* off = data_ + row*ncols_; for (unsigned col=0; coldata_ + row*newCols; const unsigned rstart = row*n; const unsigned rstop = rstart + n; for (unsigned col=0; col inline void Matrix::coarseSum( const unsigned n, const unsigned m, Matrix* result) const { if (!nrows_) throw std::runtime_error( "In npstat::Matrix::coarseSum: uninitialized matrix"); if (!(n && m)) throw std::invalid_argument( "In npstat::Matrix::coarseSum: " "number of divisions must be positive"); if (nrows_ % n) throw std::invalid_argument( "In npstat::Matrix::coarseSum: " "n is not a divisor of the number or rows"); if (ncols_ % m) throw std::invalid_argument( "In npstat::Matrix::coarseSum: " "m is not a divisor of the number or columns"); makeCoarse(n, m, result); } template inline void Matrix::coarseAverage( const unsigned n, const unsigned m, Matrix* result) const { if (!nrows_) throw std::runtime_error( "In npstat::Matrix::coarseAverage: uninitialized matrix"); const unsigned nm = n*m; if (!nm) throw std::invalid_argument( "In npstat::Matrix::coarseAverage: " "number of divisions must be positive"); if (nrows_ % n) throw std::invalid_argument( "In npstat::Matrix::coarseAverage: " "n is not a divisor of the number or rows"); if (ncols_ % m) throw std::invalid_argument( "In npstat::Matrix::coarseAverage: " "m is not a divisor of the number or columns"); makeCoarse(n, m, result); if (nm > 1U) *result /= static_cast(nm); } template inline Matrix Matrix::removeRowAndColumn( const unsigned rRow, const unsigned rColumn) const { if (!nrows_) throw std::runtime_error( "In npstat::Matrix::removeRowAndColumn: uninitialized matrix"); unsigned newRows = nrows_; unsigned newCols = ncols_; if (rRow < nrows_) --newRows; if (rColumn < ncols_) --newCols; if (!(newRows && newCols)) throw std::invalid_argument( "In npstat::Matrix::removeRowAndColumn: can not remove the only " "row and/or column"); Matrix result(newRows, newCols); unsigned fromrow = 0U; for (unsigned torow=0U; torow inline Matrix Matrix::outer(const Matrix& r) const { if (!(nrows_ && r.nrows_)) throw std::runtime_error( "In npstat::Matrix::outer: uninitialized matrix"); const unsigned nrows = nrows_*r.nrows_; const unsigned ncols = ncols_*r.ncols_; Matrix result(nrows, ncols); for (unsigned tr=0; tr inline Matrix Matrix::operator*(const Matrix& r) const { if (!(nrows_ && r.nrows_)) throw std::runtime_error( "In npstat::Matrix::operator*: uninitialized matrix"); if (ncols_ != r.nrows_) throw std::domain_error( "In npstat::Matrix::operator*: incompatible matrix dimensions"); Matrix result(nrows_, r.ncols_); const bool thisDiag = isDiagonal(); const bool rDiag = r.isDiagonal(); if (thisDiag && rDiag) { result.zeroOut(); for (unsigned row=0; row inline void Matrix::times(const Matrix& r, Matrix* result) const { if (!(nrows_ && r.nrows_)) throw std::runtime_error( "In npstat::Matrix::times: uninitialized matrix"); if (ncols_ != r.nrows_) throw std::domain_error( "In npstat::Matrix::times: incompatible matrix dimensions"); assert(result); if (result == this || result == &r) throw std::invalid_argument( "In npstat::Matrix::times: " "can not have result set to one of the arguments"); result->resize(nrows_, r.ncols_); Numeric* const rdata = result->data_; const bool thisDiag = isDiagonal(); const bool rDiag = r.isDiagonal(); if (thisDiag && rDiag) { result->zeroOut(); for (unsigned row=0; rowdiagonalityKnown_ = true; result->isDiagonal_ = true; } else if (thisDiag) { // scale each row of the matrix r for (unsigned row=0; row template inline Matrix Matrix::timesVector( const Num2* data, const unsigned dataLen) const { if (!nrows_) throw std::runtime_error( "In npstat::Matrix::timesVector: uninitialized matrix"); if (ncols_ != dataLen) throw std::domain_error( "In npstat::Matrix::timesVector: incompatible vector length"); assert(data); Matrix result(nrows_, 1U); Numeric* res = result[0]; if (isDiagonal()) { for (unsigned row=0; row template inline Numeric Matrix::timesVector( const unsigned row, const Num2* data, const unsigned dataLen) const { if (!nrows_) throw std::runtime_error( "In npstat::Matrix::timesVector: uninitialized matrix"); if (ncols_ != dataLen) throw std::domain_error( "In npstat::Matrix::timesVector: incompatible vector length"); if (row >= nrows_) throw std::invalid_argument( "In npstat::Matrix::timesVector: row number out of range"); assert(data); if (isDiagonal()) return data_[row*ncols_ + row]*data[row]; else { const Numeric* prow = data_ + row*ncols_; Numeric sum = Numeric(); for (unsigned k=0; k template inline void Matrix::timesVector( const Num2* data, const unsigned dataLen, Num3* res, const unsigned resultLen) const { if (!nrows_) throw std::runtime_error( "In npstat::Matrix::timesVector: uninitialized matrix"); if (nrows_ > resultLen) throw std::invalid_argument( "In npstat::Matrix::timesVector: " "insufficient length of the output buffer"); if (ncols_ != dataLen) throw std::domain_error( "In npstat::Matrix::timesVector: incompatible vector length"); assert(data); assert(res); if (isDiagonal()) { for (unsigned row=0; row template inline Matrix Matrix::rowMultiply( const Num2* data, const unsigned dataLen) const { if (!nrows_) throw std::runtime_error( "In npstat::Matrix::rowMultiply: uninitialized matrix"); if (nrows_ != dataLen) throw std::domain_error( "In npstat::Matrix::rowMultiply: incompatible row length"); assert(data); Matrix result(1U, ncols_); Numeric* res = result[0]; if (isDiagonal()) { for (unsigned col=0; col template inline Numeric Matrix::rowMultiply( const unsigned col, const Num2* data, const unsigned dataLen) const { if (!nrows_) throw std::runtime_error( "In npstat::Matrix::rowMultiply: uninitialized matrix"); if (nrows_ != dataLen) throw std::domain_error( "In npstat::Matrix::rowMultiply: incompatible row length"); if (col >= ncols_) throw std::invalid_argument( "In npstat::Matrix::rowMultiply: column number out of range"); assert(data); if (isDiagonal()) return data[col]*data_[col*ncols_ + col]; else { const Numeric* pcol = data_ + col; Numeric sum = Numeric(); for (unsigned k=0; k template inline void Matrix::rowMultiply( const Num2* data, const unsigned dataLen, Num3* res, const unsigned resultLen) const { if (!nrows_) throw std::runtime_error( "In npstat::Matrix::rowMultiply: uninitialized matrix"); if (ncols_ > resultLen) throw std::invalid_argument( "In npstat::Matrix::rowMultiply: " "insufficient length of the output buffer"); if (nrows_ != dataLen) throw std::domain_error( "In npstat::Matrix::rowMultiply: incompatible row length"); assert(data); assert(res); if (isDiagonal()) { for (unsigned col=0; col template inline Numeric Matrix::bilinear( const Num2* data, const unsigned dataLen) const { typedef typename PreciseType::type Precise; if (!nrows_) throw std::runtime_error( "In npstat::Matrix::bilinear: uninitialized matrix"); if (nrows_ != ncols_) throw std::domain_error( "In npstat::Matrix::bilinear: matrix is not square"); if (ncols_ != dataLen) throw std::domain_error( "In npstat::Matrix::bilinear: incompatible vector length"); assert(data); Precise sum = Precise(); if (isDiagonal()) { for (unsigned row=0; row(sum); } template inline bool Matrix::solveLinearSystem( const Numeric* rhs, const unsigned lenRhs, Numeric* solution) const { if (!nrows_) throw std::runtime_error( "In npstat::Matrix::solveLinearSystem: uninitialized matrix"); if (nrows_ != ncols_) throw std::domain_error( "In npstat::Matrix::solveLinearSystem: matrix is not square"); if (nrows_ != lenRhs) throw std::invalid_argument( "In npstat::Matrix::solveLinearSystem: " "incompatible dimensionality of the right hand side vector"); assert(rhs); assert(solution); return solve_linear_system(data_, lenRhs, rhs, solution); } template inline bool Matrix::underdeterminedLinearSystem( const Numeric* rhs, const unsigned lenRhs, const Matrix& V, Numeric* solution, const unsigned lenSolution, Numeric* resultNormSq, Matrix* pa) const { if (!nrows_) throw std::runtime_error( "In npstat::Matrix::underdeterminedLinearSystem: " "uninitialized matrix"); if (nrows_ != lenRhs) throw std::invalid_argument( "In npstat::Matrix::underdeterminedLinearSystem: " "incompatible dimensionality of the right hand side vector"); if (ncols_ != lenSolution) throw std::invalid_argument( "In npstat::Matrix::underdeterminedLinearSystem: " "incompatible dimensionality of the solution vector"); if (nrows_ >= ncols_) throw std::domain_error( "In npstat::Matrix::underdeterminedLinearSystem: " "the system is not underdetermined"); if (!V.isSquare()) throw std::invalid_argument( "In npstat::Matrix::underdeterminedLinearSystem: " "covariance matrix is not square"); if (V.nRows() != ncols_) throw std::invalid_argument( "In npstat::Matrix::underdeterminedLinearSystem: " "incompatible covariance matrix dimensionality"); assert(rhs); assert(solution); const Matrix& vct = V.timesT(*this); const Matrix& cvctinv = (*this * vct).symPDInv(); const Matrix& A = vct*cvctinv; A.timesVector(rhs, lenRhs, solution, lenSolution); if (pa) { if (!A.isCompatible(*pa)) pa->uninitialize(); *pa = A; } if (resultNormSq) { const Matrix& B = this->Ttimes(cvctinv); *resultNormSq = Numeric(); for (unsigned i=0; i inline bool Matrix::linearLeastSquares( const Numeric* rhs, const unsigned lenRhs, Numeric* solution, const unsigned lenSolution) const { if (!nrows_) throw std::runtime_error( "In npstat::Matrix::linearLeastSquares: uninitialized matrix"); if (nrows_ != lenRhs) throw std::invalid_argument( "In npstat::Matrix::linearLeastSquares: " "incompatible dimensionality of the right hand side vector"); if (ncols_ != lenSolution) throw std::invalid_argument( "In npstat::Matrix::linearLeastSquares: " "incompatible dimensionality of the solution vector"); if (nrows_ <= ncols_) throw std::domain_error( "In npstat::Matrix::linearLeastSquares: " "the system is not overdetermined"); assert(rhs); assert(solution); return linear_least_squares(data_, nrows_, ncols_, rhs, solution); } template inline bool Matrix::weightedLeastSquares( const Numeric* rhsVec, const unsigned lenRhs, const Matrix& vinv, Numeric* solution, const unsigned lenSolution, Numeric* resultChiSquare, Matrix* resultCovarianceMatrix) const { if (!nrows_) throw std::runtime_error( "In npstat::Matrix::weightedLeastSquares: uninitialized matrix"); if (nrows_ != lenRhs) throw std::invalid_argument( "In npstat::Matrix::weightedLeastSquares: " "incompatible dimensionality of the right hand side vector"); if (ncols_ != lenSolution) throw std::invalid_argument( "In npstat::Matrix::weightedLeastSquares: " "incompatible dimensionality of the solution vector"); if (vinv.nrows_ != nrows_) throw std::invalid_argument( "In npstat::Matrix::weightedLeastSquares: " "incompatible inverse covariance matrix dimensionality"); if (!vinv.isSquare()) throw std::invalid_argument( "In npstat::Matrix::weightedLeastSquares: " "inverse covariance matrix is not square"); if (nrows_ <= ncols_) throw std::domain_error( "In npstat::Matrix::weightedLeastSquares: " "the system is not overdetermined"); assert(rhsVec); assert(solution); Matrix rhs(lenRhs, 1, rhsVec); const Matrix& rcov = vinv.bilinear(*this).symPDInv(); const Matrix& sol = rcov * Ttimes(vinv * rhs); assert(sol.nrows_ == lenSolution); copyBuffer(solution, sol.data_, lenSolution); if (resultChiSquare) { rhs -= *this * sol; *resultChiSquare = vinv.bilinear(rhs.data_, lenRhs); } if (resultCovarianceMatrix) { if (!rcov.isCompatible(*resultCovarianceMatrix)) resultCovarianceMatrix->uninitialize(); *resultCovarianceMatrix = rcov; } return true; } template inline bool Matrix::constrainedLeastSquares( const Numeric* rhs1, const unsigned lenRhs1, const Matrix& B, const Numeric* rhs2, const unsigned lenRhs2, Numeric* solution, const unsigned lenSolution, Numeric* resultChiSquare, Matrix* resultCovarianceMatrix, Numeric* unconstrainedSolution, Matrix* unconstrainedCovmat, Matrix* projectionMatrix, Matrix* pA) const { if (!nrows_) throw std::runtime_error( "In npstat::Matrix::constrainedLeastSquares: uninitialized matrix"); if (nrows_ != lenRhs1) throw std::invalid_argument( "In npstat::Matrix::constrainedLeastSquares: " "incompatible dimensionality of the right hand side vector"); if (ncols_ != lenSolution) throw std::invalid_argument( "In npstat::Matrix::constrainedLeastSquares: " "incompatible dimensionality of the solution vector"); if (!B.nrows_) throw std::runtime_error( "In npstat::Matrix::constrainedLeastSquares: " "uninitialized constraint matrix"); if (B.ncols_ != lenSolution) throw std::runtime_error( "In npstat::Matrix::constrainedLeastSquares: " "incompatible constraint matrix"); if (B.nrows_ != lenRhs2) throw std::invalid_argument( "In npstat::Matrix::constrainedLeastSquares: " "incompatible dimensionality of the constraint vector"); if (nrows_ + lenRhs2 <= ncols_) throw std::domain_error( "In npstat::Matrix::constrainedLeastSquares: " "the system is not overdetermined"); assert(rhs1); assert(rhs2); assert(solution); bool status = false; if (resultCovarianceMatrix || unconstrainedSolution || unconstrainedCovmat || projectionMatrix || pA) { if (nrows_ < ncols_) throw std::domain_error( "In npstat::Matrix::constrainedLeastSquares: " "unconstrained system is underdetermined"); // Solution of the unconstrained problem const Matrix& cov0 = TtimesThis().symPDInv(); const Matrix& sol0 = cov0 * T().timesVector(rhs1, lenRhs1); assert(sol0.nrows_ == lenSolution); if (unconstrainedSolution) copyBuffer(unconstrainedSolution, sol0.data_, lenSolution); if (unconstrainedCovmat) { if (!unconstrainedCovmat->isCompatible(cov0)) unconstrainedCovmat->uninitialize(); *unconstrainedCovmat = cov0; } // Solve the constrained problem const Matrix& tmp = cov0.timesT(B); const Matrix& A = tmp * (B * tmp).symPDInv(); const Matrix& CP = Matrix(lenSolution, lenSolution, 1) - A * B; const Matrix& sol = CP * sol0 + A.timesVector(rhs2, lenRhs2); copyBuffer(solution, sol.data_, lenSolution); if (projectionMatrix) { if (!projectionMatrix->isCompatible(CP)) projectionMatrix->uninitialize(); *projectionMatrix = CP; } if (pA) { if (!pA->isCompatible(A)) pA->uninitialize(); *pA = A; } if (resultCovarianceMatrix) { // Evaluate the result covariance matrix const Matrix& rcov = cov0.bilinearT(CP); if (!rcov.isCompatible(*resultCovarianceMatrix)) resultCovarianceMatrix->uninitialize(); *resultCovarianceMatrix = rcov; } status = true; } else status = constrained_least_squares(data_, nrows_, ncols_, rhs1, B.data_, B.nrows_, B.ncols_, rhs2, solution); if (resultChiSquare) { if (status) { long double chisq = 0.0L; for (unsigned i=0; i(chisq); } else *resultChiSquare = static_cast(-1); } return status; } template inline bool Matrix::solveLinearSystems( const Matrix& rhs, Matrix* solution) const { if (!nrows_) throw std::runtime_error( "In npstat::Matrix::solveLinearSystems: uninitialized matrix"); if (nrows_ != ncols_) throw std::domain_error( "In npstat::Matrix::solveLinearSystems: matrix is not square"); if (nrows_ != rhs.nrows_) throw std::invalid_argument( "In npstat::Matrix::solveLinearSystems: " "incompatible dimensionality of the right hand side matrix"); assert(solution); if (solution == this || solution == &rhs) throw std::invalid_argument( "In npstat::Matrix::solveLinearSystems: " "can not have result set to one of the arguments"); solution->resize(rhs.nrows_, rhs.ncols_); const bool ok = solve_linear_systems(data_, rhs.nrows_, rhs.ncols_, rhs.data_, solution->data_); if (!ok) { const Numeric zero = Numeric(); const unsigned len = solution->length(); Numeric* sdata = solution->data_; for (unsigned i=0; i inline Numeric Matrix::minValue() const { if (!nrows_) throw std::runtime_error( "In npstat::Matrix::minValue: uninitialized matrix"); Numeric minval = data_[0]; for (unsigned i=1; i::less(data_[i], minval)) minval = data_[i]; return minval; } template inline std::pair Matrix::argmin() const { if (!nrows_) throw std::runtime_error( "In npstat::Matrix::minValue: uninitialized matrix"); Numeric minval = data_[0]; std::pair location(0U, 0U); for (unsigned i=1; i::less(data_[i], minval)) { minval = data_[i]; location = std::pair(i/ncols_, i%ncols_); } return location; } template inline Numeric Matrix::maxValue() const { if (!nrows_) throw std::runtime_error( "In npstat::Matrix::maxValue: uninitialized matrix"); Numeric maxval = data_[0]; for (unsigned i=1; i::less(maxval, data_[i])) maxval = data_[i]; return maxval; } template inline std::pair Matrix::argmax() const { if (!nrows_) throw std::runtime_error( "In npstat::Matrix::maxValue: uninitialized matrix"); Numeric maxval = data_[0]; std::pair location(0U, 0U); for (unsigned i=1; i::less(maxval, data_[i])) { maxval = data_[i]; location = std::pair(i/ncols_, i%ncols_); } return location; } template inline Numeric Matrix::maxAbsValue() const { if (!nrows_) throw std::runtime_error( "In npstat::Matrix::maxAbsValue: uninitialized matrix"); Numeric maxval = std::abs(data_[0]); for (unsigned i=1; i maxval) maxval = std::abs(data_[i]); return maxval; } template inline std::pair Matrix::argmaxAbs() const { if (!nrows_) throw std::runtime_error( "In npstat::Matrix::maxAbsValue: uninitialized matrix"); Numeric maxval = std::abs(data_[0]); std::pair location(0U, 0U); for (unsigned i=1; i maxval) { maxval = std::abs(data_[i]); location = std::pair(i/ncols_, i%ncols_); } return location; } template inline Numeric Matrix::tr() const { typedef typename PreciseType::type Precise; if (!nrows_) throw std::runtime_error( "In npstat::Matrix::tr: uninitialized matrix"); if (nrows_ != ncols_) throw std::domain_error( "In npstat::Matrix::tr: matrix is not square"); Precise sum = Precise(); for (unsigned row=0; row(sum); } template inline Numeric Matrix::productTr(const Matrix& r) const { typedef typename PreciseType::type Precise; if (!(nrows_ && r.nrows_)) throw std::runtime_error( "In npstat::Matrix::productTr: uninitialized matrix"); if (ncols_ != r.nrows_ || nrows_ != r.ncols_) throw std::domain_error( "In npstat::Matrix::productTr: incompatible matrix dimensions"); Precise dsum = Precise(); for (unsigned row=0; row(dsum); } template inline Numeric Matrix::det() const { if (!nrows_) throw std::runtime_error( "In npstat::Matrix::det: uninitialized matrix"); if (nrows_ != ncols_) throw std::domain_error( "In npstat::Matrix::det: matrix is not square"); Numeric dvalue = Numeric(); if (nrows_ == 1) dvalue = data_[0]; else if (nrows_ == 2) dvalue = data_[0]*data_[3] - data_[1]*data_[2]; else if (nrows_ == 3) dvalue = -(data_[2]*data_[4]*data_[6]) + data_[1]*data_[5]*data_[6] + data_[2]*data_[3]*data_[7] - data_[0]*data_[5]*data_[7] - data_[1]*data_[3]*data_[8] + data_[0]*data_[4]*data_[8]; else if (nrows_ == 4) dvalue = data_[3]*data_[6]*data_[9]*data_[12] - data_[2]*data_[7]*data_[9]*data_[12] - data_[3]*data_[5]*data_[10]*data_[12] + data_[1]*data_[7]*data_[10]*data_[12] + data_[2]*data_[5]*data_[11]*data_[12] - data_[1]*data_[6]*data_[11]*data_[12] - data_[3]*data_[6]*data_[8]*data_[13] + data_[2]*data_[7]*data_[8]*data_[13] + data_[3]*data_[4]*data_[10]*data_[13] - data_[0]*data_[7]*data_[10]*data_[13] - data_[2]*data_[4]*data_[11]*data_[13] + data_[0]*data_[6]*data_[11]*data_[13] + data_[3]*data_[5]*data_[8]*data_[14] - data_[1]*data_[7]*data_[8]*data_[14] - data_[3]*data_[4]*data_[9]*data_[14] + data_[0]*data_[7]*data_[9]*data_[14] + data_[1]*data_[4]*data_[11]*data_[14] - data_[0]*data_[5]*data_[11]*data_[14] - data_[2]*data_[5]*data_[8]*data_[15] + data_[1]*data_[6]*data_[8]*data_[15] + data_[2]*data_[4]*data_[9]*data_[15] - data_[0]*data_[6]*data_[9]*data_[15] - data_[1]*data_[4]*data_[10]*data_[15] + data_[0]*data_[5]*data_[10]*data_[15]; else dvalue = lu_decomposition_matrix_det(data_, nrows_); return dvalue; } template inline Matrix Matrix::operator*(const Numeric r) const { if (!nrows_) throw std::runtime_error( "In npstat::Matrix::operator*: uninitialized matrix"); Matrix result(nrows_, ncols_); for (unsigned i=0; i inline void Matrix::times(const Numeric r, Matrix* result) const { if (!nrows_) throw std::runtime_error( "In npstat::Matrix::times: uninitialized matrix"); assert(result); if (result != this) result->resize(nrows_, ncols_); Numeric* const rdata = result->data_; for (unsigned i=0; idiagonalityKnown_ = true; result->isDiagonal_ = true; } } template inline Matrix Matrix::operator/(const Numeric r) const { if (!nrows_) throw std::runtime_error( "In npstat::Matrix::operator/: uninitialized matrix"); if (r == static_cast(0)) throw std::domain_error( "In npstat::Matrix::operator/: division by zero"); Matrix result(nrows_, ncols_); for (unsigned i=0; i inline void Matrix::over(const Numeric r, Matrix* result) const { if (!nrows_) throw std::runtime_error( "In npstat::Matrix::over: uninitialized matrix"); if (r == static_cast(0)) throw std::domain_error( "In npstat::Matrix::over: division by zero"); assert(result); if (result != this) result->resize(nrows_, ncols_); Numeric* const rdata = result->data_; for (unsigned i=0; idiagonalityKnown_ = true; result->isDiagonal_ = true; } } template inline Matrix Matrix::operator+() const { return *this; } template inline Matrix Matrix::operator-() const { if (!nrows_) throw std::runtime_error( "In npstat::Matrix::operator-(unary): uninitialized matrix"); Matrix result(nrows_, ncols_); for (unsigned i=0; i inline Matrix& Matrix::operator*=(const Numeric r) { if (!nrows_) throw std::runtime_error( "In npstat::Matrix::operator*=: uninitialized matrix"); for (unsigned i=0; i 0 if (!(diagonalityKnown_ && isDiagonal_)) diagonalityKnown_ = false; return *this; } template inline Matrix& Matrix::operator/=(const Numeric r) { if (!nrows_) throw std::runtime_error( "In npstat::Matrix::operator/=: uninitialized matrix"); if (r == static_cast(0)) throw std::domain_error( "In npstat::Matrix::operator/=: division by zero"); for (unsigned i=0; i 0 if (!(diagonalityKnown_ && isDiagonal_)) diagonalityKnown_ = false; return *this; } template inline Matrix& Matrix::operator+=(const Matrix& r) { if (!(nrows_ && r.nrows_)) throw std::runtime_error( "In npstat::Matrix::operator+=: uninitialized matrix"); if (!(nrows_ == r.nrows_ && ncols_ == r.ncols_)) throw std::domain_error( "In npstat::Matrix::operator+=: " "incompatible argument dimensions"); for (unsigned i=0; i inline Matrix& Matrix::operator-=(const Matrix& r) { if (!(nrows_ && r.nrows_)) throw std::runtime_error( "In npstat::Matrix::operator-=: uninitialized matrix"); if (!(nrows_ == r.nrows_ && ncols_ == r.ncols_)) throw std::domain_error( "In npstat::Matrix::operator-=: " "incompatible argument dimensions"); for (unsigned i=0; i inline const char* Matrix::classname() { static const std::string name( gs::template_class_name("npstat::Matrix")); return name.c_str(); } template inline bool Matrix::write(std::ostream& os) const { if (!nrows_) throw std::runtime_error( "In npstat::Matrix::write: uninitialized matrix"); assert(data_); gs::write_pod(os, nrows_); gs::write_pod(os, ncols_); gs::write_array(os, data_, len_); return !os.fail(); } template inline void Matrix::restore(const gs::ClassId& id, std::istream& in, Matrix* m) { static const gs::ClassId current( gs::ClassId::makeId >()); current.ensureSameId(id); assert(m); m->uninitialize(); gs::read_pod(in, &m->nrows_); gs::read_pod(in, &m->ncols_); m->len_ = m->nrows_ * m->ncols_; if (m->len_ <= Len) m->data_ = m->local_; else m->data_ = new Numeric[m->len_]; gs::read_array(in, m->data_, m->len_); if (in.fail()) { m->uninitialize(); throw gs::IOReadFailure("In npstat::Matrix::restore: " "input stream failure"); } } template inline unsigned Matrix::nonZeros() const { if (!nrows_) throw std::runtime_error( "In npstat::Matrix::nonZeros: uninitialized matrix"); const Numeric null = Numeric(); unsigned cnt = 0U; for (unsigned i=0; i inline Matrix diag(const Numeric* data, const unsigned nrows) { Matrix m(nrows, nrows); m.diagFill(data, nrows); return m; } template inline Matrix diag(const Numeric* data, const unsigned nrows, const unsigned ncols) { Matrix m(nrows, ncols, 0); assert(data); const unsigned imax = std::min(nrows, ncols); for (unsigned i=0; i inline std::ostream& operator<<( std::ostream& os, const npstat::Matrix& m) { // Figure out the largest string length needed to print an element const unsigned nrows = m.nRows(); const unsigned ncols = m.nColumns(); const unsigned len = nrows*ncols; if (!len) throw std::runtime_error( "In ::operator<<: can not write uninitialized npstat::Matrix object"); const N* data = m.data(); unsigned maxlen = 0; for (unsigned i=0; i maxlen) maxlen = siz; } ++maxlen; for (unsigned i=0; i #include #include #include #include #include "geners/CPP11_auto_ptr.hh" #include "npstat/nm/OrthoPolyMethod.hh" #include "npstat/nm/Matrix.hh" #include "npstat/nm/SimpleFunctors.hh" namespace npstat { class StorablePolySeries1D; 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 */ template 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; //@} /** // Average values of the polynomials for a sample of points. // Size of the "averages" array should be at least maxdeg + 1. */ template void sampleAverages(const Numeric* coords, unsigned long lenCoords, double* averages, unsigned maxdeg) const; /** // Average values of the pairwise polynomial products for // a sample of points. The returned matrix will be symmetric // and will have the dimensions (maxdeg + 1) x (maxdeg + 1). */ template Matrix sampleProductAverages(const Numeric* coords, unsigned long lenCoords, unsigned maxdeg) 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 given function in such a way that the integral of // the truncated series on the [xmin, xmax] interval is constrained // to "integralValue". 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). The construction minimizes the ISE weighted // by the empirical density. // // The function returns the chi-square from the corresponding // constrained least squares problem. This is the sum of // (coeffs[i] - c[i])^2, 0 <= i <= maxdeg, where c[i] are // determined by the "calculateCoeffs" method. */ template double constrainedCoeffs(const Functor& fcn, double *coeffs, unsigned maxdeg, double xmin, double xmax, double integralValue = 1.0) const; /** // Build the coefficients of the orthogonal polynomial series // for the discrete weight values (that is, fcn(x_i) = w_i, // using the terminology of the "calculateCoeffs" method). This // can sometimes be useful for weighted measures. Of course, // for unweighted measures this results in just delta_{deg,0}. */ void weightCoeffs(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; /** + // Unweighted triple product integrals arranged in a matrix. + // The returned matrix will be dimensioned pairs.size() x (maxdeg+1). + */ + Matrix tripleProductMatrix( + const std::vector >& pairs, + unsigned maxdeg, 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; /** Covariance between two cov4 estimates (i.e., error on the error) */ double covCov4(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. - */ + // 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; + /** + // Covariance matrix created by multiple calls to "epsCovariance" + */ + Matrix epsCovarianceMatrix( + const std::vector >& pairs, + bool highOrder) const; + /** Construct a StorablePolySeries1D object out of this one */ CPP11_auto_ptr makeStorablePolySeries( double xmin, double xmax, const double *coeffs=0, unsigned maxdeg=0) 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_; #ifdef SWIG public: inline StorablePolySeries1D* makeStorablePolySeries_2( double xmin, double xmax, const std::vector& coeffs) const { CPP11_auto_ptr ptr; if (coeffs.empty()) ptr = this->makeStorablePolySeries(xmin, xmax); else ptr = this->makeStorablePolySeries(xmin, xmax, &coeffs[0], coeffs.size() - 1); return ptr.release(); } #endif }; } #include "npstat/nm/ContOrthoPoly1D.icc" #endif // NPSTAT_CONTORTHOPOLY1D_HH_