Page MenuHomeHEPForge
Files F7877556

No OneTemporary
Actions

View Options

Index: trunk/npstat/nm/ContOrthoPoly1D.cc
===================================================================
--- trunk/npstat/nm/ContOrthoPoly1D.cc (revision 522)
+++ trunk/npstat/nm/ContOrthoPoly1D.cc (revision 523)
@@ -1,452 +1,472 @@
#include <algorithm>
#include <cfloat>
#include "npstat/nm/ContOrthoPoly1D.hh"
#include "npstat/nm/PairCompare.hh"
inline static double kdelta(const unsigned i,
const unsigned j)
{
return i == j ? 1.0 : 0.0;
}
namespace npstat {
ContOrthoPoly1D::ContOrthoPoly1D(const unsigned maxDegree,
const std::vector<double>& inCoords,
const OrthoPolyMethod m)
: wsum_(inCoords.size()),
wsumsq_(inCoords.size()),
minX_(DBL_MAX),
maxX_(-DBL_MAX),
maxdeg_(maxDegree),
allWeightsEqual_(true)
{
if (inCoords.empty()) throw std::invalid_argument(
"In npstat::ContOrthoPoly1D constructor: empty list of coordinates");
- const long double xsum = std::accumulate(inCoords.begin(), inCoords.end(), 0.0L);
+
const unsigned long sz = inCoords.size();
+ long double xsum = std::accumulate(inCoords.begin(), inCoords.end(), 0.0L);
meanX_ = xsum/sz;
- measure_.reserve(sz);
+
+ // Improve the numerical precision of the mean x subtraction
+ xsum = 0.0L;
for (unsigned long i = 0; i < sz; ++i)
{
const double x = inCoords[i];
if (x < minX_)
minX_ = x;
if (x > maxX_)
maxX_ = x;
- measure_.push_back(MeasurePoint(x - meanX_, 1.0));
+ xsum += (x - meanX_);
}
+ const double dx = xsum/sz;
+
+ // Fill the measure
+ measure_.reserve(sz);
+ for (unsigned long i = 0; i < sz; ++i)
+ measure_.push_back(MeasurePoint(inCoords[i] - meanX_ - dx, 1.0));
+
+ meanX_ += dx;
calcRecurrenceCoeffs(m);
}
ContOrthoPoly1D::ContOrthoPoly1D(const unsigned maxDegree,
const std::vector<MeasurePoint>& 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<MeasurePoint>());
if (measure_[0].second < 0.0) throw std::invalid_argument(
"In npstat::ContOrthoPoly1D constructor: negative measure entry found");
// 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;
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<degree; ++k)
{
const long double p = (x - rCoeffs_[k].alpha)*polyk -
rCoeffs_[k].beta*polykm1;
polykm1 = polyk;
polyk = p;
}
return polyk;
}
long double ContOrthoPoly1D::normpoly(const unsigned degree,
const long double x) const
{
long double polyk = 1.0L, polykm1 = 0.0L;
for (unsigned k=0; k<degree; ++k)
{
const long double 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
{
if (deg > 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<long double,long double> 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<degmax; ++k)
{
if (k == deg1)
p1 = polyk;
if (k == deg2)
p2 = polyk;
const long double 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<long double,long double>(p1, p2);
}
long double ContOrthoPoly1D::fournormpolyprod(const unsigned deg1,
const unsigned deg2,
const unsigned deg3,
const unsigned deg4,
const long double x) const
{
long double prod = 1.0L, polyk = 1.0L, polykm1 = 0.0L;
const unsigned degmax = std::max(std::max(deg1, deg2), std::max(deg3, deg4));
for (unsigned k=0; k<degmax; ++k)
{
if (k == deg1)
prod *= polyk;
if (k == deg2)
prod *= polyk;
if (k == deg3)
prod *= polyk;
if (k == deg4)
prod *= polyk;
const long double p = ((x - rCoeffs_[k].alpha)*polyk -
rCoeffs_[k].sqrbeta*polykm1)/rCoeffs_[k+1].sqrbeta;
polykm1 = polyk;
polyk = p;
}
if (degmax == deg1)
prod *= polyk;
if (degmax == deg2)
prod *= polyk;
if (degmax == deg3)
prod *= polyk;
if (degmax == deg4)
prod *= polyk;
return prod;
}
long double ContOrthoPoly1D::normseries(const double *coeffs, const unsigned maxdeg,
const long double x) const
{
long double sum = coeffs[0];
long double polyk = 1.0L, polykm1 = 0.0L;
for (unsigned k=0; k<maxdeg; ++k)
{
const long double p = ((x - rCoeffs_[k].alpha)*polyk -
rCoeffs_[k].sqrbeta*polykm1)/rCoeffs_[k+1].sqrbeta;
sum += p*coeffs[k+1];
polykm1 = polyk;
polyk = p;
}
return sum;
}
std::pair<long double,long double>
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<long double,long double>(xsum/wsum_, sum/wsum_);
}
else
// Average x should be 0 for degree == 0
return std::pair<long double,long double>(0.0L, 1.0L);
}
void ContOrthoPoly1D::calcRecurrenceStieltjes()
{
long double prevSprod = 1.0L;
for (unsigned deg=0; deg<=maxdeg_; ++deg)
{
const std::pair<long double,long double> ip = monicInnerProducts(deg);
rCoeffs_.push_back(Recurrence(ip.first/ip.second, ip.second/prevSprod));
prevSprod = ip.second;
}
}
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_);
}
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_);
}
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<long double, long double>& p =
twonormpoly(ideg1, ideg2, measure_[i].first);
sum += measure_[i].second*p.first*p.second;
}
return sum/wsum_;
}
double ContOrthoPoly1D::jointPowerAverage(
const unsigned ideg1, const unsigned power1,
const unsigned ideg2, const unsigned power2) const
{
if (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<long double, long double>& 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 deg1, const unsigned deg2,
const unsigned deg3, const unsigned deg4) const
{
if (deg1 > maxdeg_ || deg2 > maxdeg_ ||
deg3 > maxdeg_ || deg4 > 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*fournormpolyprod(
deg1, deg2, deg3, deg4, measure_[i].first);
return sum/wsum_;
}
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)))
{
cov = (jointAverage(deg1, deg2, deg3, deg4) -
kdelta(deg1, deg2)*kdelta(deg3, deg4))/effectiveSampleSize();
if (deg1 == deg3 && deg2 == 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<double> 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<double> 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;
}
void ContOrthoPoly1D::calculateRoots(
double *roots, const unsigned deg) const
{
if (deg)
{
assert(roots);
const Matrix<double>& mat = jacobiMatrix(deg);
if (deg == 1U)
roots[0] = mat[0][0];
else
mat.tdSymEigen(roots, deg);
for (unsigned i=0; i<deg; ++i)
roots[i] += meanX_;
}
}
double ContOrthoPoly1D::integratePoly(
const unsigned deg1, const unsigned power,
const double xmin, const double xmax) const
{
if (deg1 > 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_);
}
}

File Metadata

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

Event Timeline