Index: trunk/tests/jacobiEpsTest.cc
===================================================================
--- trunk/tests/jacobiEpsTest.cc	(revision 555)
+++ trunk/tests/jacobiEpsTest.cc	(revision 556)
@@ -1,189 +1,185 @@
 #include <string>
 #include <sstream>
 #include <iostream>
 
 #include "test_utils.hh"
 
 #include "npstat/nm/ClassicalOrthoPolys1D.hh"
 #include "npstat/nm/ContOrthoPoly1D.hh"
 #include "npstat/nm/GaussLegendreQuadrature.hh"
 
 #include "npstat/stat/Distributions1D.hh"
 #include "npstat/stat/SampleAccumulator.hh"
 #include "npstat/stat/orthoPoly1DVProducts.hh"
 
 #include "npstat/rng/MersenneTwister.hh"
 
 using namespace npstat;
 using namespace std;
 
 
 inline static double kdelta(const unsigned i,
                             const unsigned j)
 {
     return i == j ? 1.0 : 0.0;
 }
 
 
 static double calculateExpectation(
     const JacobiOrthoPoly1D& jac, const unsigned nPoints,
     const unsigned m, const unsigned n, const bool highOrder)
 {
     const unsigned sumdeg = (unsigned)(jac.alpha()) + (unsigned)(jac.beta()) +
                             4U*std::max(m, n);
     const unsigned nGood = GaussLegendreQuadrature::minimalExactRule(sumdeg);
     assert(nGood);
     GaussLegendreQuadrature glq(nGood);
     return oracleEpsExpectation(jac, glq, nPoints, m, n, highOrder);
 }
 
 
 static double calculateCovariance(
     const JacobiOrthoPoly1D& jac, const unsigned nPoints,
     const unsigned m1, const unsigned n1,
     const unsigned m2, const unsigned n2,
     const bool highOrder)
 {
     const unsigned sumdeg = (unsigned)(jac.alpha()) + (unsigned)(jac.beta()) +
         8U*std::max(std::max(m1, n1), std::max(m2, n2));
     const unsigned nGood = GaussLegendreQuadrature::minimalExactRule(sumdeg);
     assert(nGood);
     GaussLegendreQuadrature glq(nGood);
     return oracleEpsCovariance(jac, glq, nPoints, m1, n1, m2, n2, highOrder);
 }
 
 
 static string eps_name(const unsigned m, const unsigned n, const bool isReal)
 {
     ostringstream os;
     if (isReal)
         os << "r_";
     os << "eps_" << m << '_' << n;
     return os.str();
 }
 
 
 template <class Acc>
 static void print_accumulator(const string& which,
                               const Acc& acc, const double expMean)
 {
     const double m = acc.mean();
     const double u = acc.meanUncertainty();
     const double nsig = (m - expMean)/u;
     
     cout << which << " average is " << m << " +- " << u
          << ", expected " << expMean << " (" << nsig << " sigma)" << endl;
 }
 
 
 int main(int argc, char const* argv[])
 {
     if (argc != 10)
     {
         cout << "Usage: " << parse_progname(argv[0])
              << " alpha beta pointsPerPseudo nPseudo m n p q highOrder" << endl;
         exit(1);
     }
 
     unsigned alpha, beta;
     if (parse_unsigned(argv[1], &alpha)) exit(1);
     if (parse_unsigned(argv[2], &beta)) exit(1);
 
     unsigned perpseudo;
     if (parse_unsigned(argv[3], &perpseudo)) exit(1);
     if (!perpseudo)
     {
         cerr << "pointsPerPseudo must be positive" << endl;
         exit(1);
     }
 
     unsigned npseudo;
     if (parse_unsigned(argv[4], &npseudo)) exit(1);
     if (npseudo < 2)
     {
         cerr << "nPseudo must be at least 2" << endl;
         exit(1);
     }
 
     unsigned m, n, p, q;
     if (parse_unsigned(argv[5], &m)) exit(1);
     if (parse_unsigned(argv[6], &n)) exit(1);
     if (parse_unsigned(argv[7], &p)) exit(1);
     if (parse_unsigned(argv[8], &q)) exit(1);
     const unsigned maxdeg = std::max(std::max(m, n), std::max(p, q));
 
     bool highOrder;
     if (parse_bool(argv[9], &highOrder)) exit(1);
 
     vector<double> rn;
     rn.reserve(perpseudo);
 
     Beta1D distro(-1.0, 2.0, beta+1, alpha+1);
     MersenneTwister rng;
     JacobiOrthoPoly1D jac(alpha, beta);
     OrthoPoly1DWeight jacWeight(jac);
 
     const unsigned sumdeg = (unsigned)(jac.alpha()) + (unsigned)(jac.beta()) + 2U*maxdeg;
     const unsigned integPoints = GaussLegendreQuadrature::minimalExactRule(sumdeg);
 
     SampleAccumulator<double> acc0, acc1, acc2, acc3;
     acc0.reserve(npseudo);
     acc1.reserve(npseudo);
     acc2.reserve(npseudo);
     acc3.reserve(npseudo);
 
     for (unsigned ips=0; ips<npseudo; ++ips)
     {
         rn.clear();
         double tmp;
         for (unsigned i=0; i<perpseudo; ++i)
         {
             distro.random(rng, &tmp);
             rn.push_back(tmp);
         }
         ContOrthoPoly1D samplePoly(maxdeg, rn);
 
         const double v0 = sampleEpsValue(jac, &rn[0], perpseudo, m, n, highOrder);
         acc0.accumulate(v0);
 
         const double v1 = sampleEpsValue(jac, &rn[0], perpseudo, p, q, highOrder);
         acc1.accumulate(v1);
 
-        const double winteg = samplePoly.integrateEWPolyProduct(
-            jacWeight, 0, 1, -1.0, 1.0, integPoints);
-        cout << "winteg is " << winteg << endl;
-
         const double eps0 = samplePoly.integrateEWPolyProduct(
             jacWeight, m, n, -1.0, 1.0, integPoints) - kdelta(m, n);
         acc2.accumulate(eps0);
 
         const double eps1 = samplePoly.integrateEWPolyProduct(
             jacWeight, p, q, -1.0, 1.0, integPoints) - kdelta(m, n);
         acc3.accumulate(eps1);
     }
 
     const std::string& name0 = eps_name(m, n, false);
     const std::string& rname0 = eps_name(m, n, true);
     const double em0 = calculateExpectation(jac, perpseudo, m, n, highOrder);
     print_accumulator(name0, acc0, em0);
     cout << endl;
     print_accumulator(rname0, acc2, em0);
     cout << endl;
 
     const std::string& name1 = eps_name(p, q, false);
     const std::string& rname1 = eps_name(p, q, true);
     const double em1 = calculateExpectation(jac, perpseudo, p, q, highOrder);
     print_accumulator(name1, acc1, em1);
     cout << endl;
     print_accumulator(rname1, acc3, em1);
     cout << endl;
 
     const double cov = calculateCovariance(jac, perpseudo, m, n, p, q, highOrder);
     cout << name0 << '/' << name1 << " covariance is " << acc0.cov(acc1)
          << ", expected " << cov << endl;
     cout << endl;
     cout << rname0 << '/' << rname1 << " covariance is " << acc2.cov(acc3)
          << ", expected " << cov << endl;
 
     return 0;
 }