Page MenuHomeHEPForge
Files F9501593

No OneTemporary
Actions

View Options

Index: trunk/npstat/nm/ContOrthoPoly1D.hh
===================================================================
--- trunk/npstat/nm/ContOrthoPoly1D.hh (revision 720)
+++ trunk/npstat/nm/ContOrthoPoly1D.hh (revision 721)
@@ -1,632 +1,640 @@
#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 <map>
#include "geners/CPP11_auto_ptr.hh"
#include "npstat/nm/OrthoPolyMethod.hh"
#include "npstat/nm/Matrix.hh"
#include "npstat/nm/SimpleFunctors.hh"
#include "npstat/nm/PreciseType.hh"
#ifdef USE_LAPACK_QUAD
#include "npstat/nm/GaussLegendreQuadratureQ.hh"
#else
#include "npstat/nm/GaussLegendreQuadrature.hh"
#endif
namespace npstat {
class StorablePolySeries1D;
class ContOrthoPoly1D
{
public:
typedef PreciseType<double>::type precise_type;
// The first element of the pair is the coordinate and the
// second element is the weight (all weights must be non-negative)
typedef std::pair<double,double> MeasurePoint;
#ifdef USE_LAPACK_QUAD
typedef GaussLegendreQuadratureQ PreciseQuadrature;
#else
typedef GaussLegendreQuadrature PreciseQuadrature;
#endif
/**
// Main constructor, with obvious arguments. The first element
// of the measure pair is the coordinate and the second element
// is the weight (all weights must be non-negative). 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.
*/
template<typename Numeric1, typename Numeric2>
ContOrthoPoly1D(unsigned maxDegree,
const std::vector<std::pair<Numeric1,Numeric2> >& measure,
OrthoPolyMethod m = OPOLY_STIELTJES);
/** Constructor which assumes that all weights are 1.0 */
template<typename Numeric>
ContOrthoPoly1D(unsigned maxDegree,
const std::vector<Numeric>& 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 precise_type sumOfWeights() const {return wsum_;}
inline precise_type sumOfWeightSquares() const {return wsumsq_;}
inline bool areAllWeightsEqual() const {return allWeightsEqual_;}
inline double minCoordinate() const {return minX_;}
inline double maxCoordinate() const {return maxX_;}
inline precise_type 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<double,double> 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.
*/
template <typename Numeric>
inline void allPolys(const double x, const unsigned maxdeg,
Numeric* polyValues) const
{getAllPolys(x - meanX_, maxdeg, polyValues);}
/**
+ // Function added for interface compatibility with
+ // AbsClassicalOrthoPoly1D
+ */
+ inline void allpoly(const long double x, long double* values,
+ const unsigned maxdeg) const
+ {getAllPolys(x - meanX_, maxdeg, values);}
+
+ /**
// Average values of the polynomials for a sample of points.
// Size of the "averages" array should be at least maxdeg + 1.
// Note that the resulting averages are not weighted by the measure.
*/
template <typename Numeric>
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).
// Note that the resulting averages are not weighted by the measure.
*/
template <typename Numeric>
Matrix<double> sampleProductAverages(const Numeric* coords,
unsigned long lenCoords,
unsigned maxdeg) const;
/** Return the value of the orthonormal polynomial series */
template <typename Numeric>
double series(const Numeric* coeffs, unsigned deg, 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, typename Numeric>
void calculateCoeffs(const Functor& fcn,
Numeric* 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 <class Functor>
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 <class Functor>
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 <class Functor>
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.
*/
precise_type 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.
// Higher precision than double is used internally.
*/
std::pair<double,double> recurrenceCoeffs(unsigned deg) const;
/**
// Generate principal minor of order n of the Jacobi matrix.
// n must not exceed "maxDegree"
*/
Matrix<precise_type> jacobiMatrix(unsigned n) const;
/**
// Roots of the polynomial with the given degree. Naturally,
// the degree argument must not exceed "maxDegree".
*/
void calculateRoots(precise_type* 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).
*/
precise_type 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)
*/
precise_type 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<precise_type> tripleProductMatrix(
const std::vector<std::pair<unsigned,unsigned> >& 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 <class Functor>
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.
*/
double epsCovariance(unsigned i, unsigned j,
unsigned m, unsigned n, bool highOrder) const;
/**
// Covariance matrix created by multiple calls to "epsCovariance"
*/
Matrix<double> epsCovarianceMatrix(
const std::vector<std::pair<unsigned,unsigned> >& pairs,
bool highOrder) const;
/** Construct a StorablePolySeries1D object out of this one */
CPP11_auto_ptr<StorablePolySeries1D> makeStorablePolySeries(
double xmin, double xmax,
const double *coeffs=0, unsigned maxdeg=0) const;
private:
static const precise_type precise_zero;
struct Recurrence
{
inline Recurrence(const precise_type a,
const precise_type b)
: alpha(a), beta(b)
{
assert(beta > 0.0);
sqrbeta = std::sqrt(beta);
}
precise_type alpha;
precise_type beta;
precise_type sqrbeta;
};
class PolyFcn;
friend class PolyFcn;
class PolyFcn : public Functor1<precise_type, precise_type>
{
public:
inline PolyFcn(const ContOrthoPoly1D& p,
const unsigned d1, const unsigned power)
: poly(p), deg1(d1), polypow(power) {}
inline precise_type operator()(const precise_type& x) const
{
precise_type p = 1.0;
if (polypow)
{
p = poly.normpoly(deg1, x);
switch (polypow)
{
case 1U:
break;
case 2U:
p *= p;
break;
default:
const precise_type myPow = polypow;
p = std::pow(p, myPow);
break;
}
}
return p;
}
private:
const ContOrthoPoly1D& poly;
unsigned deg1;
unsigned polypow;
};
class TripleProdFcn;
friend class TripleProdFcn;
class TripleProdFcn : public Functor1<precise_type, precise_type>
{
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 precise_type operator()(const precise_type& x) const
{return poly.normpolyprod(degs, 3, degmax, x);}
private:
const ContOrthoPoly1D& poly;
unsigned degs[3];
unsigned degmax;
};
template<class Funct> class EWPolyProductFcn;
template<class Funct> friend class EWPolyProductFcn;
template<class Funct>
class EWPolyProductFcn : public Functor1<precise_type, precise_type>
{
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 precise_type operator()(const precise_type& 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<nDegs_; ++i)
{
if (degs_[i] < r.degs_[i])
return true;
if (r.degs_[i] < degs_[i])
return false;
}
return false;
}
inline bool operator==(const MemoKey& r) const
{
if (nDegs_ != r.nDegs_)
return false;
for (unsigned i=0; i<nDegs_; ++i)
if (degs_[i] != r.degs_[i])
return false;
return true;
}
inline bool operator!=(const MemoKey& r) const
{return !(*this == r);}
private:
unsigned degs_[8];
unsigned nDegs_;
};
void calcMeanXandWeightSums();
void calcRecurrenceCoeffs(OrthoPolyMethod m);
void calcRecurrenceStieltjes();
void calcRecurrenceLanczos();
precise_type monicpoly(unsigned degree, precise_type x) const;
std::pair<precise_type,precise_type> monicInnerProducts(
unsigned degree) const;
precise_type normpoly(unsigned degree, precise_type x) const;
std::pair<precise_type,precise_type> twonormpoly(
unsigned deg1, unsigned deg2, precise_type x) const;
precise_type normpolyprod(const unsigned* degrees, unsigned nDeg,
unsigned degmax, precise_type x) const;
template <typename Numeric>
precise_type normseries(const Numeric* coeffs, unsigned maxdeg,
precise_type x) const;
template <typename Numeric>
void getAllPolys(double x, unsigned maxdeg, Numeric *polyValues) const;
double getCachedAverage(const MemoKey& key) const;
std::vector<MeasurePoint> measure_;
std::vector<Recurrence> rCoeffs_;
mutable std::map<MemoKey,double> cachedAverages_;
precise_type wsum_;
precise_type wsumsq_;
precise_type meanX_;
double minX_;
double maxX_;
unsigned maxdeg_;
bool allWeightsEqual_;
#ifdef SWIG
public:
inline StorablePolySeries1D* makeStorablePolySeries_2(
double xmin, double xmax, const std::vector<double>& coeffs) const
{
CPP11_auto_ptr<StorablePolySeries1D> ptr;
if (coeffs.empty())
ptr = this->makeStorablePolySeries(xmin, xmax);
else
ptr = this->makeStorablePolySeries(xmin, xmax, &coeffs[0],
coeffs.size() - 1);
return ptr.release();
}
inline double sumOfWeights_2() const
{return sumOfWeights();}
inline double sumOfWeightSquares_2() const
{return sumOfWeightSquares();}
inline double meanCoordinate_2() const
{return meanCoordinate();}
inline double empiricalKroneckerDelta_2(const unsigned deg1,
const unsigned deg2) const
{return empiricalKroneckerDelta(deg1, deg2);}
inline Matrix<double> jacobiMatrix_2(const unsigned n) const
{return jacobiMatrix(n);}
inline double integratePoly_2(const unsigned degree, const unsigned power,
const double xmin, const double xmax) const
{return integratePoly(degree, power, xmin, xmax);}
inline double integrateTripleProduct_2(const unsigned deg1, const unsigned deg2,
const unsigned deg3, const double xmin,
const double xmax) const
{return integrateTripleProduct(deg1, deg2, deg3, xmin, xmax);}
inline Matrix<double> tripleProductMatrix_2(
const std::vector<std::pair<unsigned,unsigned> >& pairs,
const unsigned maxdeg, const double xmin, const double xmax) const
{return tripleProductMatrix(pairs, maxdeg, xmin, xmax);}
#endif
};
}
#include "npstat/nm/ContOrthoPoly1D.icc"
#endif // NPSTAT_CONTORTHOPOLY1D_HH_
Index: trunk/npstat/stat/KDE1DHOSymbetaKernel.cc
===================================================================
--- trunk/npstat/stat/KDE1DHOSymbetaKernel.cc (revision 720)
+++ trunk/npstat/stat/KDE1DHOSymbetaKernel.cc (revision 721)
@@ -1,91 +1,92 @@
#include <cassert>
#include <stdexcept>
#include "npstat/nm/ClassicalOrthoPolys1D.hh"
#include "npstat/stat/KDE1DHOSymbetaKernel.hh"
#include "npstat/stat/Distributions1D.hh"
namespace npstat {
KDE1DHOSymbetaKernel::KDE1DHOSymbetaKernel(
const int i_power, const double i_filterDegree)
- : power_(i_power), filterDegree_(i_filterDegree), lastShrinkage_(1.0)
+ : power_(i_power), filterDegree_(i_filterDegree),
+ lastShrinkage_(1.0), poly_(0)
{
if (filterDegree_ < 0.0) throw std::invalid_argument(
"In npstat::KDE1DHOSymbetaKernel constructor: filterDegree"
" argument can not be negative");
if (i_power < 0)
poly_ = new HermiteProbOrthoPoly1D();
else if (i_power == 0)
poly_ = new LegendreOrthoPoly1D();
else
poly_ = new JacobiOrthoPoly1D(i_power, i_power);
maxdeg_ = floor(filterDegree_);
if (static_cast<double>(maxdeg_) != filterDegree_)
{
lastShrinkage_ = sqrt(filterDegree_ - maxdeg_);
++maxdeg_;
}
polyValues_.resize(2U*(maxdeg_ + 1U));
poly_->allpoly(0.0L, &polyValues_[0], maxdeg_);
if (i_power < 0)
{
Gauss1D g(0.0, 1.0);
xmin_ = g.quantile(0.0);
xmax_ = g.quantile(1.0);
}
else
{
xmin_ = poly_->xmin();
xmax_ = poly_->xmax();
}
}
KDE1DHOSymbetaKernel::KDE1DHOSymbetaKernel(const KDE1DHOSymbetaKernel& r)
: poly_(0)
{
copyInternals(r);
}
void KDE1DHOSymbetaKernel::copyInternals(const KDE1DHOSymbetaKernel& r)
{
assert(!poly_);
poly_ = r.poly_->clone();
power_ = r.power_;
maxdeg_ = r.maxdeg_;
filterDegree_ = r.filterDegree_;
xmin_ = r.xmin_;
xmax_ = r.xmax_;
lastShrinkage_ = r.lastShrinkage_;
polyValues_ = r.polyValues_;
}
KDE1DHOSymbetaKernel& KDE1DHOSymbetaKernel::operator=(
const KDE1DHOSymbetaKernel& r)
{
if (this != &r)
{
delete poly_;
poly_ = 0;
copyInternals(r);
}
return *this;
}
double KDE1DHOSymbetaKernel::operator()(const double x) const
{
long double sum = 1.0L;
if (filterDegree_)
{
const long double* poly0 = &polyValues_[0];
long double* polyx = &polyValues_[0] + (maxdeg_ + 1U);
poly_->allpoly(x, polyx, maxdeg_);
for (unsigned i=1U; i<maxdeg_; ++i)
sum += poly0[i]*polyx[i];
sum += poly0[maxdeg_]*polyx[maxdeg_]*lastShrinkage_;
}
return sum*poly_->weight(x);
}
}

File Metadata

Mime Type
text/x-diff
Expires
Sun, Feb 23, 2:46 PM (23 h, 3 m)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4486698
Default Alt Text
(28 KB)

Event Timeline