Index: trunk/tests/Makefile =================================================================== --- trunk/tests/Makefile (revision 556) +++ trunk/tests/Makefile (revision 557) @@ -1,118 +1,119 @@ # Set the following two variables correctly UTESTPP_DIR = /home/igv/Code/UnitTest++ FFTW_LIBDIR = /usr/local/lib OFILES_COMMON = test_utils.o OTESTS_FAST = test_ArrayND.o test_OrthoPoly1D.o test_OrthoPolyND.o \ test_Functors.o test_interpolation.o test_LocalPolyFilterND.o \ test_Distributions1D.o test_InterpolatedDistribution1D.o \ test_EmpiricalCopula.o test_BoxNDScanner.o \ test_QuadraticOrthoPolyND.o test_LocalLogisticRegression.o \ test_randomPermutation.o test_gegenbauerSeriesSum.o \ test_GaussHermiteQuadrature.o test_SpecialFunctions.o \ test_GaussLegendreQuadrature.o test_WeightedStatAccumulator.o \ test_RandomSequenceRepeater.o test_logLikelihoodPeak.o \ test_legendreSeriesSum.o test_solveCubic.o test_HistoND.o \ test_SampleAccumulator.o test_GridInterpolatedDistribution.o \ test_GridAxis.o test_Matrix.o test_StatUtils.o \ test_amiseOptimalBandwidth.o test_OrderedPointND.o \ test_CompositeDistributionND.o test_CompositeBuilder.o \ test_ArchiveIO.o test_Ntuples.o test_LocalQuantileRegression.o \ test_rectangleQuadrature.o test_RegularSampler1D.o \ test_Instantiations.o test_chebyshevSeriesSum.o \ test_dawsonIntegral.o test_NMCombinationSequencer.o \ test_CrossCovarianceAccumulator.o test_StatAccumulatorPair.o \ test_findPeak2D.o test_StatAccumulatorArr.o \ test_goldenSectionSearch.o test_empiricalQuantile.o \ test_StatAccumulator.o test_binomialCoefficient.o \ test_CircularMapper1d.o test_LinInterpolatedTableND.o \ test_StorableMultivariateFunctor.o test_NeymanPearsonWindow1D.o \ test_hermiteSeriesSum.o test_hermiteSeriesSumPhys.o \ test_DiscreteDistributions1D.o test_LocalPolyFilter1D.o \ test_definiteIntegrals.o test_ExpMapper1d.o test_MathUtils.o \ test_MemoizingSymbetaFilterProvider.o test_ResponseMatrix.o \ test_gaussianResponseMatrix.o LogLogQuadratic1D.o \ WeightedTestAccumulator.o TestAccumulator.o test_histoUtils.o \ test_LocalMultiFilter1D.o test_InterpolatedCompositeBuilder.o \ test_DensityAveScanND.o test_Distro1DBuilder.o \ test_findRootNewtonRaphson.o test_Uncertainties.o \ test_ClassicalOrthoPolys1D.o test_ContOrthoPoly1D.o \ - test_findRootUsingBisections.o + test_findRootUsingBisections.o test_orthoPoly1DVProducts.o OTESTS_COMPREHENSIVE = test_ArrayND_cmh.o test_DiscreteBernstein.o \ test_KDTree.o test_rescanArray_cmh.o test_lorpeMise1D.o # OPTIMIZE = -g OPTIMIZE = -std=c++11 -O3 INCLUDES = -I. -I.. -I$(UTESTPP_DIR)/src -I$(FFTW_LIBDIR)/../include CPPFLAGS = $(OPTIMIZE) $(INCLUDES) -Wall -W -Werror -Wno-unused-parameter LIBS = -L../npstat/.libs -lnpstat -L$(UTESTPP_DIR) -lUnitTest++ \ -L$(FFTW_LIBDIR) -L/usr/lib64 -lfftw3 -llapack -lblas \ -L/usr/local/lib -lrk -lgeners -lbz2 -lz -lkstest -lm %.o : %.cc g++ -c $(CPPFLAGS) -fPIC -MD $< -o $@ @sed -i 's,\($*\.o\)[:]*\(.*\),$@: $$\(wildcard\2\)\n\1:\2,g' $*.d # Change the "all" target below to suit your development needs. # Useful targets which can be included are: # # fast fast_run # comprehensive comprehensive_run # regression regression_run # # The tests inside the "fast" target are a subset of tests from the # "comprehensive" target all: fast fast_run # all: fast regression fast_run regression_run # all: comprehensive comprehensive_run regression regression_run fast: test_main.o $(OTESTS_FAST) $(OFILES_COMMON) rm -f $@ g++ $(OPTIMIZE) -fPIC -o $@ $^ $(LIBS) fast_run: fast ./fast comprehensive: test_main.o $(OTESTS_FAST) \ $(OTESTS_COMPREHENSIVE) $(OFILES_COMMON) rm -f $@ g++ $(OPTIMIZE) -fPIC -o $@ $^ $(LIBS) comprehensive_run: comprehensive ./comprehensive regression: regression_run: regression @ echo Running regression tests. PROGRAMS = kdtree_speed.cc quantileBinFromCdf.cc incompleteGamma.cc \ convertToSpherical.cc printPermutations.cc showIOTraits.cc \ hugeNtuple.cc hugeNtupleRead.cc histoStdev.cc effectiveDim.cc \ buildInterpolatedCheck.cc legendreRoots.cc jacobiPolyStats.cc \ - jacobiVProducts.cc jacobiEpsTest.cc gauss2DRandom.cc + jacobiVProducts.cc jacobiEpsTest.cc jacobiEpsPseudoExp.cc \ + gauss2DRandom.cc PROGRAMS += cpp11Random.cc BINARIES = $(PROGRAMS:.cc=) binaries: $(BINARIES) $(BINARIES): % : %.o $(OFILES_COMMON); g++ $(OPTIMIZE) -fPIC -o $@ $^ $(LIBS) clean: rm -f fast comprehensive $(BINARIES) \ *.out core.* *.o *.d *~ *.gsbd *.gsbmf -include test_main.d -include $(OFILES_COMMON:.o=.d) -include $(OTESTS_FAST:.o=.d) -include $(OTESTS_COMPREHENSIVE:.o=.d) -include $(PROGRAMS:.cc=.d) Index: trunk/tests/test_ContOrthoPoly1D.cc =================================================================== --- trunk/tests/test_ContOrthoPoly1D.cc (revision 556) +++ trunk/tests/test_ContOrthoPoly1D.cc (revision 557) @@ -1,92 +1,123 @@ #include "UnitTest++.h" #include "test_utils.hh" #include "npstat/nm/ContOrthoPoly1D.hh" #include "npstat/rng/MersenneTwister.hh" using namespace npstat; using namespace std; inline static int kdelta(const unsigned i, const unsigned j) { return i == j ? 1 : 0; } namespace { TEST(ContOrthoPoly1D_orthonormalization) { const double eps = 1.0e-15; const OrthoPolyMethod method[] = {OPOLY_STIELTJES, OPOLY_LANCZOS}; const unsigned nMethods = sizeof(method)/sizeof(method[0]); const unsigned npoints = 64; const unsigned maxdeg = 10; const unsigned ntries = 10; std::vector points(2U*npoints); std::vector measure; measure.reserve(npoints); MersenneTwister rng; for (unsigned imeth=0; imeth(kdelta(i, j)), d, eps); } } } TEST(ContOrthoPoly1D_cov8) { const double eps = 1.0e-12; MersenneTwister rng; const unsigned npoints = 64; const unsigned maxdeg = 6; const unsigned ntries = 5; const unsigned degtries = 100; std::vector points(npoints); for (unsigned itry=0; itry points(npoints); + + for (unsigned itry=0; itry +#include +#include + +inline std::string time_stamp() +{ + struct tm *current; + time_t now; + + time(&now); + current = localtime(&now); + + char buf[10]; + sprintf(buf, "%02i:%02i:%02i", current->tm_hour, + current->tm_min, current->tm_sec); + return std::string(buf); +} + +#endif // TIME_STAMP_H_ Index: trunk/tests/jacobiEpsPseudoExp.cc =================================================================== --- trunk/tests/jacobiEpsPseudoExp.cc (revision 0) +++ trunk/tests/jacobiEpsPseudoExp.cc (revision 557) @@ -0,0 +1,337 @@ +#include +#include +#include +#include +#include +#include + +#include +#include + +#include +#include + +#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; + } + + 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; + } + + 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 +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 rn; + rn.reserve(perpseudo); + + const unsigned integPoints = GaussLegendreQuadrature::minimalExactRule( + (unsigned)(jac.alpha()) + (unsigned)(jac.beta()) + 2U*maxdeg); + + for (unsigned ips=0; ips 0) + if (ips % reportFrequency == 0) + cout << time_stamp() << " running pseudo exp " << ips << endl; + + rn.clear(); + double tmp; + for (unsigned i=0; i 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 m=0; m<=n; ++m) + { + const double ecov = calculateCovariance( + jac, perpseudo, m, n, i, j, maxdeg, highOrder_cov); + orafile << ecov << ' '; + } + + orafile << '\n'; + return 0; +} Index: trunk/tests/test_orthoPoly1DVProducts.cc =================================================================== --- trunk/tests/test_orthoPoly1DVProducts.cc (revision 0) +++ trunk/tests/test_orthoPoly1DVProducts.cc (revision 557) @@ -0,0 +1,24 @@ +#include "UnitTest++.h" +#include "test_utils.hh" + +#include "npstat/nm/ClassicalOrthoPolys1D.hh" +#include "npstat/nm/GaussLegendreQuadrature.hh" + +#include "npstat/stat/orthoPoly1DVProducts.hh" + +using namespace npstat; +using namespace std; + +namespace { + TEST(oracleEpsCovarianceArray) + { + const unsigned maxdeg = 5; + const unsigned long nPoints = 100; + const bool highOrder = false; + + JacobiOrthoPoly1D jac(3, 4); + const unsigned nQuad = GaussLegendreQuadrature::minimalExactRule(3+4+8*maxdeg); + GaussLegendreQuadrature glq(nQuad); + oracleEpsCovarianceArray(jac, glq, nPoints, maxdeg, highOrder); + } +} Index: trunk/npstat/nm/ContOrthoPoly1D.hh =================================================================== --- trunk/npstat/nm/ContOrthoPoly1D.hh (revision 556) +++ trunk/npstat/nm/ContOrthoPoly1D.hh (revision 557) @@ -1,474 +1,483 @@ #ifndef NPSTAT_CONTORTHOPOLY1D_HH_ #define NPSTAT_CONTORTHOPOLY1D_HH_ /*! // \file ContOrthoPoly1D.hh // // \brief Continuous orthogonal polynomial series in one dimension // generated for a discrete measure // // Author: I. Volobouev // // May 2017 */ #include #include #include #include #include #include "npstat/nm/OrthoPolyMethod.hh" #include "npstat/nm/Matrix.hh" #include "npstat/nm/SimpleFunctors.hh" namespace npstat { class ContOrthoPoly1D { public: // The first element of the pair is the coordinate and the // second measure is the weight (all weights must be non-negative) typedef std::pair MeasurePoint; /** // Main constructor, with obvious arguments. Internally, the // measure will be sorted in the order of increasing weight // and the measure coordinates will be shifted so that their // weighted mean is at 0. */ ContOrthoPoly1D(unsigned maxDegree, const std::vector& measure, OrthoPolyMethod m = OPOLY_STIELTJES); /** Constructor which assumes that all weights are 1.0 */ ContOrthoPoly1D(unsigned maxDegree, const std::vector& coords, OrthoPolyMethod m = OPOLY_STIELTJES); //@{ /** A simple inspector of object properties */ inline unsigned maxDegree() const {return maxdeg_;} inline unsigned long measureLength() const {return measure_.size();} inline MeasurePoint measure(const unsigned index) const {return measure_.at(index);} inline double coordinate(const unsigned index) const {return measure_.at(index).first;} inline double weight(const unsigned index) const {return measure_.at(index).second;} inline double sumOfWeights() const {return wsum_;} inline double sumOfWeightSquares() const {return wsumsq_;} inline bool areAllWeightsEqual() const {return allWeightsEqual_;} inline double minCoordinate() const {return minX_;} inline double maxCoordinate() const {return maxX_;} inline double meanCoordinate() const {return meanX_;} //@} /** Kish's effective sample size for the measure */ double effectiveSampleSize() const; /** Return the value of one of the normalized polynomials */ double poly(unsigned deg, double x) const; /** Return the values of two of the normalized polynomials */ std::pair polyPair( unsigned deg1, unsigned deg2, double x) const; //@{ /** // Return the values of all orthonormal polynomials up to some // maximum degree. Size of the "polyValues" array should be // at least maxdeg + 1. */ void allPolys(double x, unsigned maxdeg, double *polyValues) const; void allPolys(double x, unsigned maxdeg, long double *polyValues) const; //@} /** Return the value of the orthonormal polynomial series */ double series(const double *coeffs, unsigned maxdeg, double x) 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, double *coeffs, unsigned maxdeg) const; /** // Integrated squared error between the given function and the // polynomial series. Parameter "integrationRulePoints" specifies // the parameter "npoints" for the GaussLegendreQuadrature class. // If "integrationRulePoints" is 0, the rule will be picked // automatically in such a way that it integrates a polynomial // of degree maxdeg*2 exactly. */ template double integratedSquaredError(const Functor& fcn, const double *coeffs, unsigned maxdeg, double xmin, double xmax, unsigned integrationRulePoints=0U) const; /** // Squared error between the given function and the polynomial // series, weighted according to the measure */ template double weightedSquaredError(const Functor& fcn, 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. */ double empiricalKroneckerDelta(unsigned deg1, unsigned deg2) const; /** // If the Kronecker deltas were calculated from a sample, they // would be random and would have a covariance matrix. This // is an estimate of the terms in that covariance matrix. The // covariance is between delta_{deg1,deg2} and delta_{deg3,deg4}. */ double empiricalKroneckerCovariance(unsigned deg1, unsigned deg2, unsigned deg3, unsigned deg4) const; /** // Examine the recurrence coefficients. The first element of // the returned pair is alpha, and the second element is beta. */ std::pair recurrenceCoeffs(unsigned deg) const; /** // Generate principal minor of order n of the Jacobi matrix. // n must not exceed "maxDegree" */ Matrix jacobiMatrix(unsigned n) const; /** // Roots of the polynomial with the given degree. Naturally, // the degree argument must not exceed "maxDegree". */ void calculateRoots(double *roots, unsigned degree) const; /** // Integral of an orthonormal polynomials to some power without // the weight (so this is not an inner product with the poly of // degree 0). */ double integratePoly(unsigned degree, unsigned power, double xmin, double xmax) const; /** // Integral of the product of three orthonormal polynomials // without the weight (so this is not an inner product) */ double integrateTripleProduct(unsigned deg1, unsigned deg2, unsigned deg3, double xmin, double xmax) const; /** // Integral of the product of two orthonormal polynomials // with some external weight function (e.g., oracle density). // "EW" in the method name stands for "externally weighted". // // "integrationRulePoints" argument will be passed to the // GaussLegendreQuadrature constructor and must be supported // by that class. */ template double integrateEWPolyProduct(const Functor& weight, unsigned deg1, unsigned deg2, double xmin, double xmax, unsigned integrationRulePoints) const; /** // A measure-weighted average of orthonormal poly values // to some power */ double powerAverage(unsigned deg, unsigned power) const; /** // A measure-weighted average of the product of two orthonormal // poly values to some powers */ double jointPowerAverage(unsigned deg1, unsigned power1, unsigned deg2, unsigned power2) const; /** // A measure-weighted average of the product of several // orthonormal poly values */ double jointAverage(const unsigned* degrees, unsigned nDegrees, bool degreesSortedIncreasingOrder = false) const; //@{ /** // A measure-weighted average that will be remembered internally // so that another call to this function with the same arguments // will return the answer quickly */ double cachedJointAverage(unsigned deg1, unsigned deg2, unsigned deg3, unsigned deg4) const; double cachedJointAverage(unsigned deg1, unsigned deg2, unsigned deg3, unsigned deg4, unsigned deg5, unsigned deg6) const; double cachedJointAverage(unsigned deg1, unsigned deg2, unsigned deg3, unsigned deg4, unsigned deg5, unsigned deg6, unsigned deg7, unsigned deg8) const; //@} /** Covariance between v_{ab} and v_{cd} */ double cov4(unsigned a, unsigned b, unsigned c, unsigned d) const; /** Covariance between v_{ab}*v_{cd} and v_{ef} */ double cov6(unsigned a, unsigned b, unsigned c, unsigned d, unsigned e, unsigned f) const; /** Covariance between v_{ab}*v_{cd} and v_{ef}*v_{gh} */ double cov8(unsigned a, unsigned b, unsigned c, unsigned d, unsigned e, unsigned f, unsigned g, unsigned h) const; /** // Alternative implementation of "cov8". It is slower than "cov8" // but easier to program and verify. Useful mainly for testing "cov8". */ double slowCov8(unsigned a, unsigned b, unsigned c, unsigned d, unsigned e, unsigned f, unsigned g, unsigned h) const; + /** Estimate of eps_{mn} expectation */ double epsExpectation(unsigned m, unsigned n, bool highOrder) const; - double epsCovariance(unsigned m1, unsigned n1, - unsigned m2, unsigned n2, bool highOrder) const; + /** + // Estimate of covariance between eps_{ij} and eps_{mn}. + // + // The result should be identical with "empiricalKroneckerCovariance" + // in case "highOrder" argument is "false". However, this method + // caches results of various calculations of averages and should + // perform better inside cycles over indices. + */ + double epsCovariance(unsigned i, unsigned j, + unsigned m, unsigned n, bool highOrder) const; private: struct Recurrence { inline Recurrence(const long double a, const long double b) : alpha(a), beta(b) { assert(beta > 0.0L); sqrbeta = sqrtl(beta); } long double alpha; long double beta; long double sqrbeta; }; class PolyFcn; friend class PolyFcn; class PolyFcn : public Functor1 { public: inline PolyFcn(const ContOrthoPoly1D& p, const unsigned d1, const unsigned power) : poly(p), deg1(d1), polypow(power) {} inline long double operator()(const long double& x) const { long double p = 1.0L; if (polypow) { p = poly.normpoly(deg1, x); switch (polypow) { case 1U: break; case 2U: p *= p; break; default: p = powl(p, polypow); break; } } return p; } private: const ContOrthoPoly1D& poly; unsigned deg1; unsigned polypow; }; class TripleProdFcn; friend class TripleProdFcn; class TripleProdFcn : public Functor1 { public: inline TripleProdFcn(const ContOrthoPoly1D& p, const unsigned d1, const unsigned d2, const unsigned d3) : poly(p) { degs[0] = d1; degs[1] = d2; degs[2] = d3; degmax = std::max(d1, std::max(d2, d3)); } inline long double operator()(const long double& x) const {return poly.normpolyprod(degs, 3, degmax, x);} private: const ContOrthoPoly1D& poly; unsigned degs[3]; unsigned degmax; }; template class EWPolyProductFcn; template friend class EWPolyProductFcn; template class EWPolyProductFcn : public Functor1 { public: inline EWPolyProductFcn(const ContOrthoPoly1D& p, const Funct& w, const unsigned d1, const unsigned d2) : poly(p), wf(w) { degs[0] = d1; degs[1] = d2; degmax = std::max(d1, d2); } inline long double operator()(const long double& x) const {return poly.normpolyprod(degs,2,degmax,x-poly.meanX_)*wf(x);} private: const ContOrthoPoly1D& poly; const Funct& wf; unsigned degs[2]; unsigned degmax; }; class MemoKey { public: inline MemoKey(const unsigned d0, const unsigned d1, const unsigned d2, const unsigned d3) : nDegs_(4U) { degs_[0] = d0; degs_[1] = d1; degs_[2] = d2; degs_[3] = d3; std::sort(degs_, degs_+nDegs_); } inline MemoKey(const unsigned d0, const unsigned d1, const unsigned d2, const unsigned d3, const unsigned d4, const unsigned d5) : nDegs_(6U) { degs_[0] = d0; degs_[1] = d1; degs_[2] = d2; degs_[3] = d3; degs_[4] = d4; degs_[5] = d5; std::sort(degs_, degs_+nDegs_); } inline MemoKey(const unsigned d0, const unsigned d1, const unsigned d2, const unsigned d3, const unsigned d4, const unsigned d5, const unsigned d6, const unsigned d7) : nDegs_(8U) { degs_[0] = d0; degs_[1] = d1; degs_[2] = d2; degs_[3] = d3; degs_[4] = d4; degs_[5] = d5; degs_[6] = d6; degs_[7] = d7; std::sort(degs_, degs_+nDegs_); } inline const unsigned* degrees() const {return degs_;} inline unsigned nDegrees() const {return nDegs_;} inline bool operator<(const MemoKey& r) const { if (nDegs_ < r.nDegs_) return true; if (r.nDegs_ < nDegs_) return false; for (unsigned i=0; i monicInnerProducts( unsigned degree) const; long double normpoly(unsigned degree, long double x) const; std::pair twonormpoly( unsigned deg1, unsigned deg2, long double x) const; long double normpolyprod(const unsigned* degrees, unsigned nDeg, unsigned degmax, long double x) const; long double normseries(const double *coeffs, unsigned maxdeg, long double x) const; template void getAllPolys(double x, unsigned maxdeg, Numeric *polyValues) const; double getCachedAverage(const MemoKey& key) const; std::vector measure_; std::vector rCoeffs_; mutable std::map cachedAverages_; long double wsum_; long double wsumsq_; double meanX_; double minX_; double maxX_; unsigned maxdeg_; bool allWeightsEqual_; }; } #include "npstat/nm/ContOrthoPoly1D.icc" #endif // NPSTAT_CONTORTHOPOLY1D_HH_ Index: trunk/npstat/stat/orthoPoly1DVProducts.icc =================================================================== --- trunk/npstat/stat/orthoPoly1DVProducts.icc (revision 556) +++ trunk/npstat/stat/orthoPoly1DVProducts.icc (revision 557) @@ -1,453 +1,467 @@ #include #include #include #include +#include namespace npstat { namespace Private { inline void validateNPoints(const char* where, const unsigned long nPoints) { if (!nPoints) { std::string message("In npstat::"); message += where; message += ": sample size must be positive"; throw std::invalid_argument(message); } } inline int kdelta(const UUPair ab) { return ab.first == ab.second ? 1 : 0; } inline bool isZeroPair(const UUPair pair) { return !(pair.first || pair.second); } } template double expectedVProduct(const AbsClassicalOrthoPoly1D& poly, const Quadrature& quad, const unsigned long nPoints, const UUPair ab, const UUPair cd) { Private::validateNPoints("expectedVProduct", nPoints); if (Private::isZeroPair(ab) || Private::isZeroPair(cd)) return 0.0; else return (poly.jointAverage(quad, ab.first, ab.second, cd.first, cd.second) - Private::kdelta(ab)*Private::kdelta(cd))/nPoints; } template inline double expectedVCovariance(const AbsClassicalOrthoPoly1D& poly, const Quadrature& quad, const unsigned long nPoints, const UUPair ab, const UUPair cd) { Private::validateNPoints("expectedVCovariance", nPoints); return expectedVProduct(poly, quad, nPoints, ab, cd); } template double expectedVProduct(const AbsClassicalOrthoPoly1D& poly, const Quadrature& quad, const unsigned long nPoints, const UUPair ab, const UUPair cd, const UUPair ef) { Private::validateNPoints("expectedVProduct", nPoints); if (Private::isZeroPair(ab) || Private::isZeroPair(cd) || Private::isZeroPair(ef)) return 0.0; const double pcdef = poly.jointAverage(quad, cd.first, cd.second, ef.first, ef.second); const double pabef = poly.jointAverage(quad, ab.first, ab.second, ef.first, ef.second); const double pabcd = poly.jointAverage(quad, ab.first, ab.second, cd.first, cd.second); unsigned degs[6]; degs[0] = ab.first; degs[1] = ab.second; degs[2] = cd.first; degs[3] = cd.second; degs[4] = ef.first; degs[5] = ef.second; const double pabcdef = poly.jointAverage(quad, degs, 6); const int deltaprod = Private::kdelta(ab)*Private::kdelta(cd)* Private::kdelta(ef); return (pabcdef - Private::kdelta(ab)*pcdef - Private::kdelta(cd)*pabef - Private::kdelta(ef)*pabcd + 2.0*deltaprod)/nPoints/nPoints; } template inline double expectedVCovariance(const AbsClassicalOrthoPoly1D& poly, const Quadrature& quad, const unsigned long nPoints, const UUPair ab, const UUPair cd, const UUPair ef) { Private::validateNPoints("expectedVCovariance", nPoints); return expectedVProduct(poly, quad, nPoints, ab, cd, ef); } template double expectedVProduct(const AbsClassicalOrthoPoly1D& poly, const Quadrature& quad, const unsigned long nPoints, const UUPair ab, const UUPair cd, const UUPair ef, const UUPair gh) { Private::validateNPoints("expectedVProduct", nPoints); if (Private::isZeroPair(ab) || Private::isZeroPair(cd) || Private::isZeroPair(ef) || Private::isZeroPair(gh)) return 0.0; const double pabcd = poly.jointAverage(quad, ab.first, ab.second, cd.first, cd.second); const double pefgh = poly.jointAverage(quad, ef.first, ef.second, gh.first, gh.second); const double pabef = poly.jointAverage(quad, ab.first, ab.second, ef.first, ef.second); const double pcdgh = poly.jointAverage(quad, cd.first, cd.second, gh.first, gh.second); const double pabgh = poly.jointAverage(quad, ab.first, ab.second, gh.first, gh.second); const double pcdef = poly.jointAverage(quad, cd.first, cd.second, ef.first, ef.second); unsigned degs[8]; degs[0] = ab.first; degs[1] = ab.second; degs[2] = cd.first; degs[3] = cd.second; degs[4] = ef.first; degs[5] = ef.second; degs[6] = gh.first; degs[7] = gh.second; const double pabcdefgh = poly.jointAverage(quad, degs, 8); const double pabcdef = poly.jointAverage(quad, degs, 6); degs[4] = gh.first; degs[5] = gh.second; const double pabcdgh = poly.jointAverage(quad, degs, 6); degs[2] = ef.first; degs[3] = ef.second; const double pabefgh = poly.jointAverage(quad, degs, 6); degs[0] = cd.first; degs[1] = cd.second; const double pcdefgh = poly.jointAverage(quad, degs, 6); const int deltaprod = Private::kdelta(ab)*Private::kdelta(cd)* Private::kdelta(ef)*Private::kdelta(gh); const double tmp = Private::kdelta(ef)*Private::kdelta(gh)*pabcd + Private::kdelta(cd)*Private::kdelta(gh)*pabef + Private::kdelta(cd)*Private::kdelta(ef)*pabgh + Private::kdelta(ab)*Private::kdelta(cd)*pefgh + Private::kdelta(ab)*Private::kdelta(ef)*pcdgh + Private::kdelta(ab)*Private::kdelta(gh)*pcdef; const double term1 = pabcd*pefgh + pabef*pcdgh + pabgh*pcdef + 3.0*deltaprod - tmp; const double term2 = pabcdefgh - 6.0*deltaprod - Private::kdelta(ab)*pcdefgh - Private::kdelta(cd)*pabefgh - Private::kdelta(ef)*pabcdgh - Private::kdelta(gh)*pabcdef - pabcd*pefgh - pabef*pcdgh - pabgh*pcdef + 2.0*tmp; return (term1 + term2/nPoints)/nPoints/nPoints; } template inline double expectedVCovariance(const AbsClassicalOrthoPoly1D& poly, const Quadrature& quad, const unsigned long nPoints, const UUPair ab, const UUPair cd, const UUPair ef, const UUPair gh) { Private::validateNPoints("expectedVCovariance", nPoints); return expectedVProduct(poly, quad, nPoints, ab, cd, ef, gh) - expectedVProduct(poly, quad, nPoints, ab, cd)* expectedVProduct(poly, quad, nPoints, ef, gh); } template inline double sampleVProduct(const AbsClassicalOrthoPoly1D& poly, const Numeric* coords, const unsigned long nPoints, const UUPair ab) { Private::validateNPoints("sampleVProduct", nPoints); assert(coords); if (Private::isZeroPair(ab)) return 0.0; else { unsigned degs[2]; degs[0] = ab.first; degs[1] = ab.second; return poly.jointSampleAverage(coords, nPoints, degs, 2) - Private::kdelta(ab); } } template inline double sampleVProduct(const AbsClassicalOrthoPoly1D& poly, const Numeric* coords, const unsigned long nPoints, const UUPair ab, const UUPair cd) { Private::validateNPoints("sampleVProduct", nPoints); assert(coords); if (Private::isZeroPair(ab) || Private::isZeroPair(cd)) return 0.0; else { unsigned degs1[2]; degs1[0] = ab.first; degs1[1] = ab.second; unsigned degs2[2]; degs2[0] = cd.first; degs2[1] = cd.second; const std::pair& dp = poly.twoJointSampleAverages( coords, nPoints, degs1, 2, degs2, 2); return (dp.first - Private::kdelta(ab))* (dp.second - Private::kdelta(cd)); } } template inline double sampleVProduct(const AbsClassicalOrthoPoly1D& poly, const Numeric* coords, const unsigned long nPoints, const UUPair ab, const UUPair cd, const UUPair ef) { Private::validateNPoints("sampleVProduct", nPoints); assert(coords); if (Private::isZeroPair(ab) || Private::isZeroPair(cd) || Private::isZeroPair(ef)) return 0.0; else return sampleVProduct(poly, coords, nPoints, ab, cd)* sampleVProduct(poly, coords, nPoints, ef); } template inline double sampleVProduct(const AbsClassicalOrthoPoly1D& poly, const Numeric* coords, const unsigned long nPoints, const UUPair ab, const UUPair cd, const UUPair ef, const UUPair gh) { Private::validateNPoints("sampleVProduct", nPoints); assert(coords); if (Private::isZeroPair(ab) || Private::isZeroPair(cd) || Private::isZeroPair(ef) || Private::isZeroPair(gh)) return 0.0; else return sampleVProduct(poly, coords, nPoints, ab, cd)* sampleVProduct(poly, coords, nPoints, ef, gh); } template double oracleEpsExpectation(const AbsClassicalOrthoPoly1D& poly, const Quadrature& quad, const unsigned long nPoints, const unsigned m_in, const unsigned n_in, const bool highOrder) { if (highOrder) { Private::validateNPoints("oracleEpsExpectation", nPoints); long double sum = 0.0L; if (m_in || n_in) { const unsigned m = std::min(m_in, n_in); const unsigned n = std::max(m_in, n_in); for (unsigned k=0; k<=n; ++k) { const double f = k == n ? 1.0 : (k > m ? 1.0 : 2.0); sum += f*poly.jointAverage(quad, k, k, m, n); } if (m == n) sum -= 1.0; else sum -= (poly.jointAverage(quad, m, m, m, n) + poly.jointAverage(quad, m, n, n, n))/2.0; } return sum/nPoints; } else return 0.0; } template double sampleEpsValue(const AbsClassicalOrthoPoly1D& poly, const Numeric* coords, unsigned long nPoints, const unsigned m_in, const unsigned n_in, const bool highOrder) { Private::validateNPoints("sampleEpsValue", nPoints); assert(coords); if (highOrder) { long double sum = 0.0L; if (m_in || n_in) { const unsigned m = std::min(m_in, n_in); const unsigned n = std::max(m_in, n_in); const double vmn = sampleVProduct(poly, coords, nPoints, UUPair(m, n)); if (m == n) sum -= (vmn + vmn*vmn); else { const double vnn = sampleVProduct(poly, coords, nPoints, UUPair(n, n)); const double vmm = sampleVProduct(poly, coords, nPoints, UUPair(m, m)); sum -= (vmn + 0.5*vmn*(vnn + vmm)); } for (unsigned k=0; k<=n; ++k) { const double f = k > m ? 1.0 : 2.0; const double vkmkn = sampleVProduct(poly, coords, nPoints, UUPair(k, m), UUPair(k, n)); sum += f*vkmkn; } } return sum; } else return -sampleVProduct(poly, coords, nPoints, UUPair(m_in, n_in)); } template double oracleEpsCovariance(const AbsClassicalOrthoPoly1D& poly, const Quadrature& quad, const unsigned long nPoints, const unsigned m1_in, const unsigned n1_in, const unsigned m2_in, const unsigned n2_in, const bool highOrder) { Private::validateNPoints("oracleEpsCovariance", nPoints); const bool has0 = (m1_in == 0 && n1_in == 0) || (m2_in == 0 && n2_in == 0); if (has0) return 0.0; if (highOrder) { const unsigned m1 = std::min(m1_in, n1_in); const unsigned n1 = std::max(m1_in, n1_in); const unsigned m2 = std::min(m2_in, n2_in); const unsigned n2 = std::max(m2_in, n2_in); long double sum = 0.0L; // Process the -v_{m1,n1} term sum += expectedVCovariance(poly, quad, nPoints, UUPair(m2, n2), UUPair(m1, n1)); sum += expectedVCovariance(poly, quad, nPoints, UUPair(m2, n2), UUPair(n2, n2), UUPair(m1, n1))/2.0; sum += expectedVCovariance(poly, quad, nPoints, UUPair(m2, n2), UUPair(m2, m2), UUPair(m1, n1))/2.0; for (unsigned k=0; k<=n2; ++k) { const double factor = k > m2 ? 1.0 : 2.0; sum -= factor*expectedVCovariance(poly, quad, nPoints, UUPair(k, m2), UUPair(k, n2), UUPair(m1, n1)); } // Process the term -v_{m1,n1}/2 (v_{n1,n1} + v_{m1,m1}) unsigned idx[2]; idx[0] = n1; idx[1] = m1; for (unsigned ii=0; ii<2; ++ii) { const unsigned mOrN = idx[ii]; sum += expectedVCovariance(poly, quad, nPoints, UUPair(m1, n1), UUPair(mOrN, mOrN), UUPair(m2, n2))/2.0; sum += expectedVCovariance(poly, quad, nPoints, UUPair(m1, n1), UUPair(mOrN, mOrN), UUPair(m2, n2), UUPair(n2, n2))/4.0; sum += expectedVCovariance(poly, quad, nPoints, UUPair(m1, n1), UUPair(mOrN, mOrN), UUPair(m2, n2), UUPair(m2, m2))/4.0; for (unsigned k=0; k<=n2; ++k) { const double factor = k > m2 ? 1.0 : 2.0; sum -= factor/2.0*expectedVCovariance(poly, quad, nPoints, UUPair(m1, n1), UUPair(mOrN, mOrN), UUPair(k, m2), UUPair(k, n2)); } } // Process the sum for (unsigned k1=0; k1<=n1; ++k1) { const double f1 = k1 > m1 ? 1.0 : 2.0; sum -= f1*expectedVCovariance(poly, quad, nPoints, UUPair(k1, m1), UUPair(k1, n1), UUPair(m2, n2)); sum -= f1*expectedVCovariance(poly, quad, nPoints, UUPair(k1, m1), UUPair(k1, n1), UUPair(m2, n2), UUPair(n2, n2))/2.0; sum -= f1*expectedVCovariance(poly, quad, nPoints, UUPair(k1, m1), UUPair(k1, n1), UUPair(m2, n2), UUPair(m2, m2))/2.0; for (unsigned k=0; k<=n2; ++k) { const double factor = k > m2 ? 1.0 : 2.0; sum += f1*factor*expectedVCovariance(poly, quad, nPoints, UUPair(k1, m1), UUPair(k1, n1), UUPair(k, m2), UUPair(k, n2)); } } return sum; } else return expectedVCovariance(poly, quad, nPoints, UUPair(m2_in, n2_in), UUPair(m1_in, n1_in)); } template ArrayND oracleEpsCovarianceArray( const AbsClassicalOrthoPoly1D& poly, const Quadrature& quad, const unsigned long nPoints, const unsigned maxdeg, const bool highOrder) { const unsigned mdp1 = maxdeg + 1; ArrayND arr(mdp1, mdp1, mdp1, mdp1); + unsigned long nDefined = 0; for (unsigned i=0; i<=maxdeg; ++i) for (unsigned j=0; j<=i; ++j) - for (unsigned k=0; k<=j; ++k) + for (unsigned k=0; k<=maxdeg; ++k) for (unsigned m=0; m<=k; ++m) - arr(i, j, k, m) = oracleEpsCovariance( - poly, quad, nPoints, i, j, k, m, highOrder); + if (std::make_pair(j,i) <= std::make_pair(m,k)) + { + arr(j, i, m, k) = oracleEpsCovariance( + poly, quad, nPoints, j, i, m, k, highOrder); + ++nDefined; + } - unsigned d[4]; for (unsigned i=0; i<=maxdeg; ++i) for (unsigned j=0; j<=maxdeg; ++j) for (unsigned k=0; k<=maxdeg; ++k) for (unsigned m=0; m<=maxdeg; ++m) - { - d[0] = m; - d[1] = k; - d[2] = j; - d[3] = i; - std::sort(d, d+4); - - if (!(d[0] == m && d[1] == k && d[2] == j && d[3] == i)) - arr(i, j, k, m) = arr(d[3], d[2], d[1], d[0]); - } + if (!(j<=i && m<=k && std::make_pair(j,i) <= std::make_pair(m,k))) + { + unsigned i1 = i; + unsigned j1 = j; + unsigned k1 = k; + unsigned m1 = m; + if (!(j1 <= i1)) + std::swap(i1, j1); + if (!(m1 <= k1)) + std::swap(m1, k1); + if (!(std::make_pair(j1,i1) <= std::make_pair(m1,k1))) + { + std::swap(m1, j1); + std::swap(k1, i1); + } + arr(j, i, m, k) = arr(j1, i1, m1, k1); + ++nDefined; + } + assert(nDefined == arr.length()); return arr; } }