Page MenuHomeHEPForge
Files F10881763

No OneTemporary
Actions

View Options

Index: trunk/npstat/nm/ContOrthoPoly1D.hh
===================================================================
--- trunk/npstat/nm/ContOrthoPoly1D.hh (revision 514)
+++ trunk/npstat/nm/ContOrthoPoly1D.hh (revision 515)
@@ -1,233 +1,233 @@
#ifndef NPSTAT_CONTORTHOPOLY1D_HH_
#define NPSTAT_CONTORTHOPOLY1D_HH_
/*!
// \file ContOrthoPoly1D.hh
//
// \brief Continuous orthogonal polynomial series in one dimension
// generated for a discrete measure
//
// Author: I. Volobouev
//
// May 2017
*/
#include <vector>
#include <utility>
#include <cassert>
#include <cmath>
#include "npstat/nm/OrthoPolyMethod.hh"
#include "npstat/nm/Matrix.hh"
#include "npstat/nm/SimpleFunctors.hh"
namespace npstat {
class ContOrthoPoly1D
{
public:
// The first element of the pair is the coordinate and the
// second measure is the weight (all weights must be non-negative)
typedef std::pair<double,double> MeasurePoint;
ContOrthoPoly1D();
/**
// 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<MeasurePoint>& measure,
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 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 <class Functor>
void calculateCoeffs(const Functor& fcn,
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<double,double> recurrenceCoeffs(unsigned deg) const;
/**
// Generate principal minor of order n of the Jacobi matrix.
// n must not exceed "maxDegree"
*/
Matrix<double> 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 without the weight
// (so this is not an inner product with the poly of degree 0)
*/
- long double integratePoly(unsigned degree,
- long double xmin, long double xmax) const;
+ double integratePoly(unsigned degree,
+ double xmin, double xmax) const;
/**
// Integral of the product of three orthonormal polynomials
// without the weight (so this is not an inner product)
*/
- long double integrateTripleProduct(unsigned deg1, unsigned deg2,
- unsigned deg3, long double xmin,
- long double xmax) const;
+ double integrateTripleProduct(unsigned deg1, unsigned deg2,
+ unsigned deg3, double xmin,
+ double xmax) const;
/**
// A measure-weighted average of orthonormal poly values
// to some power
*/
double powerAverage(unsigned deg, unsigned power) const;
/**
// A measure-weighted average of two orthonormal poly values
// to some powers
*/
double jointPowerAverage(unsigned deg1, unsigned power1,
unsigned deg2, unsigned power2) const;
/** A measure-weighted average of four orthonormal poly values */
double jointAverage(unsigned deg1, unsigned deg2,
unsigned deg3, unsigned deg4) 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<long double, long double>
{
public:
inline PolyFcn(const ContOrthoPoly1D& p, const unsigned d1)
: poly(p), deg1(d1) {}
inline long double operator()(const long double& x) const
{return poly.normpoly(deg1, x);}
private:
const ContOrthoPoly1D& poly;
unsigned deg1;
};
class TripleProdFcn;
friend class TripleProdFcn;
class TripleProdFcn : public Functor1<long double, long double>
{
public:
inline TripleProdFcn(const ContOrthoPoly1D& p, const unsigned d1,
const unsigned d2, const unsigned d3)
: poly(p), deg1(d1), deg2(d2), deg3(d3) {}
inline long double operator()(const long double& x) const
{return poly.fournormpolyprod(0U, deg1, deg2, deg3, x);}
private:
const ContOrthoPoly1D& poly;
unsigned deg1;
unsigned deg2;
unsigned deg3;
};
void calcRecurrenceStieltjes();
long double monicpoly(unsigned degree, long double x) const;
std::pair<long double,long double> monicInnerProducts(
unsigned degree) const;
long double normpoly(unsigned degree, long double x) const;
std::pair<long double,long double> twonormpoly(
unsigned deg1, unsigned deg2, long double x) const;
long double fournormpolyprod(unsigned deg1, unsigned deg2,
unsigned deg3, unsigned deg4,
long double x) const;
long double normseries(const double *coeffs, unsigned maxdeg,
long double x) const;
std::vector<MeasurePoint> measure_;
std::vector<Recurrence> rCoeffs_;
long double wsum_;
long double wsumsq_;
double meanX_;
double minX_;
double maxX_;
unsigned maxdeg_;
bool allWeightsEqual_;
};
}
#include "npstat/nm/ContOrthoPoly1D.icc"
#endif // NPSTAT_CONTORTHOPOLY1D_HH_
Index: trunk/npstat/nm/ContOrthoPoly1D.cc
===================================================================
--- trunk/npstat/nm/ContOrthoPoly1D.cc (revision 514)
+++ trunk/npstat/nm/ContOrthoPoly1D.cc (revision 515)
@@ -1,424 +1,424 @@
#include <algorithm>
#include <cfloat>
#include "npstat/nm/ContOrthoPoly1D.hh"
#include "npstat/nm/PairCompare.hh"
#include "npstat/nm/GaussLegendreQuadrature.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<MeasurePoint>& inMeasure,
const OrthoPolyMethod m)
: measure_(inMeasure),
wsum_(0.0L),
wsumsq_(0.0L),
minX_(DBL_MAX),
maxX_(-DBL_MAX),
maxdeg_(maxDegree)
{
// 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;
allWeightsEqual_ = true;
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_;
// Calculate the recurrence coefficients
rCoeffs_.reserve(maxdeg_ + 1);
switch (m)
{
case OPOLY_STIELTJES:
calcRecurrenceStieltjes();
break;
default:
throw std::runtime_error(
"In npstat::ContOrthoPoly1D constructor: "
"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
{
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_);
}
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_;
}
}
- long double ContOrthoPoly1D::integratePoly(
- const unsigned deg1, const long double xmin,
- const long double xmax) const
+ double ContOrthoPoly1D::integratePoly(
+ const unsigned deg1, 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 std::vector<unsigned>& nAllowed = GaussLegendreQuadrature::allowedNPonts();
const unsigned sz = nAllowed.size();
unsigned nGood = 0;
for (unsigned i=0; i<sz && !nGood; ++i)
if (nAllowed[i]*2 > deg1)
nGood = nAllowed[i];
if (!nGood) throw std::invalid_argument(
"In npstat::ContOrthoPoly1D::integratePoly: "
"poly degree is too large");
GaussLegendreQuadrature glq(nGood);
PolyFcn fcn(*this, deg1);
return glq.integrate(fcn, xmin - meanX_, xmax - meanX_);
}
- long double ContOrthoPoly1D::integrateTripleProduct(
+ double ContOrthoPoly1D::integrateTripleProduct(
const unsigned deg1, const unsigned deg2,
- const unsigned deg3, const long double xmin,
- const long double xmax) const
+ 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 std::vector<unsigned>& nAllowed = GaussLegendreQuadrature::allowedNPonts();
const unsigned sz = nAllowed.size();
unsigned nGood = 0;
for (unsigned i=0; i<sz && !nGood; ++i)
if (nAllowed[i]*2 > sumdeg)
nGood = nAllowed[i];
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
Sat, May 3, 6:49 AM (22 h, 59 m)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4983144
Default Alt Text
(24 KB)

Event Timeline