Index: trunk/tests/test_ClassicalOrthoPolys1D.cc
===================================================================
--- trunk/tests/test_ClassicalOrthoPolys1D.cc	(revision 728)
+++ trunk/tests/test_ClassicalOrthoPolys1D.cc	(revision 729)
@@ -1,422 +1,433 @@
 #include <cmath>
 #include <utility>
 #include <vector>
 #include <sstream>
 
 #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/rng/MersenneTwister.hh"
 
 #include "npstat/stat/Distributions1D.hh"
 #include "npstat/stat/DensityOrthoPoly1D.hh"
 
 using namespace npstat;
 using namespace std;
 
 namespace {
     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; ideg<sizeof(degs)/sizeof(degs[0]); ++ideg)
                 CHECK_EQUAL(poly1->poly(degs[ideg], x), He.poly(degs[ideg], x));
         }
         for (unsigned i=0; i<ntries; ++i)
         {
             const double x = test_rng()*10 - 5;
             for (unsigned ideg=0; ideg<sizeof(degs)/sizeof(degs[0]); ++ideg)
             {
                 CHECK_EQUAL(poly1->poly(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;
 
         LegendreOrthoPoly1D leg;
         CPP11_auto_ptr<StorablePolySeries1D> 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(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<OrthoPoly1DWeight> wpoly_00(
             kernel, densMaxDeg, 2*densMaxDeg+5+4, -1, 1, OPOLY_STIELTJES);
         const ClassicalOrthoPoly1DFromWeight<OrthoPoly1DWeight> 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);
+        const Matrix<double>& m00 = dpoly_00.altWeightProductAverages(
+            DensityFunctor1D(b1), ptest, densMaxDeg);
+        CHECK_EQUAL(densMaxDeg+1U, m00.nRows());
+        CHECK_EQUAL(densMaxDeg+1U, m00.nColumns());
+        Matrix<double> diag(densMaxDeg+1U, densMaxDeg+1U, 1);
+        CHECK((m00 - diag).maxAbsValue() < eps);
+        const Matrix<double>& m01 = dpoly_01.altWeightProductAverages(
+            DensityFunctor1D(b1), ptest, densMaxDeg);
+        CHECK((m01 - diag).maxAbsValue() < 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<OrthoPoly1DWeight> wpoly_20(
             kernel2, densMaxDeg, 2*densMaxDeg+6+3, -1, 1, OPOLY_STIELTJES);
         const ClassicalOrthoPoly1DFromWeight<OrthoPoly1DWeight> 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<double> points(npoints);
         for (unsigned itry=0; itry<ntries; ++itry)
         {
             rng.run(&points[0], npoints, npoints);
             for (unsigned ipt=0; ipt<npoints; ++ipt)
             {
                 points[ipt] *= 2.0;
                 points[ipt] -= 1.0;
             }
             const Matrix<double> avemat(jac.sampleProductAverages(&points[0], npoints, maxdeg));
             std::vector<double> 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<double, double>& 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<double,double> 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<double> points(perPseudo);
         std::vector<Point> wPoints(perPseudo);
         std::vector<double> coeffs(maxdeg + 1U, 100000.0);
         std::vector<double> vars(maxdeg + 1U, -100000.0);
 
         for (unsigned ips=0; ips<nPseudo; ++ips)
         {
             for (unsigned i=0; i<perPseudo; ++i)
             {
                 const double x = test_rng()*scale + loc;
                 const double w = test_rng() + 0.1;
                 points[i] = x;
                 wPoints[i] = Point(x, w);
             }
 
             leg.sampleCoeffs(&points[0], perPseudo, loc, scale,
                              &coeffs[0], maxdeg);
             leg.sampleCoeffVars(&points[0], perPseudo, loc, scale,
                                 &coeffs[0], maxdeg, &vars[0]);
             for (unsigned i=1; i<=maxdeg; ++i)
                 CHECK(vars[i] > 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/nm/AbsClassicalOrthoPoly1D.icc
===================================================================
--- trunk/npstat/nm/AbsClassicalOrthoPoly1D.icc	(revision 728)
+++ trunk/npstat/nm/AbsClassicalOrthoPoly1D.icc	(revision 729)
@@ -1,583 +1,624 @@
 #include <cassert>
 #include <vector>
 #include <stdexcept>
 #include <algorithm>
 
 namespace npstat {
     namespace Private {
         template <class SomePoly1D, class Functor>
         class ClassOrthoPolyHlp : public Functor1<long double, long double>
         {
         public:
             inline ClassOrthoPolyHlp(const SomePoly1D& p,
                                      const Functor& f, const unsigned d1)
                 : poly(p), fcn(f), deg(d1) {}
 
             inline long double operator()(const long double& x) const
                 {return fcn(x)*poly.poly(deg, x)*poly.weight(x);}
 
         private:
             const SomePoly1D& poly;
             const Functor& fcn;
             unsigned deg;
         };
     }
 
     template <class Functor, class Quadrature>
     void AbsClassicalOrthoPoly1D::calculateCoeffs(
         const Functor& fcn, const Quadrature& quad,
         double *coeffs, const unsigned maxdeg) const
     {
         assert(coeffs);
 
         const double a = xmin();
         const double b = xmax();
 
         for (unsigned i=0; i<=maxdeg; ++i)
         {
             Private::ClassOrthoPolyHlp<AbsClassicalOrthoPoly1D,Functor>
                 func(*this, fcn, i);
             coeffs[i] = quad.integrate(func, a, b);
         }
     }
 
     template <class Quadrature>
     double AbsClassicalOrthoPoly1D::empiricalKroneckerDelta(
         const Quadrature& quad, const unsigned deg1, const unsigned deg2) const
     {
         ProdFcn fcn(*this, deg1, deg2);
         return quad.integrate(fcn, xmin(), xmax());
     }
 
     template <class Quadrature>
     double AbsClassicalOrthoPoly1D::jointAverage(
         const Quadrature& quad, const unsigned deg1, const unsigned deg2,
         const unsigned deg3, const unsigned deg4) const
     {
         unsigned degs[4];
         unsigned nDegs = 0;
         if (deg1)
             degs[nDegs++] = deg1;
         if (deg2)
             degs[nDegs++] = deg2;
         if (deg3)
             degs[nDegs++] = deg3;
         if (deg4)
             degs[nDegs++] = deg4;
         return jointAverage(quad, degs, nDegs, true);
     }
 
     template <class Quadrature>
     double AbsClassicalOrthoPoly1D::jointAverage(
         const Quadrature& quad, const unsigned* degrees,
         const unsigned nDegrees, const bool checkedForZeros) const
     {
         if (nDegrees)
         {
             assert(degrees);
 
             if (nDegrees == 1U)
                 return kdelta(degrees[0], 0U);
             if (nDegrees == 2U)
                 return kdelta(degrees[0], degrees[1]);
 
             if (!checkedForZeros)
             {
                 bool allNonZero = true;
                 for (unsigned ideg=0; ideg<nDegrees && allNonZero; ++ideg)
                     allNonZero = degrees[ideg];
                 if (!allNonZero)
                 {
                     unsigned degBuf[256];
                     if (nDegrees - 1U <= sizeof(degBuf)/sizeof(degBuf[0]))
                     {
                         unsigned nonZeroDegs = 0;
                         for (unsigned ideg=0; ideg<nDegrees; ++ideg)
                             if (degrees[ideg])
                                 degBuf[nonZeroDegs++] = degrees[ideg];
                         return jointAverage(quad, degBuf, nonZeroDegs, true);
                     }
                 }
             }
 
             MultiProdFcn fcn(*this, degrees, nDegrees);
             return quad.integrate(fcn, xmin(), xmax());
         }
         else
             return 1.0;
     }
 
     template <class Numeric>
     Matrix<double> AbsClassicalOrthoPoly1D::sampleProductAverages(
         const Numeric* coords, const unsigned long lenCoords,
         const unsigned maxdeg) const
     {
         if (!lenCoords) throw std::invalid_argument(
             "In npstat::AbsClassicalOrthoPoly1D::sampleProductAverages:"
             " empty sample");
         assert(coords);
         const unsigned long nsums = (maxdeg + 1)*(unsigned long)(maxdeg + 2)/2;
         std::vector<long double> membuf(nsums + (maxdeg + 1U), 0.0L);
         long double* poly = &membuf[0];
         long double* sums = poly + (maxdeg + 1U);
         for (unsigned long i=0; i<lenCoords; ++i)
         {
             allpoly(coords[i], poly, maxdeg);
             unsigned long isum = 0;
             for (unsigned k=0; k<=maxdeg; ++k)
                 for (unsigned j=0; j<=k; ++j, ++isum)
                     sums[isum] += poly[k]*poly[j];
         }
         Matrix<double> mat(maxdeg + 1U, maxdeg + 1U);
         unsigned long isum = 0;
         for (unsigned k=0; k<=maxdeg; ++k)
             for (unsigned j=0; j<=k; ++j, ++isum)
             {
                 const double v = sums[isum]/lenCoords;
                 mat[k][j] = v;
                 if (j != k)
                     mat[j][k] = v;
             }
         return mat;
     }
 
+    template <class Functor>
+    Matrix<double> AbsClassicalOrthoPoly1D::altWeightProductAverages(
+        const Functor& fcn, const AbsIntervalQuadrature1D& quad,
+        const unsigned maxdeg) const
+    {
+        const unsigned nInteg = quad.npoints();
+        const unsigned long nsums = (maxdeg + 1)*(unsigned long)(maxdeg + 2)/2;
+        std::vector<long double> membuf(nsums + (maxdeg + 1U) + 2U*nInteg, 0.0L);
+        long double* poly = &membuf[0];
+        long double* sums = poly + (maxdeg + 1U);
+        long double* absc = sums + nsums;
+        long double* weights = absc + nInteg;
+
+        quad.getAllAbscissae(absc, nInteg);
+        quad.getAllWeights(weights, nInteg);
+        const long double midpoint = (static_cast<long double>(xmin()) + xmax())/2.0L;
+        const long double unit = (static_cast<long double>(xmax()) - xmin())/2.0L;
+
+        for (unsigned long i=0; i<nInteg; ++i)
+        {
+            const long double x = midpoint + absc[i]*unit;
+            const long double w = fcn(x)*weights[i]*unit;
+            allpoly(x, poly, maxdeg);
+            unsigned long isum = 0;
+            for (unsigned k=0; k<=maxdeg; ++k)
+                for (unsigned j=0; j<=k; ++j, ++isum)
+                    sums[isum] += poly[k]*poly[j]*w;
+        }
+        Matrix<double> mat(maxdeg + 1U, maxdeg + 1U);
+        unsigned long isum = 0;
+        for (unsigned k=0; k<=maxdeg; ++k)
+            for (unsigned j=0; j<=k; ++j, ++isum)
+            {
+                const double v = sums[isum];
+                mat[k][j] = v;
+                if (j != k)
+                    mat[j][k] = v;
+            }
+        return mat;
+    }
+
     template <class Numeric>
     void AbsClassicalOrthoPoly1D::sampleAverages(
         const Numeric* coords, const unsigned long lenCoords,
         double *averages, const unsigned maxdeg) const
     {
         if (!lenCoords) throw std::invalid_argument(
             "In npstat::AbsClassicalOrthoPoly1D::sampleAverages: empty sample");
         assert(coords);
         assert(averages);
 
         std::vector<long double> membuf(2U*(maxdeg + 1U), 0.0L);
         long double* poly = &membuf[0];
         long double* sums = poly + (maxdeg + 1U);
 
         for (unsigned long i=0; i<lenCoords; ++i)
         {
             allpoly(coords[i], poly, maxdeg);
             for (unsigned k=0; k<=maxdeg; ++k)
                 sums[k] += poly[k];
         }
         const long double ldsz = lenCoords;
         for (unsigned k=0; k<=maxdeg; ++k)
             averages[k] = sums[k]/ldsz;
     }
 
     template <class Numeric>
     void AbsClassicalOrthoPoly1D::sampleCoeffs(
         const Numeric* coords, const unsigned long lenCoords,
         const Numeric location, const Numeric scale,
         double *coeffs, const unsigned maxdeg) const
     {
         if (!lenCoords) throw std::invalid_argument(
             "In npstat::AbsClassicalOrthoPoly1D::sampleCoeffs: empty sample");
         assert(coords);
         assert(coeffs);
         const Numeric zero = Numeric();
         if (scale == zero) throw std::invalid_argument(
             "In npstat::AbsClassicalOrthoPoly1D::sampleCoeffs: zero scale");
         std::vector<long double> membuf(2U*(maxdeg + 1U), 0.0L);
         long double* poly = &membuf[0];
         long double* sums = poly + (maxdeg + 1U);
         for (unsigned long i=0; i<lenCoords; ++i)
         {
             const long double x = (coords[i] - location)/scale;
             const long double w = weight(x);
             allpoly(x, poly, maxdeg);
             for (unsigned k=0; k<=maxdeg; ++k)
                 sums[k] += poly[k]*w;
         }
         for (unsigned k=0; k<=maxdeg; ++k)
             coeffs[k] = sums[k]/scale/lenCoords;
     }
 
     template <class Numeric>
     void AbsClassicalOrthoPoly1D::sampleCoeffVars(
         const Numeric* coords, const unsigned long lenCoords,
         const Numeric location, const Numeric scale,
         const double *coeffs, const unsigned maxdeg,
         double *variances) const
     {
         if (lenCoords < 2UL) throw std::invalid_argument(
             "In npstat::AbsClassicalOrthoPoly1D::sampleCoeffVars: "
             "insufficient sample size");
         assert(coords);
         assert(coeffs);
         assert(variances);
         const Numeric zero = Numeric();
         if (scale == zero) throw std::invalid_argument(
             "In npstat::AbsClassicalOrthoPoly1D::sampleCoeffVars: zero scale");
         std::vector<long double> membuf(2U*(maxdeg + 1U), 0.0L);
         long double* poly = &membuf[0];
         long double* sums = poly + (maxdeg + 1U);
         for (unsigned long i=0; i<lenCoords; ++i)
         {
             const long double x = (coords[i] - location)/scale;
             const long double w = weight(x)/scale;
             allpoly(x, poly, maxdeg);
             for (unsigned k=0; k<=maxdeg; ++k)
             {
                 const long double delta = poly[k]*w - coeffs[k];
                 sums[k] += delta*delta;
             }
         }
         long double lm1 = lenCoords - 1UL;
         lm1 *= lenCoords;
         for (unsigned k=0; k<=maxdeg; ++k)
             variances[k] = sums[k]/lm1;
     }
 
     template <class Numeric>
     npstat::Matrix<double>
     AbsClassicalOrthoPoly1D::sampleCoeffCovariance(
         const Numeric* coords, const unsigned long lenCoords,
         const Numeric location, const Numeric scale,
         const double *coeffs, const unsigned maxdeg) const
     {
         if (lenCoords < 2UL) throw std::invalid_argument(
             "In npstat::AbsClassicalOrthoPoly1D::sampleCoeffCovariance: "
             "insufficient sample size");
         assert(coords);
         assert(coeffs);
         const Numeric zero = Numeric();
         if (scale == zero) throw std::invalid_argument(
             "In npstat::AbsClassicalOrthoPoly1D::sampleCoeffCovariance: "
             "zero scale");
         std::vector<long double> membuf(maxdeg + 1U, 0.0L);
         long double* poly = &membuf[0];
         npstat::Matrix<long double> sums(maxdeg + 1U, maxdeg + 1U, 0);
         for (unsigned long i=0; i<lenCoords; ++i)
         {
             const long double x = (coords[i] - location)/scale;
             const long double w = weight(x)/scale;
             allpoly(x, poly, maxdeg);
             for (unsigned k=0; k<=maxdeg; ++k)
             {
                 const long double deltak = poly[k]*w - coeffs[k];
                 long double* tmp = sums[k];
                 for (unsigned m=0; m<=k; ++m)
                 {
                     const long double deltam = poly[m]*w - coeffs[m];
                     tmp[m] += deltak*deltam;
                 }
             }
         }
         for (unsigned k=0; k<=maxdeg; ++k)
             for (unsigned m=k+1; m<=maxdeg; ++m)
                 sums[k][m] = sums[m][k];
         long double lm1 = lenCoords - 1UL;
         lm1 *= lenCoords;
         sums /= lm1;
         return sums;
     }
 
     template <class Numeric>
     void AbsClassicalOrthoPoly1D::sampleCoeffs(
         const Numeric* coords, const unsigned long lenCoords,
         double *coeffs, const unsigned maxdeg) const
     {
         const Numeric zero = Numeric();
         const Numeric one = static_cast<Numeric>(1);
         sampleCoeffs(coords, lenCoords, zero, one, coeffs, maxdeg);
     }
 
     template <class Numeric>
     void AbsClassicalOrthoPoly1D::sampleCoeffVars(
         const Numeric* coords, const unsigned long lenCoords,
         const double *coeffs, const unsigned maxdeg,
         double *variances) const
     {
         const Numeric zero = Numeric();
         const Numeric one = static_cast<Numeric>(1);
         sampleCoeffVars(coords, lenCoords, zero, one,
                         coeffs, maxdeg, variances);
     }
 
     template <class Numeric>
     npstat::Matrix<double>
     AbsClassicalOrthoPoly1D::sampleCoeffCovariance(
         const Numeric* coords, unsigned long lenCoords,
         const double* coeffs, unsigned maxdeg) const
     {
         const Numeric zero = Numeric();
         const Numeric one = static_cast<Numeric>(1);
         return sampleCoeffCovariance(coords, lenCoords, zero,
                                      one, coeffs, maxdeg);
     }
 
     template <class Pair>
     void AbsClassicalOrthoPoly1D::weightedSampleCoeffs(
         const Pair* points, const unsigned long numPoints,
         const typename Pair::first_type location,
         const typename Pair::first_type scale,
         double *coeffs, const unsigned maxdeg) const
     {
         typedef typename Pair::first_type Numeric;
 
         if (!numPoints) throw std::invalid_argument(
             "In npstat::AbsClassicalOrthoPoly1D::weightedSampleCoeffs: "
             "empty sample");
         assert(points);
         assert(coeffs);
         const Numeric zero = Numeric();
         if (scale == zero) throw std::invalid_argument(
             "In npstat::AbsClassicalOrthoPoly1D::weightedSampleCoeffs: "
             "zero scale");
         std::vector<long double> membuf(2U*(maxdeg + 1U), 0.0L);
         long double* poly = &membuf[0];
         long double* sums = poly + (maxdeg + 1U);
         long double wsum = 0.0L;
         for (unsigned long i=0; i<numPoints; ++i)
         {
             const long double ptw = points[i].second;
             if (ptw < 0.0L) throw std::invalid_argument(
                 "In npstat::AbsClassicalOrthoPoly1D::weightedSampleCoeffs: "
                 "all weights must be non-negative");
             wsum += ptw;
             const long double x = (points[i].first - location)/scale;
             const long double w = weight(x);
             allpoly(x, poly, maxdeg);
             for (unsigned k=0; k<=maxdeg; ++k)
                 sums[k] += poly[k]*w*ptw;
         }
         if (wsum == 0.0L) throw std::invalid_argument(
             "In npstat::AbsClassicalOrthoPoly1D::weightedSampleCoeffs: "
             "at least one weight must be positive");
         for (unsigned k=0; k<=maxdeg; ++k)
             coeffs[k] = sums[k]/wsum/scale;
     }
 
     template <class Pair>
     void AbsClassicalOrthoPoly1D::weightedSampleCoeffVars(
         const Pair* points, const unsigned long numPoints,
         typename Pair::first_type location,
         typename Pair::first_type scale,
         const  double *coeffs, const unsigned maxdeg,
         double *variances) const
     {
         typedef typename Pair::first_type Numeric;
 
         if (numPoints < 2UL) throw std::invalid_argument(
             "In npstat::AbsClassicalOrthoPoly1D::weightedSampleCoeffVars: "
             "insufficient sample size");
         assert(points);
         assert(coeffs);
         assert(variances);
         const Numeric zero = Numeric();
         if (scale == zero) throw std::invalid_argument(
             "In npstat::AbsClassicalOrthoPoly1D::weightedSampleCoeffVars: "
             "zero scale");
         std::vector<long double> membuf(2U*(maxdeg + 1U), 0.0L);
         long double* poly = &membuf[0];
         long double* sums = poly + (maxdeg + 1U);
         long double wsum = 0.0L, wsqsum = 0.0L;
         for (unsigned long i=0; i<numPoints; ++i)
         {
             const long double ptw = points[i].second;
             if (ptw < 0.0L) throw std::invalid_argument(
                 "In npstat::AbsClassicalOrthoPoly1D::weightedSampleCoeffVars: "
                 "all weights must be non-negative");
             wsum += ptw;
             wsqsum += ptw*ptw;
             const long double x = (points[i].first - location)/scale;
             const long double w = weight(x)/scale;
             allpoly(x, poly, maxdeg);
             for (unsigned k=0; k<=maxdeg; ++k)
             {
                 const long double delta = poly[k]*w - coeffs[k];
                 sums[k] += delta*delta*ptw;
             }
         }
         if (wsum == 0.0L) throw std::invalid_argument(
             "In npstat::AbsClassicalOrthoPoly1D::weightedSampleCoeffVars: "
             "at least one weight must be positive");
         long double lm1 = wsum*wsum/wsqsum - 1.0L;
         if (lm1 <= 0.0L) throw std::invalid_argument(
             "In npstat::AbsClassicalOrthoPoly1D::weightedSampleCoeffVars: "
             "insufficient effective sample size");
         lm1 *= wsum;
         for (unsigned k=0; k<=maxdeg; ++k)
             variances[k] = sums[k]/lm1;
     }
 
     template <class Pair>
     npstat::Matrix<double>
     AbsClassicalOrthoPoly1D::weightedSampleCoeffCovariance(
         const Pair* points, const unsigned long numPoints,
         typename Pair::first_type location,
         typename Pair::first_type scale,
         const double *coeffs, const unsigned maxdeg) const
     {
         typedef typename Pair::first_type Numeric;
 
         if (numPoints < 2UL) throw std::invalid_argument(
             "In npstat::AbsClassicalOrthoPoly1D::weightedSampleCoeffCovariance: "
             "insufficient sample size");
         assert(points);
         assert(coeffs);
         const Numeric zero = Numeric();
         if (scale == zero) throw std::invalid_argument(
             "In npstat::AbsClassicalOrthoPoly1D::weightedSampleCoeffCovariance: "
             "zero scale");
         std::vector<long double> membuf(maxdeg + 1U, 0.0L);
         long double* poly = &membuf[0];
         npstat::Matrix<long double> sums(maxdeg + 1U, maxdeg + 1U, 0);
         long double wsum = 0.0L, wsqsum = 0.0L;
         for (unsigned long i=0; i<numPoints; ++i)
         {
             const long double ptw = points[i].second;
             if (ptw < 0.0L) throw std::invalid_argument(
                 "In npstat::AbsClassicalOrthoPoly1D::weightedSampleCoeffCovariance: "
                 "all weights must be non-negative");
             wsum += ptw;
             wsqsum += ptw*ptw;
             const long double x = (points[i].first - location)/scale;
             const long double w = weight(x)/scale;
             allpoly(x, poly, maxdeg);
             for (unsigned k=0; k<=maxdeg; ++k)
             {
                 const long double deltak = poly[k]*w - coeffs[k];
                 long double* tmp = sums[k];
                 for (unsigned m=0; m<=k; ++m)
                 {
                     const long double deltam = poly[m]*w - coeffs[m];
                     tmp[m] += deltak*deltam*ptw;;
                 }
             }
         }
         if (wsum == 0.0L) throw std::invalid_argument(
             "In npstat::AbsClassicalOrthoPoly1D::weightedSampleCoeffCovariance: "
             "at least one weight must be positive");
         long double lm1 = wsum*wsum/wsqsum - 1.0L;
         if (lm1 <= 0.0L) throw std::invalid_argument(
             "In npstat::AbsClassicalOrthoPoly1D::weightedSampleCoeffCovariance: "
             "insufficient effective sample size");
         lm1 *= wsum;
         for (unsigned k=0; k<=maxdeg; ++k)
             for (unsigned m=k+1; m<=maxdeg; ++m)
                 sums[k][m] = sums[m][k];
         sums /= lm1;
         return sums;
     }
 
     template <class Pair>
     void AbsClassicalOrthoPoly1D::weightedSampleCoeffs(
         const Pair* points, const unsigned long numPoints,
         double *coeffs, const unsigned maxdeg) const
     {
         typedef typename Pair::first_type Numeric;
 
         const Numeric zero = Numeric();
         const Numeric one = static_cast<Numeric>(1);
         weightedSampleCoeffs(points, numPoints, zero, one, coeffs, maxdeg);
     }
 
     template <class Pair>
     void AbsClassicalOrthoPoly1D::weightedSampleCoeffVars(
         const Pair* points, const unsigned long nPoints,
         const  double *coeffs, const unsigned maxdeg,
         double *variances) const
     {
         typedef typename Pair::first_type Numeric;
 
         const Numeric zero = Numeric();
         const Numeric one = static_cast<Numeric>(1);
         weightedSampleCoeffVars(points, nPoints, zero, one,
                                 coeffs, maxdeg, variances);
     }
 
     template <class Pair>
     npstat::Matrix<double>
     AbsClassicalOrthoPoly1D::weightedSampleCoeffCovariance(
         const Pair* points, unsigned long nPoints,
         const double* coeffs, unsigned maxdeg) const
     {
         typedef typename Pair::first_type Numeric;
 
         const Numeric zero = Numeric();
         const Numeric one = static_cast<Numeric>(1);
         return weightedSampleCoeffCovariance(points, nPoints, zero,
                                              one, coeffs, maxdeg);
     }
 
     template <class Numeric>
     double AbsClassicalOrthoPoly1D::jointSampleAverage(
         const Numeric* coords, const unsigned long lenCoords,
         const unsigned* degrees, const unsigned lenDegrees) const
     {
         if (!lenCoords) throw std::invalid_argument(
             "In npstat::AbsClassicalOrthoPoly1D::jointSampleAverage: "
             "empty sample");
         assert(coords);
         if (lenDegrees)
         {
             assert(degrees);
             const unsigned maxdeg = *std::max_element(degrees, degrees+lenDegrees);
             std::vector<long double> membuf(maxdeg + 1U, 0.0L);
             long double* poly = &membuf[0];
             long double sum = 0.0L;
             for (unsigned long i=0; i<lenCoords; ++i)
             {
                 allpoly(coords[i], poly, maxdeg);
                 long double prod = 1.0L;
                 for (unsigned k=0; k<lenDegrees; ++k)
                     prod *= poly[degrees[k]];
                 sum += prod;
             }
             return sum/lenCoords;
         }
         else
             return 1.0;
     }
 
     template <class Numeric>
     std::pair<double, double> AbsClassicalOrthoPoly1D::twoJointSampleAverages(
         const Numeric* coords, const unsigned long lenCoords,
         const unsigned* degrees1, const unsigned lenDegrees1,
         const unsigned* degrees2, const unsigned lenDegrees2) const
     {
         if (!lenCoords) throw std::invalid_argument(
             "In npstat::AbsClassicalOrthoPoly1D::twoJointSampleAverages: empty sample");
         assert(coords);
 
         unsigned maxdeg1 = 0, maxdeg2 = 0;
         if (lenDegrees1)
         {
             assert(degrees1);
             maxdeg1 = *std::max_element(degrees1, degrees1+lenDegrees1);
         }
         if (lenDegrees2)
         {
             assert(degrees2);
             maxdeg2 = *std::max_element(degrees2, degrees2+lenDegrees2);
         }
         const unsigned maxdeg = std::max(maxdeg1, maxdeg2);
 
         if (lenDegrees1 || lenDegrees2)
         {
             std::vector<long double> membuf(maxdeg + 1U, 0.0L);
             long double* poly = &membuf[0];
             long double sum1 = 0.0L, sum2 = 0.0L;
             for (unsigned long i=0; i<lenCoords; ++i)
             {
                 allpoly(coords[i], poly, maxdeg);
 
                 long double prod1 = 1.0L;
                 for (unsigned k=0; k<lenDegrees1; ++k)
                     prod1 *= poly[degrees1[k]];
                 sum1 += prod1;
 
                 long double prod2 = 1.0L;
                 for (unsigned k=0; k<lenDegrees2; ++k)
                     prod2 *= poly[degrees2[k]];
                 sum2 += prod2;
             }
             return std::pair<double, double>(sum1/lenCoords, sum2/lenCoords);
         }
         else
             return std::pair<double, double>(1.0, 1.0);
     }
 }
Index: trunk/npstat/nm/AbsClassicalOrthoPoly1D.hh
===================================================================
--- trunk/npstat/nm/AbsClassicalOrthoPoly1D.hh	(revision 728)
+++ trunk/npstat/nm/AbsClassicalOrthoPoly1D.hh	(revision 729)
@@ -1,504 +1,514 @@
 #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 <utility>
 #include <algorithm>
 #include <climits>
 #include <vector>
 #include <stdexcept>
 
 #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 AbsClassicalOrthoPoly1D
     {
     public:
         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;}
 
         /** Polynomial values */
         long double poly(unsigned deg, long double x) const;
 
         /** 
         // Values of two orthonormal polynomials.
         // Faster than calling "poly" two times.
         */
         std::pair<long double,long double> 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 */
         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 <class Functor, class Quadrature>
         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 <class Numeric>
         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 <class Numeric>
         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 <class Numeric>
         npstat::Matrix<double>
         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 <class Numeric>
         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 <class Numeric>
         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 <class Numeric>
         npstat::Matrix<double>
         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 <class Pair>
         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 <class Pair>
         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 <class Pair>
         npstat::Matrix<double>
         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 <class Pair>
         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 <class Pair>
         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 <class Pair>
         npstat::Matrix<double>
         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 <class Quadrature>
         double empiricalKroneckerDelta(const Quadrature& quad,
                                        unsigned deg1, unsigned deg2) const;
 
         /**
         // A measure-weighted average of a product of four orthonormal
         // poly values
         */
         template <class Quadrature>
         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 <class Quadrature>
         double jointAverage(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 <class Numeric>
         void sampleAverages(const Numeric* coords, unsigned long lenCoords,
                             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 <class Numeric>
         Matrix<double> sampleProductAverages(const Numeric* coords,
                                              unsigned long lenCoords,
                                              unsigned maxdeg) const;
 
         /**
+        // Pairwise products of the polynomial values weighted by
+        // a different weight function. The returned matrix will be
+        // symmetric and will have the dimensions (maxdeg + 1) x (maxdeg + 1).
+        */
+        template <class Functor>
+        Matrix<double> altWeightProductAverages(const Functor& altWeight,
+            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 <class Numeric>
         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 <class Numeric>
         std::pair<double, double> twoJointSampleAverages(
             const Numeric* coords, unsigned long lenCoords,
             const unsigned* degrees1, unsigned lenDegrees1,
             const unsigned* degrees2, unsigned lenDegrees2) const;
 
         CPP11_auto_ptr<StorablePolySeries1D> makeStorablePolySeries(
             unsigned maxPolyDeg,
             const double *coeffs=0, unsigned maxdeg=0) const;
 
     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<long double,long double>
         recurrenceCoeffs(unsigned deg) const = 0;
 
     private:
         // Helper classes
         class ProdFcn : public Functor1<long double, long double>
         {
         public:
             inline ProdFcn(const AbsClassicalOrthoPoly1D& p,
                            const unsigned d1,
                            const unsigned d2)
                 : poly(p), deg1(d1), deg2(d2) {}
 
             inline long double operator()(const long double& x) const
             {
                 const std::pair<long double,long double>& 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<long double, long double>
         {
         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 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(&degs[0], degs.size(), maxdeg, x);
                 return w*p;
             }
 
         private:
             const AbsClassicalOrthoPoly1D& poly;
             std::vector<unsigned> 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::vector<double> allpoly2(const unsigned maxdeg, const double x) const
         {
             std::vector<long double> p(maxdeg+1);
             allpoly(x, &p[0], p.size());
             return std::vector<double>(p.begin(), p.end());
         }
 
         inline StorablePolySeries1D* makeStorablePolySeries_2(
             unsigned maxPolyDeg, const std::vector<double>& coeffs) const
         {
             CPP11_auto_ptr<StorablePolySeries1D> ptr;
             if (coeffs.empty())
                 ptr = this->makeStorablePolySeries(maxPolyDeg);
             else
                 ptr = this->makeStorablePolySeries(maxPolyDeg, &coeffs[0],
                                                    coeffs.size() - 1);
             return ptr.release();
         }
 #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<long double, long double>
     {
     public:
         inline explicit OrthoPoly1DWeight(const AbsClassicalOrthoPoly1D& fcn,
                                           const long double normfactor=1.0L)
             : fcn_(fcn), norm_(normfactor) {}
 
         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<long double, long double>
     {
     public:
         inline OrthoPoly1DDeg(const AbsClassicalOrthoPoly1D& fcn,
                               const unsigned degree,
                               const long double normfactor=1.0L)
             : fcn_(fcn), norm_(normfactor), 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_;
     };
 }
 
 #include "npstat/nm/AbsClassicalOrthoPoly1D.icc"
 
 #endif // NPSTAT_ABSCLASSICALORTHOPOLY1D_HH_