Index: trunk/tests/test_ClassicalOrthoPolys1D.cc =================================================================== --- trunk/tests/test_ClassicalOrthoPolys1D.cc (revision 885) +++ trunk/tests/test_ClassicalOrthoPolys1D.cc (revision 886) @@ -1,498 +1,540 @@ #include #include #include #include #include "UnitTest++.h" #include "test_utils.hh" #include "npstat/nm/ClassicalOrthoPolys1D.hh" #include "npstat/nm/GaussLegendreQuadrature.hh" #include "npstat/nm/FejerQuadrature.hh" #include "npstat/nm/StorablePolySeries1D.hh" #include "npstat/nm/ClassicalOrthoPoly1DFromWeight.hh" #include "npstat/nm/opsRootsFromJacobiMatrix.hh" #include "npstat/rng/MersenneTwister.hh" #include "npstat/stat/Distributions1D.hh" #include "npstat/stat/DensityOrthoPoly1D.hh" using namespace npstat; using namespace std; namespace { + template + void runIntegralCoeffsTest() + { + const double eps = 1.0e-12; + + Poly poly; + const double xmin = poly.xmin(); + const double xmax = poly.xmax(); + GaussLegendreQuadrature glq(10); + + for (unsigned ideg=0; ideg<=7; ++ideg) + { + std::vector coeffs(ideg+1U); + for (unsigned i=0; i<=ideg; ++i) + coeffs[i] = test_rng() - 0.5; + + std::vector intcoeffs(ideg+2U); + poly.integralCoeffs(&coeffs[0], ideg, &intcoeffs[0]); + + CPP11_auto_ptr pser = poly.makeStorablePolySeries( + ideg, &coeffs[0], ideg); + + for (unsigned itry=0; itry<10; ++itry) + { + const double y = xmin + (xmax - xmin)*test_rng(); + const double numInteg = glq.integrate(*pser, xmin, y); + const double serInteg = poly.series(&intcoeffs[0], ideg+1U, y); + CHECK_CLOSE(numInteg, serInteg, eps); + } + } + } + TEST(AbsOrthoPoly1D_cloning) { const unsigned ntries = 100U; const unsigned degs[3] = {3, 4, 5}; AbsClassicalOrthoPoly1D *poly1 = 0, *poly2 = 0; { HermiteProbOrthoPoly1D He; const AbsClassicalOrthoPoly1D& ref(He); poly1 = ref.clone(); poly2 = ref.clone(); const double x = test_rng()*10 - 5; for (unsigned ideg=0; idegpoly(degs[ideg], x), He.poly(degs[ideg], x)); } for (unsigned i=0; ipoly(degs[ideg], x), poly2->poly(degs[ideg], x)); CHECK_EQUAL(poly1->weight(x), poly2->weight(x)); } } delete poly2; delete poly1; } TEST(LegendreOrthoPoly1D_values) { const double eps = 1.0e-14; const LegendreOrthoPoly1D leg; CPP11_auto_ptr leg2(leg.makeStorablePolySeries(8)); CHECK(leg.weight(0.123) == 0.5); CHECK(leg.weight(10.0) == 0.0); unsigned deg = 0; CHECK_CLOSE(1.0, leg.poly(deg, 0.3), eps); CHECK_CLOSE(1.0, leg2->poly(deg, 0.3), eps); deg = 3; CHECK_CLOSE(-71/343.0L*sqrtl(2*deg+1.0L), leg.poly(deg, 1/7.0L), eps); CHECK_CLOSE(-71/343.0L*sqrtl(2*deg+1.0L), leg2->poly(deg, 1/7.0L), eps); deg = 7; CHECK_CLOSE(-326569/1771561.0L*sqrtl(2*deg+1.0L), leg.poly(deg, 1/11.0L), eps); CHECK_CLOSE(-326569/1771561.0L*sqrtl(2*deg+1.0L), leg2->poly(deg, 1/11.0L), eps); GaussLegendreQuadrature glq(10); for (unsigned i=0; i<5; ++i) for (unsigned j=0; j<5; ++j) CHECK_CLOSE((i == j ? 1.0 : 0.0), leg.empiricalKroneckerDelta(glq, i, j), eps); const unsigned maxdeg = 8; double coeffs[maxdeg+1]; for (unsigned i=0; i<=maxdeg; ++i) coeffs[i] = test_rng(); leg2->setCoeffs(coeffs, maxdeg); for (unsigned i=0; i<=maxdeg; ++i) { const double x = test_rng()*2 - 1; CHECK_CLOSE(leg.series(coeffs, maxdeg, x), (*leg2)(x), eps); } CHECK_CLOSE(3374/715.0L, leg.integratePoly(3, 4), eps); CHECK_CLOSE(3374/715.0L, leg2->integratePoly(3, 4), eps); unsigned degs[3] = {3, 4, 5}; CHECK_CLOSE(1.051943782704350297, leg.jointIntegral(degs, 3), eps); std::ostringstream os; CHECK(leg2->classId().write(os)); CHECK(leg2->write(os)); std::istringstream is(os.str()); gs::ClassId clid(is, 1); StorablePolySeries1D* rb = StorablePolySeries1D::read(clid, is); assert(rb); CHECK(*rb == *leg2); delete rb; } TEST(LegendreOrthoPoly1D_roots) { const double eps = 1.0e-14; const unsigned maxdeg = 6; const double roots_1[] = {0.0}; const double roots_2[] = {-0.57735026918962576, 0.57735026918962576}; const double roots_3[] = {-0.77459666924148338, 0.0, 0.77459666924148338}; const double roots_4[] = {-0.86113631159405258, -0.33998104358485626, 0.33998104358485626, 0.86113631159405258}; const double roots_5[] = {-0.90617984593866399, -0.53846931010568309, 0.0, 0.53846931010568309, 0.90617984593866399}; const double roots_6[] = {-0.93246951420315203, -0.66120938646626451, -0.23861918608319691, 0.23861918608319691, 0.66120938646626451, 0.93246951420315203}; const double* roots[maxdeg] = {roots_1, roots_2, roots_3, roots_4, roots_5, roots_6}; const LegendreOrthoPoly1D leg; for (unsigned deg=1; deg<=maxdeg; ++deg) { const double* knownRoots = roots[deg-1U]; const std::vector& jacobyRoots = opsRootsFromJacobiMatrix(leg, deg); for (unsigned i=0; i(); + } + + TEST(ShiftedLegendreOrthoPoly1D_integralCoeffs) + { + runIntegralCoeffsTest(); + } + TEST(JacobiOrthoPoly1D_derivativeCoeffs) { const unsigned deg = 5; const double h = 1.0e-5; double coeffs[deg+1U], deriv[deg]; for (unsigned itry=0; itry<10; ++itry) { const unsigned alpha = test_rng()*6.0; const unsigned beta = test_rng()*6.0; JacobiOrthoPoly1D jac(alpha, beta); JacobiOrthoPoly1D jacder(alpha+1U, beta+1U); for (unsigned ideg=0; ideg<=deg; ++ideg) coeffs[ideg] = test_rng()*2.0 - 1.0; jac.derivativeCoeffs(coeffs, deg, deriv); for (unsigned ix = 0; ix<10; ++ix) { const double x = test_rng()*1.8 - 0.9; const volatile double xplus = x + h; const volatile double xminus = x - h; const double splus = jac.series(coeffs, deg, xplus); const double sminus = jac.series(coeffs, deg, xminus); const double numderi = (splus - sminus)/(xplus - xminus); const double deri = jacder.series(deriv, deg-1U, x); CHECK_CLOSE(numderi, deri, 1.0e-7*(std::abs(deri) + 1.0)); } } } TEST(JacobiOrthoPoly1D_values) { const double eps = 1.0e-15; const unsigned densMaxDeg = 10; JacobiOrthoPoly1D jac(3, 4); const Beta1D b1(-1.0, 2.0, 5, 4); const DensityOrthoPoly1D dpoly_00(b1, densMaxDeg, 2*densMaxDeg+5+4, OPOLY_STIELTJES); const DensityOrthoPoly1D dpoly_01(b1, densMaxDeg, 2*densMaxDeg+5+4, OPOLY_LANCZOS); OrthoPoly1DWeight kernel(jac); const ClassicalOrthoPoly1DFromWeight wpoly_00( kernel, densMaxDeg, 2*densMaxDeg+5+4, -1, 1, OPOLY_STIELTJES); const ClassicalOrthoPoly1DFromWeight wpoly_01( kernel, densMaxDeg, 2*densMaxDeg+5+4, -1, 1, OPOLY_LANCZOS); CHECK_EQUAL(densMaxDeg, dpoly_00.maxDegree()); CHECK_EQUAL(densMaxDeg, dpoly_01.maxDegree()); CHECK_EQUAL(densMaxDeg, wpoly_00.maxDegree()); CHECK_EQUAL(densMaxDeg, wpoly_01.maxDegree()); long double x = 1/3.0L; CHECK_CLOSE(35/32.0L*powl(1-x,3)*powl(1+x,4), jac.weight(x), eps); for (unsigned i=0; i<100; ++i) { const double x = test_rng()*2.0 - 1.0; const double w = jac.weight(x); CHECK_CLOSE(w, dpoly_00.weight(x), eps); CHECK_CLOSE(w, dpoly_01.weight(x), eps); CHECK_CLOSE(w, wpoly_00.weight(x), eps); CHECK_CLOSE(w, wpoly_01.weight(x), eps); for (unsigned deg=0; deg<=densMaxDeg; ++deg) { const double poly = jac.poly(deg, x); CHECK_CLOSE(poly, dpoly_00.poly(deg, x), eps*20*deg); CHECK_CLOSE(poly, dpoly_01.poly(deg, x), eps*20*deg); CHECK_CLOSE(poly, wpoly_00.poly(deg, x), eps*20*deg); CHECK_CLOSE(poly, wpoly_01.poly(deg, x), eps*20*deg); } } GaussLegendreQuadrature ptest(128); std::vector aves(densMaxDeg + 1U); const Matrix& m00 = dpoly_00.extWeightProductAverages( DensityFunctor1D(b1), ptest, densMaxDeg); dpoly_00.extWeightAverages(DensityFunctor1D(b1), ptest, &aves[0], densMaxDeg); CHECK_EQUAL(densMaxDeg+1U, m00.nRows()); CHECK_EQUAL(densMaxDeg+1U, m00.nColumns()); Matrix diag(densMaxDeg+1U, densMaxDeg+1U, 1); CHECK((m00 - diag).maxAbsValue() < eps); for (unsigned i=0; i<=densMaxDeg; ++i) CHECK(abs(m00[0][i] - aves[i]) < eps); const Matrix& m01 = dpoly_01.extWeightProductAverages( DensityFunctor1D(b1), ptest, densMaxDeg); dpoly_01.extWeightAverages(DensityFunctor1D(b1), ptest, &aves[0], densMaxDeg); CHECK((m01 - diag).maxAbsValue() < eps); for (unsigned i=0; i<=densMaxDeg; ++i) CHECK(abs(m01[0][i] - aves[i]) < eps); GaussLegendreQuadrature glq(10); for (unsigned i=0; i<5; ++i) for (unsigned j=0; j<5; ++j) { CHECK_CLOSE((i == j ? 1.0 : 0.0), jac.empiricalKroneckerDelta(glq, i, j), eps); CHECK_CLOSE((i == j ? 1.0 : 0.0), dpoly_00.empiricalKroneckerDelta(glq, i, j), eps); CHECK_CLOSE((i == j ? 1.0 : 0.0), dpoly_01.empiricalKroneckerDelta(glq, i, j), eps); CHECK_CLOSE((i == j ? 1.0 : 0.0), wpoly_00.empiricalKroneckerDelta(glq, i, j), eps); CHECK_CLOSE((i == j ? 1.0 : 0.0), wpoly_01.empiricalKroneckerDelta(glq, i, j), eps); unsigned degs[2]; degs[0] = i; degs[1] = j; CHECK_CLOSE((i == j ? 1.0 : 0.0), jac.jointAverage(glq, degs, 2), eps); CHECK_CLOSE((i == j ? 1.0 : 0.0), dpoly_00.jointAverage(glq, degs, 2), eps); CHECK_CLOSE((i == j ? 1.0 : 0.0), dpoly_01.jointAverage(glq, degs, 2), eps); CHECK_CLOSE((i == j ? 1.0 : 0.0), wpoly_00.jointAverage(glq, degs, 2), eps); CHECK_CLOSE((i == j ? 1.0 : 0.0), wpoly_01.jointAverage(glq, degs, 2), eps); } JacobiOrthoPoly1D jac2(2, 5); const Beta1D b2(-1.0, 2.0, 6, 3); const DensityOrthoPoly1D dpoly_20(b2, densMaxDeg, 2*densMaxDeg+6+3, OPOLY_STIELTJES); const DensityOrthoPoly1D dpoly_21(b2, densMaxDeg, 2*densMaxDeg+6+3, OPOLY_LANCZOS); OrthoPoly1DWeight kernel2(jac2); const ClassicalOrthoPoly1DFromWeight wpoly_20( kernel2, densMaxDeg, 2*densMaxDeg+6+3, -1, 1, OPOLY_STIELTJES); const ClassicalOrthoPoly1DFromWeight wpoly_21( kernel2, densMaxDeg, 2*densMaxDeg+6+3, -1, 1, OPOLY_LANCZOS); for (unsigned i=0; i<100; ++i) { const double x = test_rng()*2.0 - 1.0; const double w = jac2.weight(x); CHECK_CLOSE(w, dpoly_20.weight(x), eps); CHECK_CLOSE(w, dpoly_21.weight(x), eps); CHECK_CLOSE(w, wpoly_20.weight(x), eps); CHECK_CLOSE(w, wpoly_21.weight(x), eps); for (unsigned deg=0; deg<=densMaxDeg; ++deg) { const double poly = jac2.poly(deg, x); CHECK_CLOSE(poly, dpoly_20.poly(deg, x), eps*100*deg); CHECK_CLOSE(poly, dpoly_21.poly(deg, x), eps*100*deg); CHECK_CLOSE(poly, wpoly_20.poly(deg, x), eps*100*deg); CHECK_CLOSE(poly, wpoly_21.poly(deg, x), eps*100*deg); } } for (unsigned i=0; i<5; ++i) for (unsigned j=0; j<5; ++j) { CHECK_CLOSE((i == j ? 1.0 : 0.0), jac2.empiricalKroneckerDelta(glq, i, j), eps); CHECK_CLOSE((i == j ? 1.0 : 0.0), dpoly_20.empiricalKroneckerDelta(glq, i, j), eps); CHECK_CLOSE((i == j ? 1.0 : 0.0), dpoly_21.empiricalKroneckerDelta(glq, i, j), eps); CHECK_CLOSE((i == j ? 1.0 : 0.0), wpoly_20.empiricalKroneckerDelta(glq, i, j), eps); CHECK_CLOSE((i == j ? 1.0 : 0.0), wpoly_21.empiricalKroneckerDelta(glq, i, j), eps); unsigned degs[2]; degs[0] = i; degs[1] = j; CHECK_CLOSE((i == j ? 1.0 : 0.0), jac2.jointAverage(glq, degs, 2), eps); CHECK_CLOSE((i == j ? 1.0 : 0.0), dpoly_20.jointAverage(glq, degs, 2), eps); CHECK_CLOSE((i == j ? 1.0 : 0.0), dpoly_21.jointAverage(glq, degs, 2), eps); CHECK_CLOSE((i == j ? 1.0 : 0.0), wpoly_20.jointAverage(glq, degs, 2), eps); CHECK_CLOSE((i == j ? 1.0 : 0.0), wpoly_21.jointAverage(glq, degs, 2), eps); } JacobiOrthoPoly1D jac3(0, 0); CHECK_CLOSE(0.5, jac3.weight(0.123), eps); CHECK(jac3.weight(10.0) == 0.0); unsigned deg = 7; CHECK_CLOSE(-326569/1771561.0L*sqrtl(2*deg+1.0L), jac3.poly(deg, 1/11.0L), eps); MersenneTwister rng; const unsigned npoints = 64; const unsigned maxdeg = 7; const unsigned ntries = 5; std::vector points(npoints); for (unsigned itry=0; itry avemat(jac.sampleProductAverages(&points[0], npoints, maxdeg)); std::vector averages(maxdeg + 1); jac.sampleAverages(&points[0], npoints, &averages[0], maxdeg); CHECK_EQUAL(maxdeg + 1, avemat.nRows()); CHECK_EQUAL(maxdeg + 1, avemat.nColumns()); for (unsigned i=0; i<=maxdeg; ++i) { CHECK_CLOSE(averages[i], avemat[i][0], eps); CHECK_CLOSE(averages[i], avemat[0][i], eps); for (unsigned j=0; j<=maxdeg; ++j) { unsigned degs[2]; degs[0] = i; degs[1] = j; const double ave = jac.jointSampleAverage(&points[0], npoints, degs, 2); CHECK_CLOSE(ave, avemat[i][j], eps); CHECK_CLOSE(ave, avemat[j][i], eps); } } // Test "twoJointSampleAverages" method const double degmax = maxdeg + 0.99; unsigned degs0[3]; for (unsigned i=0; i<3; ++i) degs0[i] = rng()*degmax; unsigned degs1[4]; for (unsigned i=0; i<4; ++i) degs1[i] = rng()*degmax; const double ave0 = jac.jointSampleAverage(&points[0], npoints, degs0, 3); const double ave1 = jac.jointSampleAverage(&points[0], npoints, degs1, 4); const std::pair& ave2 = jac.twoJointSampleAverages( &points[0], npoints, degs0, 3, degs1, 4); CHECK_CLOSE(ave0, ave2.first, eps); CHECK_CLOSE(ave1, ave2.second, eps); } } TEST(ChebyshevOrthoPoly2nd_values) { const double epsVal = 1.0e-15; const double epsKron = 1.0e-7; ChebyshevOrthoPoly2nd cheb; unsigned deg = 3; CHECK_CLOSE(-0.54810495626822157434L, cheb.poly(deg, 1/7.0L), epsVal); deg = 7; CHECK_CLOSE(-0.66835314371696127673L, cheb.poly(deg, 1/11.0L), epsVal); GaussLegendreQuadrature glq(1024); for (unsigned i=0; i<5; ++i) for (unsigned j=0; j<5; ++j) CHECK_CLOSE((i == j ? 1.0 : 0.0), cheb.empiricalKroneckerDelta(glq, i, j), epsKron); } TEST(ChebyshevOrthoPoly1st_values) { const double epsVal = 1.0e-15; const double epsKron = 1.0e-3; const double norm = 1.0/sqrt(2.0); ChebyshevOrthoPoly1st cheb; unsigned deg = 0; CHECK_CLOSE(1.0, cheb.poly(deg, 0.3), epsVal); deg = 3; CHECK_CLOSE(-0.41690962099125364431L/norm, cheb.poly(deg, 1/7.0L), epsVal); deg = 7; CHECK_CLOSE(-0.59498215518301758629L/norm, cheb.poly(deg, 1/11.0L), epsVal); FejerQuadrature fej(1000); for (unsigned i=0; i<5; ++i) for (unsigned j=0; j<5; ++j) CHECK_CLOSE((i == j ? 1.0 : 0.0), cheb.empiricalKroneckerDelta(fej, i, j), epsKron); } TEST(HermiteProbOrthoPoly1D_values) { const double eps = 1.0e-15; const double epsKron = 1.0e-7; HermiteProbOrthoPoly1D He; long double x = 1/3.0L; CHECK_CLOSE(He.weight(x), exp(-x*x/2)/sqrt(2.0*M_PI), eps); CHECK_CLOSE(He.poly(0, x), 1.0L, eps); CHECK_CLOSE(He.poly(1, x), 1/3.0L, eps); CHECK_CLOSE(He.poly(2, x), -0.62853936105470891058L, eps); CHECK_CLOSE(He.poly(3, x), -0.39312798340964586761L, eps); CHECK_CLOSE(He.poly(4, x), 0.47880972338354304389, eps); GaussLegendreQuadrature glq(1024); for (unsigned i=0; i<5; ++i) for (unsigned j=0; j<5; ++j) CHECK_CLOSE((i == j ? 1.0 : 0.0), He.empiricalKroneckerDelta(glq, i, j), epsKron); } TEST(HermiteProbOrthoPoly1D_dens_1) { const double eps = 1.0e-15; const double epsKron = 1.0e-7; Gauss1D g(0.0, 1.0); const DensityOrthoPoly1D He(g, 5, 1000, OPOLY_STIELTJES); long double x = 1/3.0L; CHECK_CLOSE(He.weight(x), exp(-x*x/2)/sqrt(2.0*M_PI), eps); CHECK_CLOSE(He.poly(0, x), 1.0L, eps); CHECK_CLOSE(He.poly(1, x), 1/3.0L, eps); CHECK_CLOSE(He.poly(2, x), -0.62853936105470891058L, eps); CHECK_CLOSE(He.poly(3, x), -0.39312798340964586761L, eps); CHECK_CLOSE(He.poly(4, x), 0.47880972338354304389, eps); GaussLegendreQuadrature glq(1024); for (unsigned i=0; i<5; ++i) for (unsigned j=0; j<5; ++j) CHECK_CLOSE((i == j ? 1.0 : 0.0), He.empiricalKroneckerDelta(glq, i, j), epsKron); } TEST(HermiteProbOrthoPoly1D_dens_2) { const double eps = 1.0e-15; const double epsKron = 1.0e-7; Gauss1D g(0.0, 1.0); const DensityOrthoPoly1D He(g, 5, 1000, OPOLY_LANCZOS); long double x = 1/3.0L; CHECK_CLOSE(He.weight(x), exp(-x*x/2)/sqrt(2.0*M_PI), eps); CHECK_CLOSE(He.poly(0, x), 1.0L, eps); CHECK_CLOSE(He.poly(1, x), 1/3.0L, eps); CHECK_CLOSE(He.poly(2, x), -0.62853936105470891058L, eps); CHECK_CLOSE(He.poly(3, x), -0.39312798340964586761L, eps); CHECK_CLOSE(He.poly(4, x), 0.47880972338354304389, eps); GaussLegendreQuadrature glq(1024); for (unsigned i=0; i<5; ++i) for (unsigned j=0; j<5; ++j) CHECK_CLOSE((i == j ? 1.0 : 0.0), He.empiricalKroneckerDelta(glq, i, j), epsKron); } TEST(AbsClassicalOrthoPoly1D_sampleCoeffVars) { typedef std::pair Point; const unsigned maxdeg = 10; const unsigned nPseudo = 100; const unsigned perPseudo = 100; const double loc = 2.0; const double scale = 8.0; ShiftedLegendreOrthoPoly1D leg; std::vector points(perPseudo); std::vector wPoints(perPseudo); std::vector coeffs(maxdeg + 1U, 100000.0); std::vector vars(maxdeg + 1U, -100000.0); for (unsigned ips=0; ips 0.0); leg.weightedSampleCoeffs(&wPoints[0], perPseudo, loc, scale, &coeffs[0], maxdeg); leg.weightedSampleCoeffVars(&wPoints[0], perPseudo, loc, scale, &coeffs[0], maxdeg, &vars[0]); for (unsigned i=1; i<=maxdeg; ++i) CHECK(vars[i] > 0.0); } } } Index: trunk/npstat/swig/nm/ClassicalOrthoPolys1D.i =================================================================== --- trunk/npstat/swig/nm/ClassicalOrthoPolys1D.i (revision 885) +++ trunk/npstat/swig/nm/ClassicalOrthoPolys1D.i (revision 886) @@ -1,34 +1,38 @@ %include nm/AbsClassicalOrthoPoly1D.i %include numpy.i %{ #include "npstat/nm/ClassicalOrthoPolys1D.hh" %} %apply (double* IN_ARRAY1, int DIM1) { (const double *c, const unsigned degreep1) }; %feature("notabstract") npstat::LegendreOrthoPoly1D; %ignore npstat::LegendreOrthoPoly1D::weight; +%ignore npstat::LegendreOrthoPoly1D::integralCoeffs; +%rename(integralCoeffs) npstat::LegendreOrthoPoly1D::integralCoeffs2; %feature("notabstract") npstat::ShiftedLegendreOrthoPoly1D; %ignore npstat::ShiftedLegendreOrthoPoly1D::weight; +%ignore npstat::ShiftedLegendreOrthoPoly1D::integralCoeffs; +%rename(integralCoeffs) npstat::ShiftedLegendreOrthoPoly1D::integralCoeffs2; %feature("notabstract") npstat::JacobiOrthoPoly1D; %ignore npstat::JacobiOrthoPoly1D::weight; %ignore npstat::JacobiOrthoPoly1D::derivativeCoeffs; %rename(derivativeCoeffs) npstat::JacobiOrthoPoly1D::derivativeCoeffs2; %feature("notabstract") npstat::ChebyshevOrthoPoly1st; %ignore npstat::ChebyshevOrthoPoly1st::weight; %feature("notabstract") npstat::ChebyshevOrthoPoly2nd; %ignore npstat::ChebyshevOrthoPoly2nd::weight; %feature("notabstract") npstat::HermiteProbOrthoPoly1D; %ignore npstat::HermiteProbOrthoPoly1D::weight; %include "npstat/nm/ClassicalOrthoPolys1D.hh" %clear (const double *c, const unsigned degreep1); Index: trunk/npstat/nm/ClassicalOrthoPolys1D.cc =================================================================== --- trunk/npstat/nm/ClassicalOrthoPolys1D.cc (revision 885) +++ trunk/npstat/nm/ClassicalOrthoPolys1D.cc (revision 886) @@ -1,115 +1,159 @@ +#include #include #include +#include #include "npstat/nm/ClassicalOrthoPolys1D.hh" #include "npstat/nm/SpecialFunctions.hh" namespace npstat { + void LegendreOrthoPoly1D::integralCoeffs( + const double* coeffs, const unsigned degree, + double* integralCoeffs) + { + const double sqr3 = 1.7320508075688773; + + assert(coeffs); + assert(integralCoeffs); + + integralCoeffs[0] = coeffs[0]; + integralCoeffs[1] = coeffs[0]/sqr3; + std::fill(integralCoeffs+2U, integralCoeffs+(degree+2U), 0.0); + + for (unsigned deg=1U; deg<=degree; ++deg) + { + const double denom = sqrt(2.0*deg + 1); + integralCoeffs[deg+1U] += coeffs[deg]/sqrt(2.0*deg + 3)/denom; + integralCoeffs[deg-1U] -= coeffs[deg]/sqrt(2.0*deg - 1)/denom; + } + } + + void ShiftedLegendreOrthoPoly1D::integralCoeffs( + const double* coeffs, const unsigned degree, + double* integralCoeffs) + { + const double twosqr3 = 3.4641016151377546; + + assert(coeffs); + assert(integralCoeffs); + + integralCoeffs[0] = coeffs[0]/2.0; + integralCoeffs[1] = coeffs[0]/twosqr3; + std::fill(integralCoeffs+2U, integralCoeffs+(degree+2U), 0.0); + + for (unsigned deg=1U; deg<=degree; ++deg) + { + const double denom = 2.0*sqrt(2.0*deg + 1); + integralCoeffs[deg+1U] += coeffs[deg]/sqrt(2.0*deg + 3)/denom; + integralCoeffs[deg-1U] -= coeffs[deg]/sqrt(2.0*deg - 1)/denom; + } + } + JacobiOrthoPoly1D::JacobiOrthoPoly1D(const double a, const double b) : alpha_(a), beta_(b) { if (!(alpha_ > -1.0 && beta_ > -1.0)) throw std::invalid_argument( "In npstat::JacobiOrthoPoly1D constructor: invalid arguments"); norm_ = powl(2.0L, alpha_+beta_+1)*Gamma(alpha_+1)*Gamma(beta_+1)/ Gamma(alpha_+beta_+2); } long double JacobiOrthoPoly1D::weight(const long double x) const { if (x < -1.0L || x > 1.0L) return 0.0L; else { if (x == 1.0L && alpha_ < 0.0L) throw std::invalid_argument( "In npstat::JacobiOrthoPoly1D::weight: invalid x=1 argument"); if (x == -1.0L && beta_ < 0.0L) throw std::invalid_argument( "In npstat::JacobiOrthoPoly1D::weight: invalid x=-1 argument"); // Don't call "powl" if we can avoid it because this is a slow function long double p; if (alpha_ == beta_) { const long double tmp = 1.0L - x*x; // Fourth power is a special case, often used in density estimation if (alpha_ == 4.0) { const long double tmpsq = tmp*tmp; p = tmpsq*tmpsq; } else if (alpha_ == 0.0) p = 1.0L; else if (alpha_ == 1.0) p = tmp; else if (alpha_ == 2.0) p = tmp*tmp; else if (alpha_ == 3.0) p = tmp*tmp*tmp; else if (alpha_ == 5.0) { const long double tmpsq = tmp*tmp; p = tmpsq*tmpsq*tmp; } else if (alpha_ == 6.0) { const long double tmpsq = tmp*tmp; p = tmpsq*tmpsq*tmpsq; } else p = powl(tmp, alpha_); } else p = powl(1.0L - x, alpha_)*powl(1.0L + x, beta_); return p/norm_; } } std::pair JacobiOrthoPoly1D::recurrenceCoeffs(const unsigned k) const { long double a, b = 1.0L; if (k) { const long double tmp = 2*k + alpha_ + beta_; a = (beta_ - alpha_)*(beta_ + alpha_)/tmp/(tmp + 2); if (k == 1U) b = 4*k*(k+alpha_)*(k+beta_)/tmp/tmp/(tmp + 1); else b = 4*k*(k+alpha_)*(k+beta_)*(k+alpha_+beta_)/ tmp/tmp/(tmp + 1)/(tmp - 1); } else a = (beta_ - alpha_)/(beta_ + alpha_ + 2); return std::pair(a, sqrtl(b)); } void JacobiOrthoPoly1D::derivativeCoeffs( const double* coeffs, const unsigned maxdeg, double* derivativeCoeffs) const { if (maxdeg) { assert(coeffs); assert(derivativeCoeffs); assert(derivativeCoeffs != coeffs); const double aplusbplusone = alpha_ + beta_+ 1.0; const double commonf = 0.25/(alpha_+1.0)/(beta_+1.0)*(aplusbplusone+1.0)*(aplusbplusone+2.0); for (unsigned ideg=maxdeg; ideg>0U; --ideg) derivativeCoeffs[ideg-1U] = coeffs[ideg]*sqrt(ideg*(aplusbplusone + ideg)*commonf); } } long double HermiteProbOrthoPoly1D::weight(const long double x) const { const long double norm = 0.39894228040143267793994606L; // 1/sqrt(2 Pi) return norm*expl(-0.5*x*x); } std::pair HermiteProbOrthoPoly1D::recurrenceCoeffs(const unsigned k) const { return std::pair(0.0L, sqrtl(k)); } } Index: trunk/npstat/nm/AbsClassicalOrthoPoly1D.hh =================================================================== --- trunk/npstat/nm/AbsClassicalOrthoPoly1D.hh (revision 885) +++ trunk/npstat/nm/AbsClassicalOrthoPoly1D.hh (revision 886) @@ -1,630 +1,630 @@ #ifndef NPSTAT_ABSCLASSICALORTHOPOLY1D_HH_ #define NPSTAT_ABSCLASSICALORTHOPOLY1D_HH_ /*! // \file AbsClassicalOrthoPoly1D.hh // // \brief Base class for classical continuous orthonormal polynomials // // Author: I. Volobouev // // May 2017 */ #include #include #include #include #include #include "geners/CPP11_auto_ptr.hh" #include "npstat/nm/Matrix.hh" #include "npstat/nm/SimpleFunctors.hh" #include "npstat/nm/AbsIntervalQuadrature1D.hh" namespace npstat { class StorablePolySeries1D; class OrthoPoly1DDeg; class OrthoPoly1DWeight; class AbsClassicalOrthoPoly1D { public: typedef OrthoPoly1DDeg degree_functor; typedef OrthoPoly1DWeight weight_functor; inline virtual ~AbsClassicalOrthoPoly1D() {} /** Virtual copy constructor */ virtual AbsClassicalOrthoPoly1D* clone() const = 0; /** The weight function should be normalized */ virtual long double weight(long double x) const = 0; //@{ /** Support of the polynomial */ virtual double xmin() const = 0; virtual double xmax() const = 0; //@} /** Maximum polynomial degree supported */ inline virtual unsigned maxDegree() const {return UINT_MAX - 1U;} /** // Generate principal minor of order n of the Jacobi matrix. // n must not exceed the output of "maxDegree()" method. */ Matrix jacobiMatrix(unsigned n) const; /** Polynomial values */ long double poly(unsigned deg, long double x) const; /** // Values of two orthonormal polynomials. // Faster than calling "poly" two times. */ std::pair twopoly( unsigned deg1, unsigned deg2, long double x) const; /** // Values of all orthonormal polynomials up to some degree. // Faster than calling "poly" multiple times. The size of // the "values" array should be at least maxdeg + 1. */ void allpoly(long double x, long double* values, unsigned maxdeg) const; /** Polynomial series */ double series(const double* coeffs, unsigned maxdeg, double x) const; - /** Integral of the series */ + /** Integral of the series from xmin() to xmax() */ double integrateSeries(const double* coeffs, unsigned maxdeg) 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 void calculateCoeffs(const Functor& fcn, const Quadrature& quad, double* coeffs, unsigned maxdeg) const; /** // Build the coefficients of the orthogonal polynomial series // for the given sample of points (empirical density 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 void sampleCoeffs(const Numeric* coords, unsigned long lenCoords, double* coeffs, unsigned maxdeg) const; /** // Estimate the variances of the coefficients returned by // the "sampleCoeffs" function. The "coeffs" array should be // at least maxdeg + 1 elements long and should be filled // by a previous call to "sampleCoeffs". The "variances" array // should have at least maxdeg + 1 elements. It will contain // the respective variances upon return. */ template void sampleCoeffVars(const Numeric* coords, unsigned long lenCoords, const double* coeffs, unsigned maxdeg, double* variances) const; /** // Estimate the covariances of the coefficients returned by // the "sampleCoeffs" function. The "coeffs" array should be // at least maxdeg + 1 elements long and should be filled // by a previous call to "sampleCoeffs". The returned matrix // will be dimensioned (maxdeg + 1) x (maxdeg + 1). */ template npstat::Matrix sampleCoeffCovariance(const Numeric* coords, unsigned long lenCoords, const double* coeffs, unsigned maxdeg) const; /** // Build the coefficients of the orthogonal polynomial series // for the given sample of points (empirical density 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). // Before calculating the coefficients, the coordinates are shifted // and scaled according to x_new = (x_original - location)/scale. // The resulting coefficients are also divided by scale. */ template void sampleCoeffs(const Numeric* coords, unsigned long lenCoords, Numeric location, Numeric scale, double* coeffs, unsigned maxdeg) const; /** // Estimate the variances of the coefficients returned by // the "sampleCoeffs" function. The "coeffs" array should be // at least maxdeg + 1 elements long and should be filled // by a previous call to "sampleCoeffs" with the same location // and scale. The "variances" array should have at least maxdeg + 1 // elements. It will contain the respective variances upon return. */ template void sampleCoeffVars(const Numeric* coords, unsigned long lenCoords, Numeric location, Numeric scale, const double* coeffs, unsigned maxdeg, double* variances) const; /** // Estimate the covariances of the coefficients returned by // the "sampleCoeffs" function. The "coeffs" array should be // at least maxdeg + 1 elements long and should be filled // by a previous call to "sampleCoeffs" with the same location // and scale. The returned matrix will be dimensioned // (maxdeg + 1) x (maxdeg + 1). */ template npstat::Matrix sampleCoeffCovariance(const Numeric* coords, unsigned long lenCoords, Numeric location, Numeric scale, const double* coeffs, unsigned maxdeg) const; /** // Build the coefficients of the orthogonal polynomial series for // the given sample of weighted points (empirical density function). // The first element of the pair will be treated as the coordinate // and the second element will be treated as weight. Weights must // not be negative. 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 void weightedSampleCoeffs(const Pair* points, unsigned long numPoints, double* coeffs, unsigned maxdeg) const; /** // Estimate the variances of the coefficients returned by the // "weightedSampleCoeffs" function. The "coeffs" array should // be at least maxdeg + 1 elements long and should be filled // by a previous call to "weightedSampleCoeffs". The "variances" // array should have at least maxdeg + 1 elements. It will contain // the respective variances upon return. This code assumes that // weights and coordinates are statistically independent from // each other. */ template void weightedSampleCoeffVars(const Pair* points, unsigned long nPoints, const double* coeffs, unsigned maxdeg, double* variances) const; /** // Estimate the covariances of the coefficients returned by the // "weightedSampleCoeffs" function. The "coeffs" array should // be at least maxdeg + 1 elements long and should be filled // by a previous call to "weightedSampleCoeffs". The returned // matrix will be dimensioned (maxdeg + 1) x (maxdeg + 1). This // code assumes that weights and coordinates are statistically // independent from each other. */ template npstat::Matrix weightedSampleCoeffCovariance(const Pair* points, unsigned long nPoints, const double* coeffs, unsigned maxdeg) const; /** // Build the coefficients of the orthogonal polynomial series for // the given sample of weighted points (empirical density function). // The first element of the pair will be treated as the coordinate // and the second element will be treated as weight. Weights must // not be negative. 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). Before calculating the coefficients, the // coordinates are shifted and scaled according to // x_new = (x_original - location)/scale. The resulting coefficients // are also divided by scale. */ template void weightedSampleCoeffs(const Pair* points, unsigned long numPoints, typename Pair::first_type location, typename Pair::first_type scale, double* coeffs, unsigned maxdeg) const; /** // Estimate the variances of the coefficients returned by the // "weightedSampleCoeffs" function. The "coeffs" array should // be at least maxdeg + 1 elements long and should be filled // by a previous call to "weightedSampleCoeffs" with the same // location and scale. The "variances" array should have at least // maxdeg + 1 elements. It will contain the respective variances // upon return. This code assumes that weights and coordinates // are statistically independent from each other. */ template void weightedSampleCoeffVars(const Pair* points, unsigned long nPoints, typename Pair::first_type location, typename Pair::first_type scale, const double* coeffs, unsigned maxdeg, double* variances) const; /** // Estimate the covariances of the coefficients returned by the // "weightedSampleCoeffs" function. The "coeffs" array should // be at least maxdeg + 1 elements long and should be filled // by a previous call to "weightedSampleCoeffs" with the same // location and scale. The returned matrix will be dimensioned // (maxdeg + 1) x (maxdeg + 1). This code assumes that weights // and coordinates are statistically independent from each other. */ template npstat::Matrix weightedSampleCoeffCovariance(const Pair* points, unsigned long nPoints, typename Pair::first_type location, typename Pair::first_type scale, 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. */ template double empiricalKroneckerDelta(const Quadrature& quad, unsigned deg1, unsigned deg2) const; /** // A measure-weighted average of a product of four orthonormal // poly values */ template double jointAverage(const Quadrature& q, unsigned deg1, unsigned deg2, unsigned deg3, unsigned deg4) const; /** // A measure-weighted average of a product of multiple orthonormal // poly values. "checkedForZeros" argument should be set to "true" // if we are sure that there are no zero elements in the "degrees" // array. */ template double jointAverage(const Quadrature& q, const unsigned* degrees, unsigned nDegrees, bool checkedForZeros = false) const; /** // Direct unweighted integration of of a product of four orthonormal // poly values. The integration function should support the "integrate" // method without limits. Intended for use with GaussHermiteQuadrature // and such. */ template double directQuadrature(const Quadrature& q, unsigned deg1, unsigned deg2, unsigned deg3, unsigned deg4) const; /** // Direct unweighted integration of of a product of multiple orthonormal // poly values. The integration function should support the "integrate" // method without limits. Intended for use with GaussHermiteQuadrature // and such. "checkedForZeros" argument should be set to "true" if we are // sure that there are no zero elements in the "degrees" array. */ template double directQuadrature(const Quadrature& q, const unsigned* degrees, unsigned nDegrees, bool checkedForZeros = false) const; /** // An unweighted integral of the orthonormal polynomial with the // given degree to some power. For this method, it is assumed // that the polynomials are supported on a closed interval (without // such an assumption unweighted integrals do not make much sense) // and that Gauss-Legendre quadratures can be used. */ double integratePoly(unsigned degree, unsigned power) const; /** // An unweighted integral of a product of multiple orthonormal // polynomials. For this method, it is assumed that the polynomials // are supported on a closed interval (without such an assumption // unweighted integrals do not make much sense) and that Gauss-Legendre // quadratures can be used. */ double jointIntegral(const unsigned* degrees, unsigned nDegrees) const; /** // Unweighted averages of the polynomial values for the given sample. // The length of array "averages" should be at least maxdeg + 1. */ template void sampleAverages(const Numeric* coords, unsigned long lenCoords, double* averages, unsigned maxdeg) const; /** // Averages of the polynomial values for the given sample of // weighted points (not weighting by the polynomial weight function). // The first element of the pair will be treated as the coordinate // and the second element will be treated as weight. Weights must not // be negative. The length of array "averages" should be at least // maxdeg + 1. */ template void weightedPointsAverages(const Pair* points, unsigned long numPoints, double* averages, unsigned maxdeg) const; /** // Averages of the polynomial values weighted by an external // weight function. The length of array "averages" should be // at least maxdeg + 1. */ template void extWeightAverages(const Functor& extWeight, const AbsIntervalQuadrature1D& quad, double* averages, unsigned maxdeg) const; /** // Unweighted averages of the pairwise products of the polynomial // values for the given sample. The returned matrix will be symmetric // and will have the dimensions (maxdeg + 1) x (maxdeg + 1). */ template Matrix sampleProductAverages(const Numeric* coords, unsigned long lenCoords, unsigned maxdeg) const; /** // Pairwise products of the polynomial values weighted by // an external weight function. The returned matrix will be // symmetric and will have the dimensions (maxdeg + 1) x (maxdeg + 1). */ template Matrix extWeightProductAverages(const Functor& extWeight, const AbsIntervalQuadrature1D& quad, unsigned maxdeg) const; /** // Unweighted average of a product of polynomial values for the // given sample. "degrees" is the collection of polynomial degrees. // Polynomials of these degrees will be included in the product. */ template double jointSampleAverage(const Numeric* coords, unsigned long lenCoords, const unsigned* degrees, unsigned lenDegrees) const; // Similar to the previous method but calculates two averages // simultaneously template std::pair twoJointSampleAverages( const Numeric* coords, unsigned long lenCoords, const unsigned* degrees1, unsigned lenDegrees1, const unsigned* degrees2, unsigned lenDegrees2) const; CPP11_auto_ptr makeStorablePolySeries( unsigned maxPolyDeg, const double *coeffs=0, unsigned maxdeg=0, double shift=0.0) const; // Functions needed for compatibility with certain template codes inline long double location() const {return 0.0;} inline long double scale() const {return 1.0;} protected: inline static int kdelta(const unsigned i, const unsigned j) {return i == j ? 1 : 0;} // Recurrence relationship function should return alpha // as the first element of the pair and the square root // of beta as the second virtual std::pair recurrenceCoeffs(unsigned deg) const = 0; private: // Helper classes class ProdFcn : public Functor1 { public: inline ProdFcn(const AbsClassicalOrthoPoly1D& p, const unsigned d1, const unsigned d2) : poly(p), deg1(d1), deg2(d2) {} inline virtual ~ProdFcn() {} inline long double operator()(const long double& x) const { const std::pair& p = poly.twopoly(deg1, deg2, x); return p.first*p.second*poly.weight(x); } private: const AbsClassicalOrthoPoly1D& poly; unsigned deg1; unsigned deg2; }; class MultiProdFcn; friend class MultiProdFcn; class MultiProdFcn : public Functor1 { public: inline MultiProdFcn(const AbsClassicalOrthoPoly1D& p, const unsigned* degrees, unsigned lenDegrees, const bool includeWeight=true) : poly(p), degs(degrees, degrees+lenDegrees), maxdeg(0), weighted(includeWeight) { if (lenDegrees) maxdeg = *std::max_element(degrees, degrees+lenDegrees); } inline virtual ~MultiProdFcn() {} inline long double operator()(const long double& x) const { long double w = 1.0L, p = 1.0L; if (weighted) w = poly.weight(x); if (maxdeg) p = poly.normpolyprod(°s[0], degs.size(), maxdeg, x); return w*p; } private: const AbsClassicalOrthoPoly1D& poly; std::vector degs; unsigned maxdeg; bool weighted; }; // For calling the function below, maxdeg should be calculated as // *std::max_element(degrees, degrees+nDegrees) long double normpolyprod(const unsigned* degrees, unsigned nDegrees, unsigned maxdeg, long double x) const; #ifdef SWIG public: inline double weight2(const double x) const {return weight(x);} inline double poly2(const unsigned deg, const double x) const {return poly(deg, x);} inline std::pair twopoly2(const unsigned deg1, const unsigned deg2, const double x) const { const std::pair& p = twopoly(deg1, deg2, x); return std::pair(p.first, p.second); } inline std::vector allpoly2(const unsigned maxdeg, const double x) const { std::vector p(maxdeg+1); allpoly(x, &p[0], p.size()); return std::vector(p.begin(), p.end()); } inline StorablePolySeries1D* makeStorablePolySeries_2( unsigned maxPolyDeg, const std::vector& coeffs, double shift=0.0) const { CPP11_auto_ptr ptr; if (coeffs.empty()) ptr = this->makeStorablePolySeries(maxPolyDeg, 0, 0, shift); else ptr = this->makeStorablePolySeries(maxPolyDeg, &coeffs[0], coeffs.size() - 1, shift); return ptr.release(); } inline Matrix jacobiMatrix_2(const unsigned n) const {return jacobiMatrix(n);} inline double location2() const {return location();} inline double scale2() const {return scale();} #endif // SWIG }; /** // A functor for the weight function of the given ortho poly system. // The poly system is not copied, only a reference is used. It is // a responsibility of the user to make sure that the lifetime of // the poly system object exceeds the lifetime of the functor. */ class OrthoPoly1DWeight : public Functor1 { public: inline OrthoPoly1DWeight(const AbsClassicalOrthoPoly1D& fcn, const long double normfactor) : fcn_(fcn), norm_(normfactor) {} // Separate constructor here because swig gets confused // by long doubles with defaults inline explicit OrthoPoly1DWeight(const AbsClassicalOrthoPoly1D& fcn) : fcn_(fcn), norm_(1.0L) {} inline virtual ~OrthoPoly1DWeight() {} inline virtual long double operator()(const long double& a) const {return norm_*fcn_.weight(a);} private: OrthoPoly1DWeight(); const AbsClassicalOrthoPoly1D& fcn_; long double norm_; }; /** // A functor for a particular degree polynomial of the given ortho // poly system. The poly system is not copied, only a reference is used. // It is a responsibility of the user to make sure that the lifetime of // the poly system object exceeds the lifetime of the functor. */ class OrthoPoly1DDeg : public Functor1 { public: inline OrthoPoly1DDeg(const AbsClassicalOrthoPoly1D& fcn, const unsigned degree, const long double normfactor) : fcn_(fcn), norm_(normfactor), deg_(degree) { if (deg_ > fcn_.maxDegree()) throw std::invalid_argument( "In npstat::OrthoPoly1DDeg::constructor: " "degree argument is out of range"); } // Separate constructor here because swig gets confused // by long doubles with defaults inline OrthoPoly1DDeg(const AbsClassicalOrthoPoly1D& fcn, const unsigned degree) : fcn_(fcn), norm_(1.0L), deg_(degree) { if (deg_ > fcn_.maxDegree()) throw std::invalid_argument( "In npstat::OrthoPoly1DDeg::constructor: " "degree argument is out of range"); } inline virtual ~OrthoPoly1DDeg() {} inline virtual long double operator()(const long double& a) const {return norm_*fcn_.poly(deg_, a);} private: OrthoPoly1DDeg(); const AbsClassicalOrthoPoly1D& fcn_; long double norm_; unsigned deg_; }; class UnweightedTwoPolyProd : public Functor1 { public: inline UnweightedTwoPolyProd(const AbsClassicalOrthoPoly1D& fcn, const unsigned deg1, const unsigned deg2) : fcn_(fcn), degree1_(deg1), degree2_(deg2) {} inline virtual ~UnweightedTwoPolyProd() {} inline virtual long double operator()(const long double& a) const { std::pair p = fcn_.twopoly( degree1_, degree2_, a); return p.first * p.second; } private: const AbsClassicalOrthoPoly1D& fcn_; unsigned degree1_; unsigned degree2_; }; } #include "npstat/nm/AbsClassicalOrthoPoly1D.icc" #endif // NPSTAT_ABSCLASSICALORTHOPOLY1D_HH_ Index: trunk/npstat/nm/ClassicalOrthoPolys1D.hh =================================================================== --- trunk/npstat/nm/ClassicalOrthoPolys1D.hh (revision 885) +++ trunk/npstat/nm/ClassicalOrthoPolys1D.hh (revision 886) @@ -1,201 +1,239 @@ #ifndef NPSTAT_CLASSICALORTHOPOLYS1D_HH_ #define NPSTAT_CLASSICALORTHOPOLYS1D_HH_ /*! // \file ClassicalOrthoPolys1D.hh // // \brief Orthonormal versions of some classical orthogonal polynomials // // Author: I. Volobouev // // July 2018 */ #include #include #include "npstat/nm/AbsClassicalOrthoPoly1D.hh" namespace npstat { /** Orthonormal Legendre polynomials on [-1, 1] */ class LegendreOrthoPoly1D : public AbsClassicalOrthoPoly1D { public: inline virtual LegendreOrthoPoly1D* clone() const {return new LegendreOrthoPoly1D(*this);} inline virtual ~LegendreOrthoPoly1D() {} inline long double weight(const long double x) const {return (x >= -1.0L && x <= 1.0L) ? 0.5L : 0.0L;} inline double xmin() const {return -1.0L;} inline double xmax() const {return 1.0L;} + /** + // Derive the series coefficients for the integral of the + // argument series from -1 to x. The array of coefficients + // must be at least degree+1 long, and the buffer for the + // integral coefficients must be at least degree+2 long. + */ + void integralCoeffs(const double* coeffs, unsigned degree, + double* integralCoeffs); protected: inline std::pair recurrenceCoeffs(const unsigned k) const { long double sqb = 1.0L; if (k) sqb = 1.0L/sqrtl(4.0L - 1.0L/k/k); return std::pair(0.0L, sqb); } + +#ifdef SWIG + public: + inline std::vector + integralCoeffs2(const double *c, const unsigned degreep1) const + { + std::vector tmp(degreep1+1U); + integralCoeffs(c, degreep1-1U, &tmp[0]); + return tmp; + } +#endif }; /** Orthonormal Legendre polynomials on [0, 1] (that is, shifted) */ class ShiftedLegendreOrthoPoly1D : public AbsClassicalOrthoPoly1D { public: inline virtual ShiftedLegendreOrthoPoly1D* clone() const {return new ShiftedLegendreOrthoPoly1D(*this);} inline virtual ~ShiftedLegendreOrthoPoly1D() {} inline long double weight(const long double x) const {return (x >= 0.0L && x <= 1.0L) ? 1.0L : 0.0L;} inline double xmin() const {return 0.0L;} inline double xmax() const {return 1.0L;} + /** + // Derive the series coefficients for the integral of the + // argument series from 0 to x. The array of coefficients + // must be at least degree+1 long, and the buffer for the + // integral coefficients must be at least degree+2 long. + */ + void integralCoeffs(const double* coeffs, unsigned degree, + double* integralCoeffs); protected: inline std::pair recurrenceCoeffs(const unsigned k) const { long double sqb = 1.0L; if (k) sqb = 0.5L/sqrtl(4.0L - 1.0L/k/k); return std::pair(0.5L, sqb); } + +#ifdef SWIG + public: + inline std::vector + integralCoeffs2(const double *c, const unsigned degreep1) const + { + std::vector tmp(degreep1+1U); + integralCoeffs(c, degreep1-1U, &tmp[0]); + return tmp; + } +#endif }; /** // The weight function for Jacobi polynomials is proportional // to (1 - x)^alpha (1 + x)^beta */ class JacobiOrthoPoly1D : public AbsClassicalOrthoPoly1D { public: JacobiOrthoPoly1D(double alpha, double beta); inline virtual JacobiOrthoPoly1D* clone() const {return new JacobiOrthoPoly1D(*this);} inline virtual ~JacobiOrthoPoly1D() {} inline double alpha() const {return alpha_;} inline double beta() const {return beta_;} long double weight(long double x) const; inline double xmin() const {return -1.0L;} inline double xmax() const {return 1.0L;} /** // The following function generates coefficients for derivatives // of the Jacobi polynomial series. The derivative series should be // evaluated using a JacobiOrthoPoly1D(alpha_+1, beta_+1) object and // with maxdeg argument decreased by 1. The array "coeffs" should // have at least maxdeg+1 elements, while the "derivativeCoeffs" // buffer should have at least maxdeg elements. */ void derivativeCoeffs(const double* coeffs, unsigned maxdeg, double* derivativeCoeffs) const; protected: std::pair recurrenceCoeffs(const unsigned k) const; private: double alpha_; double beta_; long double norm_; #ifdef SWIG public: inline std::vector derivativeCoeffs2(const double *c, const unsigned degreep1) const { std::vector tmp(degreep1-1U); derivativeCoeffs(c, degreep1-1U, tmp.empty() ? (double *)0 : &tmp[0]); return tmp; } #endif }; /** Orthonormal Chebyshev polynomials of the 1st kind */ class ChebyshevOrthoPoly1st : public AbsClassicalOrthoPoly1D { public: inline virtual ChebyshevOrthoPoly1st* clone() const {return new ChebyshevOrthoPoly1st(*this);} inline virtual ~ChebyshevOrthoPoly1st() {} inline long double weight(const long double x) const { const long double Pi = 3.14159265358979323846264338L; return (x > -1.0L && x < 1.0L) ? 1.0/sqrtl(1.0L - x*x)/Pi : 0.0L; } inline double xmin() const {return -1.0L;} inline double xmax() const {return 1.0L;} protected: inline std::pair recurrenceCoeffs(const unsigned k) const { long double sqb = 0.5L; if (k == 0U) sqb = 1.772453850905516027298167483L; // sqrt(Pi) else if (k == 1U) sqb = 0.7071067811865475244008443621L; // 1/sqrt(2) return std::pair(0.0L, sqb); } }; /** Orthonormal Chebyshev polynomials of the 2nd kind */ class ChebyshevOrthoPoly2nd : public AbsClassicalOrthoPoly1D { public: inline virtual ChebyshevOrthoPoly2nd* clone() const {return new ChebyshevOrthoPoly2nd(*this);} inline virtual ~ChebyshevOrthoPoly2nd() {} inline long double weight(const long double x) const { const long double PiOver2 = 1.57079632679489661923132169L; return (x > -1.0L && x < 1.0L) ? sqrtl(1.0L - x*x)/PiOver2 : 0.0L; } inline double xmin() const {return -1.0L;} inline double xmax() const {return 1.0L;} protected: inline std::pair recurrenceCoeffs(const unsigned k) const { long double sqb = 0.5L; if (!k) sqb = 1.25331413731550025120788264L; // sqrt(Pi/2) return std::pair(0.0L, sqb); } }; - /** Orthonormal(!) probabilists' Hermite polynomials */ + /** Orthonormal(!) probabilist's Hermite polynomials */ class HermiteProbOrthoPoly1D : public AbsClassicalOrthoPoly1D { public: inline virtual HermiteProbOrthoPoly1D* clone() const {return new HermiteProbOrthoPoly1D(*this);} inline virtual ~HermiteProbOrthoPoly1D() {} long double weight(long double x) const; inline double xmin() const {return -xmax();} inline double xmax() const {return 150.0;} protected: std::pair recurrenceCoeffs(const unsigned k) const; }; } #endif // NPSTAT_CLASSICALORTHOPOLYS1D_HH_