Index: trunk/tests/jacobiEpsPseudoExp.cc
===================================================================
--- trunk/tests/jacobiEpsPseudoExp.cc	(revision 560)
+++ trunk/tests/jacobiEpsPseudoExp.cc	(revision 561)
@@ -1,337 +1,325 @@
 #include <string>
 #include <sstream>
 #include <iostream>
 #include <fstream>
 #include <cerrno>
 #include <utility>
 
 #include <string.h>
 #include <unistd.h>
 
 #include <sys/stat.h>
 #include <sys/types.h>
 
 #include "test_utils.hh"
 #include "time_stamp.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 unsigned maxdeg, const bool highOrder)
 {
     const unsigned sumdeg = (unsigned)(jac.alpha()) + (unsigned)(jac.beta()) + 4U*maxdeg;
     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 unsigned maxdeg, const bool highOrder)
 {
     const unsigned sumdeg = (unsigned)(jac.alpha()) + (unsigned)(jac.beta()) + 8U*maxdeg;
     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 char* prefix)
 {
     ostringstream os;
     os << prefix << "_eps_" << m << '_' << n;
     return os.str();
 }
 
 
 static string cov_columns(const unsigned maxdeg)
 {
     ostringstream os;
     unsigned colnum = 0;
 
     for (unsigned j=0; j<=maxdeg; ++j)
         for (unsigned i=0; i<=j; ++i)
-            for (unsigned n=0; n<=i; ++n)
-                for (unsigned m=0; m<=n; ++m, ++colnum)
-                {
-                    if (colnum)
-                        os << ' ';
-                    os << "cov_" << m << '_' << n << '_' << i << '_' << j;
-                }
+            for (unsigned n=0; n<=maxdeg; ++n)
+                for (unsigned m=0; m<=n; ++m)
+                    if (std::make_pair(m,n) <= std::make_pair(i,j))
+                    {
+                        if (colnum++)
+                            os << ' ';
+                        os << "cov_" << m << '_' << n << '_' << i << '_' << j;
+                    }
 
     return os.str();
 }
 
 
 static string pseudo_exp_columns(const unsigned maxdeg)
 {
     ostringstream os;
     unsigned colnum = 0;
 
     for (unsigned j=0; j<=maxdeg; ++j)
         for (unsigned i=0; i<=j; ++i, ++colnum)
         {
             if (colnum)
                 os << ' ';
             os << eps_name(i, j, "real") << ' '
                << eps_name(i, j, "ora") << ' '
-               << eps_name(i, j, "est") << ' '
-               << "pull_" << i << '_' << j;
+               << eps_name(i, j, "est");
         }
 
     os << ' ' << cov_columns(maxdeg);
     return os.str();
 }
 
 
 static string ora_columns(const unsigned maxdeg)
 {
     ostringstream os;
     unsigned colnum = 0;
 
     for (unsigned j=0; j<=maxdeg; ++j)
         for (unsigned i=0; i<=j; ++i, ++colnum)
         {
             if (colnum)
                 os << ' ';
             os << "eps_" << i << '_' << j;
         }
 
     os << ' ' << cov_columns(maxdeg);
     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;
 }
 
 
 #define dump_parameter(str, value) do {             \
     str << "set parameters(" << #value << ") " << value << '\n';\
 } while(0);
 
 
 int main(int argc, char const* argv[])
 {
     const unsigned reportFrequency = 100;
 
     if (argc != 9)
     {
         cout << "Usage: " << parse_progname(argv[0])
              << " alpha beta pointsPerPseudo nPseudo maxdeg highOrder_eps highOrder_cov dirname" << 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 maxdeg;
     if (parse_unsigned(argv[5], &maxdeg)) exit(1);
     if (!maxdeg)
     {
         cerr << "maxdeg must be positive" << endl;
         exit(1);
     }
 
     bool highOrder_eps, highOrder_cov;
     if (parse_bool(argv[6], &highOrder_eps)) exit(1);
     if (parse_bool(argv[7], &highOrder_cov)) exit(1);
 
     // Make a new directory and go there
     const char* dirname = argv[8];
     if (mkdir(dirname, S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH))
     {
         cerr << "Failed to make directory \"" << dirname << "\": "
              << strerror(errno) << endl;
         exit(1);
     }
     if (chdir(dirname))
     {
         cerr << "Failed to change working directory to \"" << dirname << "\": "
              << strerror(errno) << endl;
         exit(1);
     }
 
     // Print the command used and parsed argument values
     {
         ofstream cmdfile("parameters.tcl");
         assert(cmdfile.is_open());
         cmdfile << "# " << argv[0];
         for (int i=1; i<argc; ++i)
             cmdfile << ' ' << argv[i];
         cmdfile << '\n';
 
         cmdfile << "set parameters(executable) \"" << argv[0] << "\"\n";
         dump_parameter(cmdfile, alpha);
         dump_parameter(cmdfile, beta);
         dump_parameter(cmdfile, perpseudo);
         dump_parameter(cmdfile, npseudo);
         dump_parameter(cmdfile, maxdeg);
         dump_parameter(cmdfile, highOrder_eps);
         dump_parameter(cmdfile, highOrder_cov);
         cmdfile << "set parameters(dirname) \"" << dirname << "\"\n";
     }
 
     // File for saving results of pseudo experiments
     ofstream outfile("pseudo.txt");
     assert(outfile.is_open());
     outfile << "# " << pseudo_exp_columns(maxdeg) << '\n';
     outfile.precision(10);
 
     // Various objects needed for calculations
     Beta1D distro(-1.0, 2.0, beta+1, alpha+1);
     MersenneTwister rng;
     JacobiOrthoPoly1D jac(alpha, beta);
     OrthoPoly1DWeight jacWeight(jac);
     vector<double> rn;
     rn.reserve(perpseudo);
 
     const unsigned integPoints = GaussLegendreQuadrature::minimalExactRule(
         (unsigned)(jac.alpha()) + (unsigned)(jac.beta()) + 2U*maxdeg);
 
     for (unsigned ips=0; ips<npseudo; ++ips)
     {
         if (reportFrequency > 0)
             if (ips % reportFrequency == 0)
                 cout << time_stamp() << " running pseudo exp " << ips << endl;
 
         rn.clear();
         double tmp;
         for (unsigned i=0; i<perpseudo; ++i)
         {
             distro.random(rng, &tmp);
             rn.push_back(tmp);
         }
 
         // Calculate actual and estimated eps values
         ContOrthoPoly1D samplePoly(maxdeg, rn);
 
         for (unsigned j=0; j<=maxdeg; ++j)
             for (unsigned i=0; i<=j; ++i)
             {
                 // Actual eps value for this sample (orthonormality uncertainty)
                 const double real_eps = samplePoly.integrateEWPolyProduct(
                     jacWeight, i, j, -1.0, 1.0, integPoints) - kdelta(i, j);
                 outfile << real_eps << ' ';
 
                 // Oracle eps approximation for this sample
                 const double ora_eps = sampleEpsValue(
                     jac, &rn[0], perpseudo, i, j, highOrder_eps);
                 outfile << ora_eps << ' ';
 
                 // Chaos poly eps approximation for this sample
                 const double est_eps = samplePoly.epsExpectation(i, j, highOrder_eps);
                 outfile << est_eps << ' ';
-
-                // Calculate the pull
-                const double cov = samplePoly.epsCovariance(i, j, i, j, highOrder_cov);
-                if (i || j)
-                {
-                    assert(cov > 0.0);
-                    outfile << (real_eps - est_eps)/sqrt(cov) << ' ';
-                }
-                else
-                {
-                    assert(cov == 0.0);
-                    outfile << 0.0 << ' ';
-                }
             }
 
         // Covariances estimated from the sample
         for (unsigned j=0; j<=maxdeg; ++j)
             for (unsigned i=0; i<=j; ++i)
                 for (unsigned n=0; n<=maxdeg; ++n)
                    for (unsigned m=0; m<=n; ++m)
                         // Calculate only if (m,n) <= (i,j)
                         if (std::make_pair(m,n) <= std::make_pair(i,j))
                         {
                             const double cov = samplePoly.epsCovariance(
                                 m, n, i, j, highOrder_cov);
                             outfile << cov << ' ';
                         }
 
         outfile << '\n';
     }
 
     // File for saving oracle predictions
     ofstream orafile("oracle.txt");
     assert(orafile.is_open());
     orafile << "# " << ora_columns(maxdeg) << '\n';
     orafile.precision(10);
 
     // Expected eps values with O(N^{-1}) or O(N^{-3/2}) approximation
     for (unsigned j=0; j<=maxdeg; ++j)
         for (unsigned i=0; i<=j; ++i)
         {
             const double exp_eps = calculateExpectation(
                 jac, perpseudo, i, j, maxdeg, highOrder_eps);
             orafile << exp_eps << ' ';
         }
 
     // Expected eps covariances
     for (unsigned j=0; j<=maxdeg; ++j)
         for (unsigned i=0; i<=j; ++i)
-            for (unsigned n=0; n<=i; ++n)
+            for (unsigned n=0; n<=maxdeg; ++n)
                 for (unsigned m=0; m<=n; ++m)
-                {
-                    const double ecov = calculateCovariance(
-                        jac, perpseudo, m, n, i, j, maxdeg, highOrder_cov);
-                    orafile << ecov << ' ';
-                }
+                    if (std::make_pair(m,n) <= std::make_pair(i,j))
+                    {
+                        const double ecov = calculateCovariance(
+                            jac, perpseudo, m, n, i, j, maxdeg, highOrder_cov);
+                        orafile << ecov << ' ';
+                    }
 
     orafile << '\n';
     return 0;
 }