Page MenuHomeHEPForge
Files F19251508

test_ContOrthoPoly1D.cc
No OneTemporary
Actions

Size
13 KB
Referenced Files
None
Subscribers
None

test_ContOrthoPoly1D.cc
View Options

#include "UnitTest++.h"
#include "test_utils.hh"
#include "npstat/nm/ContOrthoPoly1D.hh"
#include "npstat/nm/ClassicalOrthoPolys1D.hh"
#include "npstat/nm/FejerQuadrature.hh"
#include "npstat/nm/StorablePolySeries1D.hh"
#include "npstat/rng/MersenneTwister.hh"
#include "npstat/stat/Distributions1D.hh"
using namespace npstat;
using namespace std;
inline static int kdelta(const unsigned i, const unsigned j)
{
return i == j ? 1 : 0;
}
namespace {
TEST(ContOrthoPoly1D_orthonormalization)
{
const double eps = 1.0e-15;
const OrthoPolyMethod method[] = {OPOLY_STIELTJES, OPOLY_LANCZOS};
const unsigned nMethods = sizeof(method)/sizeof(method[0]);
const unsigned npoints = 64;
const unsigned maxdeg = 10;
const unsigned ntries = 10;
std::vector<double> points(2U*npoints);
std::vector<ContOrthoPoly1D::MeasurePoint> measure;
measure.reserve(npoints);
std::vector<double> coeffs(maxdeg + 1);
MersenneTwister rng;
for (unsigned imeth=0; imeth<nMethods; ++imeth)
for (unsigned itry=0; itry<ntries; ++itry)
{
const double xmin = rng()*5.0 - 6.0;
const double xmax = rng()*2.0 + 1.0;
const double range = xmax - xmin;
rng.run(&points[0], 2U*npoints, 2U*npoints);
measure.clear();
for (unsigned i=0; i<npoints; ++i)
{
ContOrthoPoly1D::MeasurePoint p(points[2*i]*range + xmin, points[2*i+1]);
measure.push_back(p);
}
ContOrthoPoly1D poly(maxdeg, measure, method[imeth]);
for (unsigned i=0; i<=maxdeg; ++i)
for (unsigned j=0; j<=maxdeg; ++j)
{
const double d = poly.empiricalKroneckerDelta(i, j);
CHECK_CLOSE(static_cast<double>(kdelta(i, j)), d, eps);
}
measure.clear();
for (unsigned i=0; i<npoints; ++i)
{
ContOrthoPoly1D::MeasurePoint p(points[2*i]*range + xmin, 1.0);
measure.push_back(p);
}
ContOrthoPoly1D poly2(maxdeg, measure, method[imeth]);
poly2.weightCoeffs(&coeffs[0], maxdeg);
for (unsigned i=0; i<=maxdeg; ++i)
CHECK_CLOSE(static_cast<double>(kdelta(i, 0)), coeffs[i], eps);
}
}
#ifdef USE_LAPACK_QUAD
TEST(ContOrthoPoly1D_orthonormalization_quad_STIELTJES)
{
const lapack_double eps = FLT128_EPSILON*100.0;
const OrthoPolyMethod method = OPOLY_STIELTJES;
const unsigned npoints = 64;
const unsigned maxdeg = 10;
const unsigned ntries = 5;
const double xmin = -1.0;
const double xmax = 1.0;
const double range = xmax - xmin;
std::vector<double> points(2U*npoints);
std::vector<ContOrthoPoly1D::MeasurePoint> measure;
measure.reserve(npoints);
MersenneTwister rng;
for (unsigned itry=0; itry<ntries; ++itry)
{
rng.run(&points[0], 2U*npoints, 2U*npoints);
measure.clear();
for (unsigned i=0; i<npoints; ++i)
{
ContOrthoPoly1D::MeasurePoint p(points[2*i]*range + xmin, points[2*i+1]);
measure.push_back(p);
}
ContOrthoPoly1D poly(maxdeg, measure, method);
for (unsigned i=0; i<=maxdeg; ++i)
for (unsigned j=0; j<=maxdeg; ++j)
{
const lapack_double d = poly.empiricalKroneckerDelta(i, j);
CHECK_CLOSE(static_cast<lapack_double>(kdelta(i, j)), d, eps);
}
}
}
TEST(ContOrthoPoly1D_orthonormalization_quad_LANCZOS)
{
const lapack_double eps = FLT128_EPSILON*100.0;
const OrthoPolyMethod method = OPOLY_LANCZOS;
const unsigned npoints = 64;
const unsigned maxdeg = 10;
const unsigned ntries = 5;
const double xmin = -1.0;
const double xmax = 1.0;
const double range = xmax - xmin;
std::vector<double> points(2U*npoints);
std::vector<ContOrthoPoly1D::MeasurePoint> measure;
measure.reserve(npoints);
MersenneTwister rng;
for (unsigned itry=0; itry<ntries; ++itry)
{
rng.run(&points[0], 2U*npoints, 2U*npoints);
measure.clear();
for (unsigned i=0; i<npoints; ++i)
{
ContOrthoPoly1D::MeasurePoint p(points[2*i]*range + xmin, points[2*i+1]);
measure.push_back(p);
}
ContOrthoPoly1D poly(maxdeg, measure, method);
for (unsigned i=0; i<=maxdeg; ++i)
for (unsigned j=0; j<=maxdeg; ++j)
{
const lapack_double d = poly.empiricalKroneckerDelta(i, j);
CHECK_CLOSE(static_cast<lapack_double>(kdelta(i, j)), d, eps);
}
}
}
#endif
TEST(ContOrthoPoly1D_weightCoeffs)
{
const double eps = 1.0e-7;
const unsigned maxdeg = 10;
const unsigned ntries = 10;
const unsigned npoints = maxdeg + 1U;
std::vector<double> points(2U*npoints);
std::vector<ContOrthoPoly1D::MeasurePoint> measure;
measure.reserve(npoints);
std::vector<double> coeffs(maxdeg + 1);
MersenneTwister rng;
for (unsigned itry=0; itry<ntries; ++itry)
{
const double xmin = rng()*5.0 - 6.0;
const double xmax = rng()*2.0 + 1.0;
const double range = xmax - xmin;
rng.run(&points[0], 2U*npoints, 2U*npoints);
measure.clear();
for (unsigned i=0; i<npoints; ++i)
{
ContOrthoPoly1D::MeasurePoint p(points[2*i]*range + xmin, points[2*i+1]);
measure.push_back(p);
}
ContOrthoPoly1D poly(maxdeg, measure, OPOLY_LANCZOS);
poly.weightCoeffs(&coeffs[0], maxdeg);
CPP11_auto_ptr<StorablePolySeries1D> poly2(
poly.makeStorablePolySeries(xmin, xmax));
for (unsigned deg=0; deg<=maxdeg; ++deg)
for (unsigned ipow=0; ipow<4; ++ipow)
{
const double i1 = poly.integratePoly(deg, ipow, xmin, xmax);
const double i2 = poly2->integratePoly(deg, ipow);
if (fabs(i2) > 1.0e-10)
CHECK_CLOSE(1.0, i1/i2, 1.0e-10);
}
for (unsigned i=0; i<npoints; ++i)
{
const double x = measure[i].first;
const double w = measure[i].second;
const double fvalue = poly.series(&coeffs[0], maxdeg, x);
const double f2 = poly2->series(&coeffs[0], maxdeg, x);
CHECK_CLOSE(w, fvalue, 1000*eps);
CHECK_CLOSE(fvalue, f2, 10*eps);
}
}
}
TEST(ContOrthoPoly1D_cov8)
{
const double eps = 1.0e-12;
MersenneTwister rng;
const unsigned npoints = 64;
const unsigned maxdeg = 6;
const unsigned ntries = 5;
const unsigned degtries = 100;
std::vector<double> points(npoints);
for (unsigned itry=0; itry<ntries; ++itry)
{
rng.run(&points[0], npoints, npoints);
ContOrthoPoly1D poly(maxdeg, points);
for (unsigned ideg=0; ideg<degtries; ++ideg)
{
unsigned d[8];
bool has0 = true;
while (has0)
{
for (unsigned i=0; i<8; ++i)
d[i] = (maxdeg + 0.99)*rng();
has0 = (d[0] == 0 && d[1] == 0) ||
(d[2] == 0 && d[3] == 0) ||
(d[4] == 0 && d[5] == 0) ||
(d[6] == 0 && d[7] == 0);
}
const double cov8 = poly.cov8(d[0],d[1],d[2],d[3],d[4],d[5],d[6],d[7]);
const double check = poly.slowCov8(d[0],d[1],d[2],d[3],d[4],d[5],d[6],d[7]);
CHECK_CLOSE(check, cov8, eps);
}
}
}
TEST(ContOrthoPoly1D_integrateEWPolyProduct)
{
const double eps = 1.0e-12;
MersenneTwister rng;
const unsigned npoints = 64;
const unsigned maxdeg = 6;
const unsigned ntries = 5;
const unsigned nInteg = 64;
std::vector<double> points(npoints);
Beta1D beta(0.0, 1.0, 2.0, 2.0);
GaussLegendreQuadrature glq(nInteg);
for (unsigned itry=0; itry<ntries; ++itry)
{
rng.run(&points[0], npoints, npoints);
ContOrthoPoly1D poly(maxdeg, points);
const Matrix<double>& mat = poly.extWeightProductAverages(
DensityFunctor1D(beta), glq, maxdeg, 0.0, 1.0);
for (unsigned i=0; i<=maxdeg; ++i)
for (unsigned j=0; j<=maxdeg; ++j)
{
const double integ = poly.integrateEWPolyProduct(
DensityFunctor1D(beta), i, j, 0.0, 1.0, nInteg);
CHECK_CLOSE(integ, mat[i][j], eps);
}
}
}
TEST(ContOrthoPoly1D_epsCovariance)
{
const double eps = 1.0e-15;
const unsigned npoints = 64;
const unsigned maxdeg = 6;
const unsigned ntries = 5;
MersenneTwister rng;
std::vector<double> points(npoints);
for (unsigned itry=0; itry<ntries; ++itry)
{
rng.run(&points[0], npoints, npoints);
ContOrthoPoly1D poly(maxdeg, points);
for (unsigned i=0; i<=maxdeg; ++i)
for (unsigned j=0; j<=i; ++j)
{
CHECK(poly.epsExpectation(i, j, false) == 0.0);
for (unsigned k=0; k<=maxdeg; ++k)
for (unsigned m=0; m<=k; ++m)
{
const double eps1 = poly.empiricalKroneckerCovariance(i, j, k, m);
const double eps2 = poly.epsCovariance(i, j, k, m, false);
CHECK_CLOSE(eps1, eps2, eps);
}
}
}
}
// Test that we can reproduce Jacobi polynomials using ContOrthoPoly1D
TEST(ContOrthoPoly1D_jacobi_1)
{
const double eps = 1.0e-13;
const unsigned alpha = 3;
const unsigned beta = 4;
const unsigned maxdeg = 5;
const unsigned npoints = 64;
JacobiOrthoPoly1D jac(alpha, beta);
OrthoPoly1DWeight weight(jac);
OrthoPoly1DDeg poly3(jac, 3);
const unsigned npt = FejerQuadrature::minimalExactRule(2*maxdeg+alpha+beta);
FejerQuadrature fej(npt);
ContOrthoPoly1D poly(maxdeg, fej.weightedIntegrationPoints(weight, -1.0, 1.0), OPOLY_LANCZOS);
for (unsigned i=0; i<npoints; ++i)
{
const double x = test_rng()*2.0 - 1.0;
for (unsigned deg=0; deg<=maxdeg; ++deg)
CHECK_CLOSE(jac.poly(deg, x), poly.poly(deg, x), eps);
CHECK_CLOSE(jac.poly(3, x), poly3(x), eps);
}
// Test "sampleAverages" method
const unsigned nrandom = 100;
std::vector<double> r(nrandom);
for (unsigned i=0; i<nrandom; ++i)
r[i] = -1.0 + 2.0*test_rng();
std::vector<double> result(maxdeg + 1);
poly.sampleAverages(&r[0], r.size(), &result[0], maxdeg);
Matrix<long double> mat(maxdeg + 1, maxdeg + 1, 0);
std::vector<long double> acc(maxdeg + 1, 0.0L);
std::vector<long double> tmp(maxdeg + 1);
for (unsigned i=0; i<nrandom; ++i)
{
poly.allPolys(r[i], maxdeg, &tmp[0]);
for (unsigned deg=0; deg<=maxdeg; ++deg)
{
acc[deg] += tmp[deg];
for (unsigned k=0; k<=maxdeg; ++k)
{
const double prod = tmp[deg]*tmp[k];
mat[deg][k] += prod;
}
}
}
for (unsigned deg=0; deg<=maxdeg; ++deg)
{
acc[deg] /= nrandom;
CHECK_CLOSE(static_cast<double>(acc[deg]), result[deg], eps);
}
// Test "pairwiseSampleAverages" method
mat /= nrandom;
Matrix<double> averages = poly.sampleProductAverages(
&r[0], r.size(), maxdeg);
Matrix<double> refmat(mat);
CHECK((averages - refmat).maxAbsValue() < eps);
}
TEST(ContOrthoPoly1D_jacobi_2)
{
typedef ContOrthoPoly1D::PreciseQuadrature Quadrature;
const double eps = 1.0e-13;
const unsigned alpha = 3;
const unsigned beta = 4;
const unsigned maxdeg = 5;
const unsigned npoints = 64;
JacobiOrthoPoly1D jac(alpha, beta);
OrthoPoly1DWeight weight(jac);
const unsigned npt = Quadrature::minimalExactRule(2*maxdeg+alpha+beta);
Quadrature glq(npt);
ContOrthoPoly1D poly(maxdeg, glq.weightedIntegrationPoints(weight, -1.0, 1.0), OPOLY_LANCZOS);
for (unsigned i=0; i<npoints; ++i)
{
const double x = test_rng()*2.0 - 1.0;
for (unsigned deg=0; deg<=maxdeg; ++deg)
CHECK_CLOSE(jac.poly(deg, x), poly.poly(deg, x), eps);
}
}
}

File Metadata

Mime Type
text/x-c
Expires
Tue, Sep 30, 6:01 AM (1 d, 14 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
6566411
Default Alt Text
test_ContOrthoPoly1D.cc (13 KB)

Event Timeline

HEPForge.org · .