Page MenuHomeHEPForge
Files F7877621

No OneTemporary
Actions

View Options

Index: trunk/npstat/nm/ContOrthoPoly1D.cc
===================================================================
--- trunk/npstat/nm/ContOrthoPoly1D.cc (revision 909)
+++ trunk/npstat/nm/ContOrthoPoly1D.cc (revision 910)
@@ -1,978 +1,979 @@
#include "npstat/nm/ContOrthoPoly1D.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 {
const ContOrthoPoly1D::precise_type ContOrthoPoly1D::precise_zero = 0.0;
void ContOrthoPoly1D::calcMeanXandWeightSums()
{
// Calculate wsum_, wsumsq_, and meanX_
const unsigned long mSize = measure_.size();
precise_type xsum = precise_zero;
if (allWeightsEqual_)
{
const precise_type lw0 = measure_[0].second;
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;
}
meanX_ = xsum/mSize;
}
else
{
for (unsigned long i = 0; i < mSize; ++i)
{
const double x = measure_[i].first;
const precise_type w = measure_[i].second;
wsum_ += w;
wsumsq_ += w*w;
xsum += w*x;
if (x < minX_)
minX_ = x;
if (x > maxX_)
maxX_ = x;
}
assert(wsum_ > precise_zero);
meanX_ = xsum/wsum_;
}
// Shift the measure
double shift = meanX_;
for (unsigned long i = 0; i < mSize; ++i)
measure_[i].first -= shift;
// Try to improve the mean
xsum = precise_zero;
precise_type dx;
if (allWeightsEqual_)
{
for (unsigned long i = 0; i < mSize; ++i)
xsum += measure_[i].first;
dx = xsum/mSize;
}
else
{
for (unsigned long i = 0; i < mSize; ++i)
{
const precise_type w = measure_[i].second;
xsum += w*measure_[i].first;
}
dx = xsum/wsum_;
}
// Shift the measure and the mean
meanX_ += dx;
shift = dx;
for (unsigned long i = 0; i < mSize; ++i)
measure_[i].first -= shift;
}
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_;
}
ContOrthoPoly1D::precise_type ContOrthoPoly1D::monicpoly(
const unsigned degree, const precise_type x) const
{
assert(degree <= maxdeg_);
precise_type polyk = 1.0, polykm1 = precise_zero;
for (unsigned k=0; k<degree; ++k)
{
const precise_type p = (x - rCoeffs_[k].alpha)*polyk -
rCoeffs_[k].beta*polykm1;
polykm1 = polyk;
polyk = p;
}
return polyk;
}
ContOrthoPoly1D::precise_type ContOrthoPoly1D::normpoly(
const unsigned degree, const precise_type x) const
{
assert(degree <= maxdeg_);
precise_type polyk = 1.0, polykm1 = precise_zero;
for (unsigned k=0; k<degree; ++k)
{
const precise_type p = ((x - rCoeffs_[k].alpha)*polyk -
rCoeffs_[k].sqrbeta*polykm1)/rCoeffs_[k+1].sqrbeta;
polykm1 = polyk;
polyk = p;
}
return polyk;
}
double ContOrthoPoly1D::powerAverage(const unsigned deg,
const unsigned power) const
{
switch (power)
{
case 0:
return 1.0;
case 1:
return kdelta(deg, 0U);
case 2:
return 1.0;
default:
{
if (deg > maxdeg_) throw std::invalid_argument(
"In npstat::ContOrthoPoly1D::powerAverage: "
"degree argument is out of range");
const unsigned long mSize = measure_.size();
precise_type sum = precise_zero;
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::precise_type,ContOrthoPoly1D::precise_type>
ContOrthoPoly1D::twonormpoly(
const unsigned deg1, const unsigned deg2, const precise_type x) const
{
precise_type p1 = precise_zero, p2 = precise_zero, polyk = 1.0, polykm1 = precise_zero;
const unsigned degmax = std::max(deg1, deg2);
for (unsigned k=0; k<degmax; ++k)
{
if (k == deg1)
p1 = polyk;
if (k == deg2)
p2 = polyk;
const precise_type p = ((x - rCoeffs_[k].alpha)*polyk -
rCoeffs_[k].sqrbeta*polykm1)/rCoeffs_[k+1].sqrbeta;
polykm1 = polyk;
polyk = p;
}
if (deg1 == degmax)
p1 = polyk;
if (deg2 == degmax)
p2 = polyk;
return std::pair<precise_type,precise_type>(p1, p2);
}
ContOrthoPoly1D::precise_type ContOrthoPoly1D::normpolyprod(
const unsigned* degrees, const unsigned nDeg,
const unsigned degmax, const precise_type x) const
{
precise_type prod = 1.0, polyk = 1.0, polykm1 = precise_zero;
for (unsigned k=0; k<degmax; ++k)
{
for (unsigned i=0; i<nDeg; ++i)
if (k == degrees[i])
prod *= polyk;
const precise_type p = ((x - rCoeffs_[k].alpha)*polyk -
rCoeffs_[k].sqrbeta*polykm1)/rCoeffs_[k+1].sqrbeta;
polykm1 = polyk;
polyk = p;
}
for (unsigned i=0; i<nDeg; ++i)
if (degmax == degrees[i])
prod *= polyk;
return prod;
}
std::pair<ContOrthoPoly1D::precise_type,ContOrthoPoly1D::precise_type>
ContOrthoPoly1D::monicInnerProducts(const unsigned degree) const
{
precise_type sum = precise_zero, xsum = precise_zero;
const unsigned long mSize = measure_.size();
for (unsigned long i = 0; i < mSize; ++i)
{
const precise_type x = measure_[i].first;
const precise_type p = monicpoly(degree, x);
const precise_type pprod = p*p*measure_[i].second;
sum += pprod;
xsum += x*pprod;
}
return std::pair<precise_type,precise_type>(xsum/wsum_, sum/wsum_);
}
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<precise_type> scalarProducts(maxdeg+1U, precise_zero);
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 precise_type f = w;
precise_type polyk = 1.0, polykm1 = precise_zero;
for (unsigned k=0; k<maxdeg; ++k)
{
scalarProducts[k] += polyk*f*w;
const precise_type poly = ((x - rCoeffs_[k].alpha)*polyk -
rCoeffs_[k].sqrbeta*polykm1)/rCoeffs_[k+1].sqrbeta;
polykm1 = polyk;
polyk = poly;
}
scalarProducts[maxdeg] += polyk*f*w;
}
for (unsigned deg=0; deg<=maxdeg; ++deg)
coeffs[deg] = scalarProducts[deg]/wsum_;
}
void ContOrthoPoly1D::calcRecurrenceStieltjes()
{
precise_type prevSprod = 1.0;
for (unsigned deg=0; deg<=maxdeg_; ++deg)
{
const std::pair<precise_type,precise_type> ip = monicInnerProducts(deg);
rCoeffs_.push_back(Recurrence(ip.first/ip.second, ip.second/prevSprod));
prevSprod = ip.second;
}
}
void ContOrthoPoly1D::calcRecurrenceLanczos()
{
const unsigned long mSize = measure_.size();
std::vector<precise_type> dmem(mSize*2UL);
precise_type* p0 = &dmem[0];
for (unsigned long i=0; i<mSize; ++i)
p0[i] = measure_[i].first;
precise_type* p1 = &dmem[mSize];
clearBuffer(p1, mSize);
p1[0] = measure_[0].second/wsum_;
const unsigned long mm1 = mSize - 1;
for (unsigned long n=0; n<mm1; ++n)
{
const unsigned long np1 = n + 1UL;
precise_type xlam = measure_[np1].first;
precise_type pn = measure_[np1].second/wsum_;
precise_type gam = 1.0, sig = 0.0, t = 0.0;
for (unsigned long k=0; k<=np1; ++k)
{
const precise_type rho = p1[k] + pn;
const precise_type tmp = gam*rho;
precise_type tsig = sig;
if (rho <= 0.0)
{
gam = 1.0;
sig = 0.0;
}
else
{
gam = p1[k]/rho;
sig = pn/rho;
}
const precise_type tk = sig*(p0[k] - xlam) - gam*t;
p0[k] -= tk - t;
t = tk;
if (sig < 0.0)
pn = tsig*p1[k];
else
pn = t*t/sig;
tsig = sig;
p1[k] = tmp;
}
}
for (unsigned deg=0; deg<=maxdeg_; ++deg)
rCoeffs_.push_back(Recurrence(p0[deg], p1[deg]));
}
double ContOrthoPoly1D::poly(const unsigned deg, const double x) const
{
if (deg > maxdeg_) throw std::invalid_argument(
"In npstat::ContOrthoPoly1D::poly: degree argument is out of range");
return normpoly(deg, x - meanX_);
}
std::pair<double,double> 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<precise_type,precise_type>& ld =
twonormpoly(deg1, deg2, x - meanX_);
return std::pair<double,double>(ld.first, ld.second);
}
ContOrthoPoly1D::precise_type 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");
precise_type sum = precise_zero;
const unsigned long mSize = measure_.size();
for (unsigned long i = 0; i < mSize; ++i)
{
const std::pair<precise_type, precise_type>& p =
twonormpoly(ideg1, ideg2, measure_[i].first);
sum += measure_[i].second*p.first*p.second;
}
return sum/wsum_;
}
Matrix<ContOrthoPoly1D::precise_type>
ContOrthoPoly1D::empiricalKroneckerMatrix(const unsigned maxdeg) const
{
if (maxdeg > maxdeg_) throw std::invalid_argument(
"In npstat::ContOrthoPoly1D::empiricalKroneckerMatrix: "
"matrix size is out of range");
Matrix<ContOrthoPoly1D::precise_type> kd(maxdeg+1U, maxdeg+1U, 0);
std::vector<ContOrthoPoly1D::precise_type> polyv(maxdeg+1U);
ContOrthoPoly1D::precise_type* polyBuf = &polyv[0];
const unsigned long mSize = measure_.size();
for (unsigned long ipt = 0; ipt < mSize; ++ipt)
{
const double weight = measure_[ipt].second;
getAllPolys(measure_[ipt].first, maxdeg, polyBuf);
for (unsigned i=0; i<=maxdeg; ++i)
{
const ContOrthoPoly1D::precise_type polyI = polyBuf[i];
ContOrthoPoly1D::precise_type* row = kd[i];
for (unsigned j=i; j<=maxdeg; ++j)
row[j] += weight*polyI*polyBuf[j];
}
}
for (unsigned i=0; i<maxdeg; ++i)
for (unsigned j=i+1U; j<=maxdeg; ++j)
kd[j][i] = kd[i][j];
- return kd/wsum_;
+ kd /= wsum_;
+ return kd;
}
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");
precise_type sum = precise_zero;
const unsigned long mSize = measure_.size();
for (unsigned long i = 0; i < mSize; ++i)
{
const std::pair<precise_type, precise_type>& 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<nDegrees && allNonZero; ++ideg)
allNonZero = degrees[ideg];
if (!allNonZero)
{
unsigned degBuf[256];
if (nDegrees - 1U <= sizeof(degBuf)/sizeof(degBuf[0]))
{
unsigned nonZeroDegs = 0;
for (unsigned ideg=0; ideg<nDegrees; ++ideg)
if (degrees[ideg])
degBuf[nonZeroDegs++] = degrees[ideg];
return jointAverage(degBuf, nonZeroDegs, sorted);
}
}
}
unsigned degmax;
if (sorted)
degmax = degrees[nDegrees-1U];
else
degmax = *std::max_element(degrees, degrees+nDegrees);
if (degmax > maxdeg_) throw std::invalid_argument(
"In npstat::ContOrthoPoly1D::jointAverage: "
"degree argument is out of range");
precise_type sum = precise_zero;
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<double,double>
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<double,double>(rc.alpha, rc.beta);
}
Matrix<ContOrthoPoly1D::precise_type>
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<precise_type> mat(n, n, 0);
unsigned ip1 = 1;
for (unsigned i=0; i<n; ++i, ++ip1)
{
mat[i][i] = rCoeffs_[i].alpha;
if (i)
mat[i][i-1] = rCoeffs_[i].sqrbeta;
if (ip1 < n)
mat[i][ip1] = rCoeffs_[ip1].sqrbeta;
}
return mat;
}
ContOrthoPoly1D::precise_type ContOrthoPoly1D::integratePoly(
const unsigned deg1, const unsigned power,
const double xmin, const double xmax) const
{
if (!deg1)
return xmax - xmin;
if (deg1 > maxdeg_) throw std::invalid_argument(
"In npstat::ContOrthoPoly1D::integratePoly: "
"degree argument is out of range");
const unsigned nGood = PreciseQuadrature::minimalExactRule(deg1*power);
if (!nGood) throw std::invalid_argument(
"In npstat::ContOrthoPoly1D::integratePoly: "
"product of poly degree and power is too large");
PreciseQuadrature glq(nGood);
PolyFcn fcn(*this, deg1, power);
return glq.integrate(fcn, xmin - meanX_, xmax - meanX_);
}
ContOrthoPoly1D::precise_type 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 = PreciseQuadrature::minimalExactRule(sumdeg);
if (!nGood) throw std::invalid_argument(
"In npstat::ContOrthoPoly1D::integrateTripleProduct: "
"sum of the polynomial degrees is too large");
PreciseQuadrature glq(nGood);
TripleProdFcn fcn(*this, deg1, deg2, deg3);
return glq.integrate(fcn, xmin - meanX_, xmax - meanX_);
}
Matrix<ContOrthoPoly1D::precise_type> ContOrthoPoly1D::tripleProductMatrix(
const std::vector<std::pair<unsigned,unsigned> >& pairs,
const unsigned maxdeg, const double xmin, const double xmax) const
{
typedef npstat::Triple<unsigned,unsigned,unsigned> Key3;
typedef std::map<Key3,precise_type> 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<precise_type> design0(nPairs, nPolys);
TripleProductMap tmap;
// Avoid calculating the same integral more than once
for (unsigned ipair=0; ipair<nPairs; ++ipair)
{
const unsigned n = pairs[ipair].first;
const unsigned m = pairs[ipair].second;
for (unsigned ipoly=0; ipoly<nPolys; ++ipoly)
{
unsigned degs[3];
degs[0] = m;
degs[1] = n;
degs[2] = ipoly;
std::sort(degs, degs+3);
const Key3 key(degs[0], degs[1], degs[2]);
TripleProductMap::iterator it = tmap.find(key);
if (it == tmap.end())
{
const precise_type value = integrateTripleProduct(
n, m, ipoly, xmin, xmax);
std::pair<TripleProductMap::iterator,bool> 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<MemoKey,double>::const_iterator it = cachedAverages_.find(key);
if (it == cachedAverages_.end())
{
value = jointAverage(key.degrees(), key.nDegrees(), true);
cachedAverages_.insert(std::pair<MemoKey,double>(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)
{
precise_type sum = precise_zero;
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);
precise_type sum = precise_zero;
// 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<double> ContOrthoPoly1D::epsCovarianceMatrix(
const std::vector<std::pair<unsigned,unsigned> >& 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<double> cov(nPairs, nPairs);
for (unsigned ipair1=0; ipair1<nPairs; ++ipair1)
{
const unsigned n = pairs[ipair1].first;
const unsigned m = pairs[ipair1].second;
for (unsigned ipair2=0; ipair2<=ipair1; ++ipair2)
{
const double value = epsCovariance(
n, m, pairs[ipair2].first, pairs[ipair2].second, highOrder);
cov[ipair1][ipair2] = value;
if (ipair1 != ipair2)
cov[ipair2][ipair1] = value;
}
}
return cov;
}
CPP11_auto_ptr<StorablePolySeries1D> ContOrthoPoly1D::makeStorablePolySeries(
const double i_xmin, const double i_xmax,
const double *coeffs, const unsigned maxdeg) const
{
const unsigned sz = rCoeffs_.size();
std::vector<std::pair<precise_type,precise_type> > rc(sz);
for (unsigned i=0; i<sz; ++i)
{
const Recurrence& r(rCoeffs_[i]);
rc[i].first = r.alpha;
rc[i].second = r.sqrbeta;
}
return CPP11_auto_ptr<StorablePolySeries1D>(
new StorablePolySeries1D(rc, i_xmin, i_xmax, meanX_, coeffs, maxdeg));
}
}

File Metadata

Mime Type
text/x-diff
Expires
Tue, Nov 19, 3:55 PM (1 d, 19 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
3805051
Default Alt Text
(37 KB)

Event Timeline