Index: trunk/tests/test_ContOrthoPoly1D.cc
===================================================================
--- trunk/tests/test_ContOrthoPoly1D.cc	(revision 864)
+++ trunk/tests/test_ContOrthoPoly1D.cc	(revision 865)
@@ -1,406 +1,406 @@
 #include "UnitTest++.h"
 #include "test_utils.hh"
 
 #include "npstat/nm/ContOrthoPoly1D.hh"
 #include "npstat/nm/ClassicalOrthoPolys1D.hh"
 #include "npstat/nm/FejerQuadrature.hh"
 #include "npstat/nm/StorablePolySeries1D.hh"
 
 #include "npstat/rng/MersenneTwister.hh"
 
 #include "npstat/stat/Distributions1D.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<double> points(2U*npoints);
         std::vector<ContOrthoPoly1D::MeasurePoint> measure;
         measure.reserve(npoints);
         std::vector<double> coeffs(maxdeg + 1);
 
         MersenneTwister rng;
 
         for (unsigned imeth=0; imeth<nMethods; ++imeth)
             for (unsigned itry=0; itry<ntries; ++itry)
             {
                 const double xmin = rng()*5.0 - 6.0;
                 const double xmax = rng()*2.0 + 1.0;
                 const double range = xmax - xmin;
                 rng.run(&points[0], 2U*npoints, 2U*npoints);
 
                 measure.clear();
                 for (unsigned i=0; i<npoints; ++i)
                 {
                     ContOrthoPoly1D::MeasurePoint p(points[2*i]*range + xmin, points[2*i+1]);
                     measure.push_back(p);
                 }
                 ContOrthoPoly1D poly(maxdeg, measure, method[imeth]);
 
                 for (unsigned i=0; i<=maxdeg; ++i)
                     for (unsigned j=0; j<=maxdeg; ++j)
                     {
                         const double d = poly.empiricalKroneckerDelta(i, j);
                         CHECK_CLOSE(static_cast<double>(kdelta(i, j)), d, eps);
                     }
 
                 measure.clear();
                 for (unsigned i=0; i<npoints; ++i)
                 {
                     ContOrthoPoly1D::MeasurePoint p(points[2*i]*range + xmin, 1.0);
                     measure.push_back(p);
                 }
                 ContOrthoPoly1D poly2(maxdeg, measure, method[imeth]);
 
                 poly2.weightCoeffs(&coeffs[0], maxdeg);
                 for (unsigned i=0; i<=maxdeg; ++i)
                     CHECK_CLOSE(static_cast<double>(kdelta(i, 0)), coeffs[i], eps);
             }
     }
 
 #ifdef USE_LAPACK_QUAD
     TEST(ContOrthoPoly1D_orthonormalization_quad_STIELTJES)
     {
         const lapack_double eps = FLT128_EPSILON*100.0;
 
         const OrthoPolyMethod method = OPOLY_STIELTJES;
 
         const unsigned npoints = 64;
         const unsigned maxdeg = 10;
         const unsigned ntries = 5;
         const double xmin = -1.0;
         const double xmax = 1.0;
         const double range = xmax - xmin;
 
         std::vector<double> points(2U*npoints);
         std::vector<ContOrthoPoly1D::MeasurePoint> measure;
         measure.reserve(npoints);
 
         MersenneTwister rng;
 
         for (unsigned itry=0; itry<ntries; ++itry)
         {
             rng.run(&points[0], 2U*npoints, 2U*npoints);
 
             measure.clear();
             for (unsigned i=0; i<npoints; ++i)
             {
                 ContOrthoPoly1D::MeasurePoint p(points[2*i]*range + xmin, points[2*i+1]);
                 measure.push_back(p);
             }
             ContOrthoPoly1D poly(maxdeg, measure, method);
 
             for (unsigned i=0; i<=maxdeg; ++i)
                 for (unsigned j=0; j<=maxdeg; ++j)
                 {
                     const lapack_double d = poly.empiricalKroneckerDelta(i, j);
                     CHECK_CLOSE(static_cast<lapack_double>(kdelta(i, j)), d, eps);
                 }
         }
     }
 
     TEST(ContOrthoPoly1D_orthonormalization_quad_LANCZOS)
     {
         const lapack_double eps = FLT128_EPSILON*100.0;
 
         const OrthoPolyMethod method = OPOLY_LANCZOS;
             
         const unsigned npoints = 64;
         const unsigned maxdeg = 10;
         const unsigned ntries = 5;
         const double xmin = -1.0;
         const double xmax = 1.0;
         const double range = xmax - xmin;
 
         std::vector<double> points(2U*npoints);
         std::vector<ContOrthoPoly1D::MeasurePoint> measure;
         measure.reserve(npoints);
 
         MersenneTwister rng;
 
         for (unsigned itry=0; itry<ntries; ++itry)
         {
             rng.run(&points[0], 2U*npoints, 2U*npoints);
 
             measure.clear();
             for (unsigned i=0; i<npoints; ++i)
             {
                 ContOrthoPoly1D::MeasurePoint p(points[2*i]*range + xmin, points[2*i+1]);
                 measure.push_back(p);
             }
             ContOrthoPoly1D poly(maxdeg, measure, method);
 
             for (unsigned i=0; i<=maxdeg; ++i)
                 for (unsigned j=0; j<=maxdeg; ++j)
                 {
                     const lapack_double d = poly.empiricalKroneckerDelta(i, j);
                     CHECK_CLOSE(static_cast<lapack_double>(kdelta(i, j)), d, eps);
                 }
         }
     }
 #endif
 
     TEST(ContOrthoPoly1D_weightCoeffs)
     {
         const double eps = 1.0e-7;
 
         const unsigned maxdeg = 10;
         const unsigned ntries = 10;
         const unsigned npoints = maxdeg + 1U;
 
         std::vector<double> points(2U*npoints);
         std::vector<ContOrthoPoly1D::MeasurePoint> measure;
         measure.reserve(npoints);
         std::vector<double> coeffs(maxdeg + 1);
 
         MersenneTwister rng;
 
         for (unsigned itry=0; itry<ntries; ++itry)
         {
             const double xmin = rng()*5.0 - 6.0;
             const double xmax = rng()*2.0 + 1.0;
             const double range = xmax - xmin;
             rng.run(&points[0], 2U*npoints, 2U*npoints);
 
             measure.clear();
             for (unsigned i=0; i<npoints; ++i)
             {
                 ContOrthoPoly1D::MeasurePoint p(points[2*i]*range + xmin, points[2*i+1]);
                 measure.push_back(p);
             }
             ContOrthoPoly1D poly(maxdeg, measure, OPOLY_LANCZOS);
             poly.weightCoeffs(&coeffs[0], maxdeg);
             CPP11_auto_ptr<StorablePolySeries1D> poly2(
                 poly.makeStorablePolySeries(xmin, xmax));
 
             for (unsigned deg=0; deg<=maxdeg; ++deg)
                 for (unsigned ipow=0; ipow<4; ++ipow)
                 {
                     const double i1 = poly.integratePoly(deg, ipow, xmin, xmax);
                     const double i2 = poly2->integratePoly(deg, ipow);
                     if (fabs(i2) > 1.0e-10)
                         CHECK_CLOSE(1.0, i1/i2, 1.0e-10);
                 }
 
             for (unsigned i=0; i<npoints; ++i)
             {
                 const double x = measure[i].first;
                 const double w = measure[i].second;
                 const double fvalue = poly.series(&coeffs[0], maxdeg, x);
                 const double f2 = poly2->series(&coeffs[0], maxdeg, x);
                 CHECK_CLOSE(w, fvalue, 1000*eps);
-                CHECK_CLOSE(fvalue, f2, 20*eps);
+                CHECK_CLOSE(fvalue, f2, 50*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<double> points(npoints);
 
         for (unsigned itry=0; itry<ntries; ++itry)
         {
             rng.run(&points[0], npoints, npoints);
             ContOrthoPoly1D poly(maxdeg, points);
 
             for (unsigned ideg=0; ideg<degtries; ++ideg)
             {
                 unsigned d[8];
                 bool has0 = true;
                 while (has0)
                 {
                     for (unsigned i=0; i<8; ++i)
                         d[i] = (maxdeg + 0.99)*rng();
                     has0 = (d[0] == 0 && d[1] == 0) ||
                            (d[2] == 0 && d[3] == 0) ||
                            (d[4] == 0 && d[5] == 0) ||
                            (d[6] == 0 && d[7] == 0);
                 }
                 const double cov8 = poly.cov8(d[0],d[1],d[2],d[3],d[4],d[5],d[6],d[7]);
                 const double check = poly.slowCov8(d[0],d[1],d[2],d[3],d[4],d[5],d[6],d[7]);
                 CHECK_CLOSE(check, cov8, eps);
             }
         }
     }
 
     TEST(ContOrthoPoly1D_integrateEWPolyProduct)
     {
         const double eps = 1.0e-12;
 
         MersenneTwister rng;
         const unsigned npoints = 64;
         const unsigned maxdeg = 6;
         const unsigned ntries = 5;
         const unsigned nInteg = 64;
         std::vector<double> points(npoints);
         std::vector<double> aves(maxdeg + 1);
         Beta1D beta(0.0, 1.0, 2.0, 2.0);
         GaussLegendreQuadrature glq(nInteg);
  
         for (unsigned itry=0; itry<ntries; ++itry)
         {
             rng.run(&points[0], npoints, npoints);
             ContOrthoPoly1D poly(maxdeg, points);
 
             DensityFunctor1D df(beta);
             const Matrix<double>& mat = poly.extWeightProductAverages(
                 df, glq, maxdeg, 0.0, 1.0);
             poly.extWeightAverages(df, glq, maxdeg, 0.0, 1.0, &aves[0]);
 
             for (unsigned i=0; i<=maxdeg; ++i)
                 for (unsigned j=0; j<=maxdeg; ++j)
                 {
                     const double integ = poly.integrateEWPolyProduct(
                         DensityFunctor1D(beta), i, j, 0.0, 1.0, nInteg);
                     CHECK_CLOSE(integ, mat[i][j], eps);
                 }
 
             for (unsigned j=0; j<=maxdeg; ++j)
                 CHECK_CLOSE(aves[j], mat[0][j], eps);
         }
     }
 
     TEST(ContOrthoPoly1D_epsCovariance)
     {
         const double eps = 1.0e-15;
         const unsigned npoints = 64;
         const unsigned maxdeg = 6;
         const unsigned ntries = 5;
 
         MersenneTwister rng;
         std::vector<double> points(npoints);
 
         for (unsigned itry=0; itry<ntries; ++itry)
         {
             rng.run(&points[0], npoints, npoints);
             ContOrthoPoly1D poly(maxdeg, points);
 
             for (unsigned i=0; i<=maxdeg; ++i)
                 for (unsigned j=0; j<=i; ++j)
                 {
                     CHECK(poly.epsExpectation(i, j, false) == 0.0);
 
                     for (unsigned k=0; k<=maxdeg; ++k)
                         for (unsigned m=0; m<=k; ++m)
                         {
                             const double eps1 = poly.empiricalKroneckerCovariance(i, j, k, m);
                             const double eps2 = poly.epsCovariance(i, j, k, m, false);
                             CHECK_CLOSE(eps1, eps2, eps);
                         }
                 }
         }
     }
 
     // Test that we can reproduce Jacobi polynomials using ContOrthoPoly1D
     TEST(ContOrthoPoly1D_jacobi_1)
     {
         const double eps = 1.0e-13;
 
         const unsigned alpha = 3;
         const unsigned beta = 4;
         const unsigned maxdeg = 5;
         const unsigned npoints = 64;
 
         JacobiOrthoPoly1D jac(alpha, beta);
         OrthoPoly1DWeight weight(jac);
         OrthoPoly1DDeg poly3(jac, 3);
 
         const unsigned npt = FejerQuadrature::minimalExactRule(2*maxdeg+alpha+beta);
         FejerQuadrature fej(npt);
         ContOrthoPoly1D poly(maxdeg, fej.weightedIntegrationPoints(weight, -1.0, 1.0), OPOLY_LANCZOS);
 
         for (unsigned i=0; i<npoints; ++i)
         {
             const double x = test_rng()*2.0 - 1.0;
             for (unsigned deg=0; deg<=maxdeg; ++deg)
                 CHECK_CLOSE(jac.poly(deg, x), poly.poly(deg, x), eps);
             CHECK_CLOSE(jac.poly(3, x), poly3(x), eps);
         }
 
         // Test "sampleAverages" method
         const unsigned nrandom = 100;
         std::vector<double> r(nrandom);
         for (unsigned i=0; i<nrandom; ++i)
             r[i] = -1.0 + 2.0*test_rng();
         std::vector<double> result(maxdeg + 1);
         poly.sampleAverages(&r[0], r.size(), &result[0], maxdeg);
 
         Matrix<long double> mat(maxdeg + 1, maxdeg + 1, 0);
         std::vector<long double> acc(maxdeg + 1, 0.0L);
         std::vector<long double> tmp(maxdeg + 1);
         for (unsigned i=0; i<nrandom; ++i)
         {
             poly.allPolys(r[i], maxdeg, &tmp[0]);
             for (unsigned deg=0; deg<=maxdeg; ++deg)
             {
                 acc[deg] += tmp[deg];
                 for (unsigned k=0; k<=maxdeg; ++k)
                 {
                     const double prod = tmp[deg]*tmp[k];
                     mat[deg][k] += prod;
                 }
             }
         }
         for (unsigned deg=0; deg<=maxdeg; ++deg)
         {
             acc[deg] /= nrandom;
             CHECK_CLOSE(static_cast<double>(acc[deg]), result[deg], eps);
         }
 
         // Test "pairwiseSampleAverages" method
         mat /= nrandom;
         Matrix<double> averages = poly.sampleProductAverages(
             &r[0], r.size(), maxdeg);
         Matrix<double> refmat(mat);
         CHECK((averages - refmat).maxAbsValue() < eps);
     }
 
     TEST(ContOrthoPoly1D_jacobi_2)
     {
         typedef ContOrthoPoly1D::PreciseQuadrature Quadrature;
 
         const double eps = 1.0e-13;
 
         const unsigned alpha = 3;
         const unsigned beta = 4;
         const unsigned maxdeg = 5;
         const unsigned npoints = 64;
 
         JacobiOrthoPoly1D jac(alpha, beta);
         OrthoPoly1DWeight weight(jac);
 
         const unsigned npt = Quadrature::minimalExactRule(2*maxdeg+alpha+beta);
         Quadrature glq(npt);
         ContOrthoPoly1D poly(maxdeg, glq.weightedIntegrationPoints(weight, -1.0, 1.0), OPOLY_LANCZOS);
 
         for (unsigned i=0; i<npoints; ++i)
         {
             const double x = test_rng()*2.0 - 1.0;
             for (unsigned deg=0; deg<=maxdeg; ++deg)
                 CHECK_CLOSE(jac.poly(deg, x), poly.poly(deg, x), eps);
         }
     }
 }
Index: trunk/npstat/nm/SimpleFunctors.hh
===================================================================
--- trunk/npstat/nm/SimpleFunctors.hh	(revision 864)
+++ trunk/npstat/nm/SimpleFunctors.hh	(revision 865)
@@ -1,966 +1,994 @@
 #ifndef NPSTAT_SIMPLEFUNCTORS_HH_
 #define NPSTAT_SIMPLEFUNCTORS_HH_
 
 /*!
 // \file SimpleFunctors.hh
 //
 // \brief Interface definitions and concrete simple functors for
 //        a variety of functor-based calculations
 //
 // Author: I. Volobouev
 //
 // March 2009
 */
 
 #include <cmath>
 #include <cfloat>
 #include <cassert>
 #include <utility>
 #include <stdexcept>
 
 // The following include is needed for std::sqrt(lapack_double)
 #include "npstat/nm/lapack_double.h"
 
 namespace npstat {
     /** Base class for a functor that takes no arguments */
     template <typename Result>
     struct Functor0
     {
         typedef Result result_type;
 
         inline virtual ~Functor0() {}
         virtual Result operator()() const = 0;
     };
 
     /** Copyable reference to Functor0 */
     template <typename Result>
     class Functor0RefHelper : public Functor0<Result>
     {
     public:
         inline Functor0RefHelper(const Functor0<Result>& ref) : ref_(ref) {}
         inline virtual ~Functor0RefHelper() {}
         inline virtual Result operator()() const {return ref_();}
 
     private:
         const Functor0<Result>& ref_;
     };
 
     /** Convenience function for creating Functor0RefHelper instances */
     template <typename Result>
     inline Functor0RefHelper<Result> Functor0Ref(const Functor0<Result>& ref)
     {
         return Functor0RefHelper<Result>(ref);
     }
 
     /** Base class for a functor that takes a single argument */
     template <typename Result, typename Arg1>
     struct Functor1
     {
         typedef Result result_type;
         typedef Arg1 first_argument_type;
 
         inline virtual ~Functor1() {}
         virtual Result operator()(const Arg1&) const = 0;
     };
 
     /** Copyable reference to Functor1 */
     template <typename Result, typename Arg1>
     class Functor1RefHelper : public Functor1<Result, Arg1>
     {
     public:
         inline Functor1RefHelper(const Functor1<Result, Arg1>& ref) : ref_(ref) {}
         inline virtual ~Functor1RefHelper() {}
         inline virtual Result operator()(const Arg1& x1) const {return ref_(x1);}
 
     private:
         const Functor1<Result, Arg1>& ref_;
     };
 
     /** Convenience function for creating Functor1RefHelper instances */
     template <typename Result, typename Arg1>
     inline Functor1RefHelper<Result,Arg1> Functor1Ref(const Functor1<Result,Arg1>& ref)
     {
         return Functor1RefHelper<Result,Arg1>(ref);
     }
 
     /** Base class for a functor that takes two arguments */
     template <typename Result, typename Arg1, typename Arg2>
     struct Functor2
     {
         typedef Result result_type;
         typedef Arg1 first_argument_type;
         typedef Arg2 second_argument_type;
 
         inline virtual ~Functor2() {}
         virtual Result operator()(const Arg1&, const Arg2&) const = 0;
     };
 
     /** Copyable reference to Functor2 */
     template <typename Result, typename Arg1, typename Arg2>
     class Functor2RefHelper : public Functor2<Result, Arg1, Arg2>
     {
     public:
         inline Functor2RefHelper(const Functor2<Result, Arg1, Arg2>& ref) : ref_(ref) {}
         inline virtual ~Functor2RefHelper() {}
         inline virtual Result operator()(const Arg1& x1, const Arg2& x2) const
             {return ref_(x1, x2);}
 
     private:
         const Functor2<Result, Arg1, Arg2>& ref_;
     };
 
     /** Convenience function for creating Functor2RefHelper instances */
     template <typename Result, typename Arg1, typename Arg2>
     inline Functor2RefHelper<Result,Arg1,Arg2> Functor2Ref(
         const Functor2<Result,Arg1,Arg2>& ref)
     {
         return Functor2RefHelper<Result,Arg1,Arg2>(ref);
     }
 
     /** Base class for a functor that takes three arguments */
     template <typename Result, typename Arg1, typename Arg2, typename Arg3>
     struct Functor3
     {
         typedef Result result_type;
         typedef Arg1 first_argument_type;
         typedef Arg2 second_argument_type;
         typedef Arg3 third_argument_type;
 
         inline virtual ~Functor3() {}
         virtual Result operator()(const Arg1&,const Arg2&,const Arg3&) const=0;
     };
 
     /** Copyable reference to Functor3 */
     template <typename Result, typename Arg1, typename Arg2, typename Arg3>
     class Functor3RefHelper : public Functor3<Result, Arg1, Arg2, Arg3>
     {
     public:
         inline Functor3RefHelper(const Functor3<Result, Arg1, Arg2, Arg3>& ref) : ref_(ref) {}
         inline virtual ~Functor3RefHelper() {}
         inline virtual Result operator()(const Arg1& x1, const Arg2& x2, const Arg3& x3) const
             {return ref_(x1, x2, x3);}
 
     private:
         const Functor3<Result, Arg1, Arg2, Arg3>& ref_;
     };
 
     /** Convenience function for creating Functor3RefHelper instances */
     template <typename Result, typename Arg1, typename Arg2, typename Arg3>
     inline Functor3RefHelper<Result,Arg1,Arg2,Arg3> Functor3Ref(
         const Functor3<Result,Arg1,Arg2,Arg3>& ref)
     {
         return Functor3RefHelper<Result,Arg1,Arg2,Arg3>(ref);
     }
 
     /** Base class for a functor that returns a pair */
     template <typename Numeric>
     struct PairFunctor
     {
         typedef std::pair<Numeric,Numeric> result_type;
         typedef Numeric first_argument_type;
 
         inline virtual ~PairFunctor() {}
         virtual std::pair<Numeric,Numeric> operator()(const Numeric&) const=0;
     };
 
     /** Copyable reference to PairFunctor */
     template <typename Numeric>
     class PairFunctorRefHelper : public PairFunctor<Numeric>
     {
     public:
         inline PairFunctorRefHelper(const PairFunctor<Numeric>& ref) : ref_(ref) {}
         inline virtual ~PairFunctorRefHelper() {}
         inline virtual std::pair<Numeric,Numeric> operator()(const Numeric& c) const
             {return ref_(c);}
 
     private:
         const PairFunctor<Numeric>& ref_;
     };
 
     /** Convenience function for creating PairFunctorRefHelper instances */
     template <typename Numeric>
     inline PairFunctorRefHelper<Numeric> PairFunctorRef(
         const PairFunctor<Numeric>& ref)
     {
         return PairFunctorRefHelper<Numeric>(ref);
     }
 
     /** A simple functor which returns a copy of its argument */
     template <typename Result>
     struct Same : public Functor1<Result, Result>
     {
         inline virtual ~Same() {}
 
         inline Result operator()(const Result& a) const {return a;}
 
         inline bool operator==(const Same&) const {return true;}
         inline bool operator!=(const Same&) const {return false;}
     };
 
     /** A simple functor which adds a constant to its argument */
     template <typename Result>
     class Shift : public Functor1<Result, Result>
     {
     public:
         inline Shift(const Result& v) : value_(v) {}
         inline virtual ~Shift() {}
 
         inline Result operator()(const Result& a) const {return a + value_;}
 
         inline bool operator==(const Shift& r) const
             {return value_ == r.value_;}
         inline bool operator!=(const Shift& r) const
             {return !(*this == r);}
 
     private:
         Shift();
         Result value_;
     };
 
     /** A simple functor which returns a reference to its argument */
     template <typename Result>
     struct SameRef : public Functor1<const Result&, Result>
     {
         inline virtual ~SameRef() {}
 
         inline const Result& operator()(const Result& a) const {return a;}
 
         inline bool operator==(const SameRef&) const {return true;}
         inline bool operator!=(const SameRef&) const {return false;}
     };
 
     /**
     // Simple functor which ignores is arguments and instead
     // builds the result using the default constructor of the result type
     */
     template <typename Result>
     struct DefaultConstructor0 : public Functor0<Result>
     {
         inline virtual ~DefaultConstructor0() {}
 
         inline Result operator()() const {return Result();}
 
         inline bool operator==(const DefaultConstructor0&) const {return true;}
         inline bool operator!=(const DefaultConstructor0&) const {return false;}
     };
 
     /**
     // Simple functor which ignores is arguments and instead
     // builds the result using the default constructor of the result type
     */
     template <typename Result, typename Arg1>
     struct DefaultConstructor1 : public Functor1<Result, Arg1>
     {
         inline virtual ~DefaultConstructor1() {}
 
         inline Result operator()(const Arg1&) const {return Result();}
 
         inline bool operator==(const DefaultConstructor1&) const {return true;}
         inline bool operator!=(const DefaultConstructor1&) const {return false;}
     };
 
     /**
     // Simple functor which ignores is arguments and instead
     // builds the result using the default constructor of the result type
     */
     template <typename Result, typename Arg1, typename Arg2>
     struct DefaultConstructor2 : public Functor2<Result, Arg1, Arg2>
     {
         inline virtual ~DefaultConstructor2() {}
 
         inline Result operator()(const Arg1&, const Arg2&) const
             {return Result();}
 
         inline bool operator==(const DefaultConstructor2&) const {return true;}
         inline bool operator!=(const DefaultConstructor2&) const {return false;}
     };
 
     /**
     // Simple functor which ignores is arguments and instead
     // builds the result using the default constructor of the result type
     */
     template <typename Result, typename Arg1, typename Arg2, typename Arg3>
     struct DefaultConstructor3 : public Functor3<Result, Arg1, Arg2, Arg3>
     {
         inline virtual ~DefaultConstructor3() {}
 
         inline Result operator()(const Arg1&, const Arg2&, const Arg3&) const
             {return Result();}
 
         inline bool operator==(const DefaultConstructor3&) const {return true;}
         inline bool operator!=(const DefaultConstructor3&) const {return false;}
     };
 
     /** Simple functor which returns a constant */
     template <typename Result>
     class ConstValue0 : public Functor0<Result>
     {
     public:
         inline ConstValue0(const Result& v) : value_(v) {}
         inline virtual ~ConstValue0() {}
 
         inline Result operator()() const {return value_;}
 
         inline bool operator==(const ConstValue0& r) const
             {return value_ == r.value_;}
         inline bool operator!=(const ConstValue0& r) const
             {return !(*this == r);}
 
     private:
         ConstValue0();
         Result value_;
     };
 
     /** Simple functor which returns a constant */
     template <typename Result, typename Arg1>
     class ConstValue1 : public Functor1<Result, Arg1>
     {
     public:
         inline ConstValue1(const Result& v) : value_(v) {}
         inline virtual ~ConstValue1() {}
 
         inline Result operator()(const Arg1&) const {return value_;}
 
         inline bool operator==(const ConstValue1& r) const
             {return value_ == r.value_;}
         inline bool operator!=(const ConstValue1& r) const
             {return !(*this == r);}
 
     private:
         ConstValue1();
         Result value_;
     };
 
     /** Simple functor which returns a constant */
     template <typename Result, typename Arg1, typename Arg2>
     class ConstValue2 : public Functor2<Result, Arg1, Arg2>
     {
     public:
         inline ConstValue2(const Result& v) : value_(v) {}
         inline virtual ~ConstValue2() {}
 
         inline Result operator()(const Arg1&, const Arg2&) const
             {return value_;}
 
         inline bool operator==(const ConstValue2& r) const
             {return value_ == r.value_;}
         inline bool operator!=(const ConstValue2& r) const
             {return !(*this == r);}
 
     private:
         ConstValue2();
         Result value_;
     };
 
     /** Simple functor which returns a constant */
     template <typename Result, typename Arg1, typename Arg2, typename Arg3>
     class ConstValue3 : public Functor3<Result, Arg1, Arg2, Arg3>
     {
     public:
         inline ConstValue3(const Result& v) : value_(v) {}
         inline virtual ~ConstValue3() {}
 
         inline Result operator()(const Arg1&, const Arg2&, const Arg3&) const
             {return value_;}
 
         inline bool operator==(const ConstValue3& r) const
             {return value_ == r.value_;}
         inline bool operator!=(const ConstValue3& r) const
             {return !(*this == r);}
 
     private:
         ConstValue3();
         Result value_;
     };
 
     /**
     // Sometimes it becomes necessary to perform an explicit cast for proper
     // overload resolution of a converting copy constructor
     */
     template <typename Result, typename Arg1, typename CastType>
     struct CastingCopyConstructor : public Functor1<Result, Arg1>
     {
         inline virtual ~CastingCopyConstructor() {}
 
         inline Result operator()(const Arg1& a) const
             {return Result(static_cast<CastType>(a));}
 
         inline bool operator==(const CastingCopyConstructor&) const
             {return true;}
         inline bool operator!=(const CastingCopyConstructor&) const
             {return false;}
     };
 
     /**
     // Adaptation for using functors without arguments with simple
     // cmath-style functions
     */
     template <typename Result>
     struct FcnFunctor0 : public Functor0<Result>
     {
         inline explicit FcnFunctor0(Result (*fcn)()) : fcn_(fcn) {}
         inline virtual ~FcnFunctor0() {}
 
         inline Result operator()() const {return fcn_();}
 
         inline bool operator==(const FcnFunctor0& r) const
             {return fcn_ == r.fcn_;}
         inline bool operator!=(const FcnFunctor0& r) const
             {return !(*this == r);}
 
     private:
         FcnFunctor0();
         Result (*fcn_)();
     };
 
     /**
     // Adaptation for using single-argument functors with simple
     // cmath-style functions
     */
     template <typename Result, typename Arg1>
     struct FcnFunctor1 : public Functor1<Result, Arg1>
     {
         inline explicit FcnFunctor1(Result (*fcn)(Arg1)) : fcn_(fcn) {}
         inline virtual ~FcnFunctor1() {}
 
         inline Result operator()(const Arg1& a) const {return fcn_(a);}
 
         inline bool operator==(const FcnFunctor1& r) const
             {return fcn_ == r.fcn_;}
         inline bool operator!=(const FcnFunctor1& r) const
             {return !(*this == r);}
 
     private:
         FcnFunctor1();
         Result (*fcn_)(Arg1);
     };
 
     /**
     // Adaptation for using two-argument functors with simple
     // cmath-style functions
     */
     template <typename Result, typename Arg1, typename Arg2>
     struct FcnFunctor2 : public Functor2<Result, Arg1, Arg2>
     {
         inline explicit FcnFunctor2(Result (*fcn)(Arg1, Arg2)) : fcn_(fcn) {}
         inline virtual ~FcnFunctor2() {}
 
         inline Result operator()(const Arg1& x, const Arg2& y) const
             {return fcn_(x, y);}
 
         inline bool operator==(const FcnFunctor2& r) const
             {return fcn_ == r.fcn_;}
         inline bool operator!=(const FcnFunctor2& r) const
             {return !(*this == r);}
 
     private:
         FcnFunctor2();
         Result (*fcn_)(Arg1, Arg2);
     };
 
     /**
     // Adaptation for using three-argument functors with simple
     // cmath-style functions
     */
     template <typename Result, typename Arg1, typename Arg2, typename Arg3>
     struct FcnFunctor3 : public Functor3<Result, Arg1, Arg2, Arg3>
     {
         inline explicit FcnFunctor3(Result (*fcn)(Arg1,Arg2,Arg3)):fcn_(fcn) {}
         inline virtual ~FcnFunctor3() {}
 
         inline Result operator()(const Arg1&x,const Arg2&y,const Arg3&z) const
             {return fcn_(x, y, z);}
 
         inline bool operator==(const FcnFunctor3& r) const
             {return fcn_ == r.fcn_;}
         inline bool operator!=(const FcnFunctor3& r) const
             {return !(*this == r);}
 
     private:
         FcnFunctor3();
         Result (*fcn_)(Arg1, Arg2, Arg3);
     };
 
     template <
         typename Result1, typename Arg1,
         typename Result2, typename Arg2,
         class Operator
     >
     struct ComboFunctor1Helper : public Functor1<Result1, Arg1>
     {
         inline ComboFunctor1Helper(const Functor1<Result1, Arg1>& f1,
                                    const Functor1<Result2, Arg2>& f2,
                                    const Operator& op)
             : f1_(f1), f2_(f2), op_(op) {}
 
         inline virtual ~ComboFunctor1Helper() {}
 
         inline virtual Result1 operator()(const Arg1& x) const
         {
             return op_(f1_(x), static_cast<Result1>(f2_(static_cast<Arg2>(x))));
         }
 
     private:
         const Functor1<Result1, Arg1>& f1_;
         const Functor1<Result2, Arg2>& f2_;
         const Operator op_;
     };
 
     /**
     // Create a functor which operates on the results produced
     // by two other functors. Note that only references to f1
     // and f2 are stored. Lifetimes of f1 and f2 must be longer
     // than the lifetime of this functor. The combo functor will
     // have the result type and the argument type of f1.
     */
     template <
         typename Result1, typename Arg1,
         typename Result2, typename Arg2,
         template<typename> class Operator
     >
     inline ComboFunctor1Helper<Result1, Arg1, Result2, Arg2, Operator<Result1> >
     ComboFunctor1(const Functor1<Result1, Arg1>& f1,
                   const Functor1<Result2, Arg2>& f2,
                   const Operator<Result1>& op)
     {
         return ComboFunctor1Helper<Result1, Arg1, Result2, Arg2, Operator<Result1> >(f1, f2, op);
     }
 
     /**
     // Functor which extracts a given element from a random access linear
     // container without bounds checking
     */
     template <class Container, class Result = typename Container::value_type>
     struct Element1D : public Functor1<Result, Container>
     {
         inline explicit Element1D(const unsigned long index) : idx(index) {}
         inline virtual ~Element1D() {}
 
         inline Result operator()(const Container& c) const {return c[idx];}
 
         inline bool operator==(const Element1D& r) const
             {return idx == r.idx;}
         inline bool operator!=(const Element1D& r) const
             {return !(*this == r);}
 
     private:
         Element1D();
         unsigned long idx;
     };
 
     /**
     // Functor which extracts a given element from a random access linear
     // container with bounds checking
     */
     template <class Container, class Result = typename Container::value_type>
     struct Element1DAt : public Functor1<Result, Container>
     {
         inline explicit Element1DAt(const unsigned long index) : idx(index) {}
         inline virtual ~Element1DAt() {}
 
         inline Result operator()(const Container& c) const {return c.at(idx);}
 
         inline bool operator==(const Element1DAt& r) const
             {return idx == r.idx;}
         inline bool operator!=(const Element1DAt& r) const
             {return !(*this == r);}
 
     private:
         Element1DAt();
         unsigned long idx;
     };
 
     /** 
     // Left assignment functor. Works just like normal binary
     // assignment operator in places where functor is needed.
     */
     template <typename T1, typename T2>
     struct assign_left
     {
         inline T1& operator()(T1& left, const T2& right) const
             {return left = right;}
     };
 
     /** 
     // Right assignment functor. Reverses the assignment direction
     // in comparison with the normal binary assignment operator.
     */
     template <typename T1, typename T2>
     struct assign_right
     {
         inline T2& operator()(const T1& left, T2& right) const
             {return right = left;}
     };
 
     /** In-place addition on the left side */
     template <typename T1, typename T2>
     struct pluseq_left
     {
         inline T1& operator()(T1& left, const T2& right) const
             {return left += right;}
     };
 
     /** In-place addition on the right side */
     template <typename T1, typename T2>
     struct pluseq_right
     {
         inline T2& operator()(const T1& left, T2& right) const
             {return right += left;}
     };
 
     /** 
     // In-place addition on the left side preceded by
     // multiplication of the right argument by a double
     */
     template <typename T1, typename T2>
     struct addmul_left
     {
         inline explicit addmul_left(const double weight) : w_(weight) {}
 
         inline T1& operator()(T1& left, const T2& right) const
             {return left += w_*right;}
 
     private:
         addmul_left();
         double w_;
     };
 
     /** 
     // In-place addition on the right side preceded by
     // multiplication of the left argument by a double
     */
     template <typename T1, typename T2>
     struct addmul_right
     {
         inline explicit addmul_right(const double weight) : w_(weight) {}
 
         inline T1& operator()(T1& left, const T2& right) const
             {return right += w_*left;}
 
     private:
         addmul_right();
         double w_;
     };
 
     /** In-place subtraction on the left side */
     template <typename T1, typename T2>
     struct minuseq_left
     {
         inline T1& operator()(T1& left, const T2& right) const
             {return left -= right;}
     };
 
     /** In-place subtraction on the right side */
     template <typename T1, typename T2>
     struct minuseq_right
     {
         inline T2& operator()(const T1& left, T2& right) const
             {return right -= left;}
     };
 
     /** In-place multiplication on the left side */
     template <typename T1, typename T2>
     struct multeq_left
     {
         inline T1& operator()(T1& left, const T2& right) const
             {return left *= right;}
     };
 
     /** In-place multiplication on the right side */
     template <typename T1, typename T2>
     struct multeq_right
     {
         inline T2& operator()(const T1& left, T2& right) const
             {return right *= left;}
     };
 
     /** In-place division on the left side without checking for division by 0 */
     template <typename T1, typename T2>
     struct diveq_left
     {
         inline T1& operator()(T1& left, const T2& right) const
             {return left /= right;}
     };
 
     /** In-place division on the right side without checking for division by 0 */
     template <typename T1, typename T2>
     struct diveq_right
     {
         inline T2& operator()(const T1& left, T2& right) const
             {return right /= left;}
     };
 
     /** In-place division on the left side. Allow 0/0 = const. */
     template <typename T1, typename T2>
     struct diveq_left_0by0isC
     {
         inline diveq_left_0by0isC() : 
             C(T1()), leftZero(T1()), rightZero(T2()) {}
         inline explicit diveq_left_0by0isC(const T1& value) :
             C(value), leftZero(T1()), rightZero(T2()) {}
 
         inline T1& operator()(T1& left, const T2& right) const
         {
             if (right == rightZero)
                 if (left == leftZero)
                 {
                     left = C;
                     return left;
                 }
             return left /= right;
         }
 
     private:
         T1 C;
         T1 leftZero;
         T2 rightZero;
     };
 
     /** In-place division on the right side. Allow 0/0 = const. */
     template <typename T1, typename T2>
     struct diveq_right_0by0isC
     {
         inline diveq_right_0by0isC() :
             C(T2()), leftZero(T1()), rightZero(T2())  {}
         inline explicit diveq_right_0by0isC(const T2& value) :
             C(value), leftZero(T1()), rightZero(T2()) {}
 
         inline T2& operator()(const T1& left, T2& right) const
         {
             if (left == leftZero)
                 if (right == rightZero)
                 {
                     right = C;
                     return right;
                 }
             return right /= left;
         }
 
     private:
         T2 C;
         T1 leftZero;
         T2 rightZero;
     };
 
     /** Left assignment functor preceded by a static cast */
     template <typename T1, typename T2, typename T3=T1>
     struct scast_assign_left
     {
         inline T1& operator()(T1& left, const T2& right) const
             {return left = static_cast<T3>(right);}
     };
 
     /** Right assignment functor preceded by a static cast */
     template <typename T1, typename T2, typename T3=T2>
     struct scast_assign_right
     {
         inline T2& operator()(const T1& left, T2& right) const
             {return right = static_cast<T3>(left);}
     };
 
     /** In-place addition on the left side preceded by a static cast */
     template <typename T1, typename T2, typename T3=T1>
     struct scast_pluseq_left
     {
         inline T1& operator()(T1& left, const T2& right) const
             {return left += static_cast<T3>(right);}
     };
 
     /** In-place addition on the right side preceded by a static cast */
     template <typename T1, typename T2, typename T3=T2>
     struct scast_pluseq_right
     {
         inline T2& operator()(const T1& left, T2& right) const
             {return right += static_cast<T3>(left);}
     };
 
     /** In-place subtraction on the left side preceded by a static cast */
     template <typename T1, typename T2, typename T3=T1>
     struct scast_minuseq_left
     {
         inline T1& operator()(T1& left, const T2& right) const
             {return left -= static_cast<T3>(right);}
     };
 
     /** In-place subtraction on the right side preceded by a static cast */
     template <typename T1, typename T2, typename T3=T2>
     struct scast_minuseq_right
     {
         inline T2& operator()(const T1& left, T2& right) const
             {return right -= static_cast<T3>(left);}
     };
 
     /** Square root functor */
     template <typename T1>
     struct SquareRoot : public Functor1<T1, T1>
     {
         inline virtual ~SquareRoot() {}
 
         inline virtual T1 operator()(const T1& x) const
         {
             static const T1 zero = static_cast<T1>(0);
 
             if (!(x >= zero)) throw std::domain_error(
                 "In npstat::SquareRoot::operator(): "
                 "argument must be non-negative");
             return std::sqrt(x);
         }
     };
 
     /** Inverse square root functor */
     template <typename T1>
     struct InverseSquareRoot : public Functor1<T1, T1>
     {
         inline virtual ~InverseSquareRoot() {}
 
         inline virtual T1 operator()(const T1& x) const
         {
             static const T1 zero = static_cast<T1>(0);
             static const T1 one = static_cast<T1>(1);
 
             if (!(x > zero)) throw std::domain_error(
                 "In npstat::InverseSquareRoot::operator(): "
                 "argument must be positive");
             return one/std::sqrt(x);
         }
     };
 
+    /** Useful functor for calculating ISE between two densities */
+    template <class F1, class F2>
+    class DeltaSquaredFcnHelper : public Functor1<double, double>
+    {
+    public:
+        inline DeltaSquaredFcnHelper(const F1& f1, const F2& f2)
+            : f1_(f1), f2_(f2) {}
+
+        inline virtual ~DeltaSquaredFcnHelper() {}
+
+        inline virtual double operator()(const double& x) const
+        {
+            const double delta = f1_(x) - f2_(x);
+            return delta*delta;
+        }
+
+    private:
+        const F1& f1_;
+        const F2& f2_;
+    };
+
+    template <class F1, class F2>
+    inline DeltaSquaredFcnHelper<F1,F2>
+    DeltaSquaredFcn(const F1& f1, const F2& f2)
+    {
+        return DeltaSquaredFcnHelper<F1,F2>(f1, f2);
+    }
+
     /**
     // Useful functor for calculating ISE between a function and
     // a polynomial series
     */
     template <class Poly, class Functor>
     class DeltaSquaredSerFcn : public Functor1<double, double>
     {
     public:
         // Note that the series coefficients are not copied
         inline DeltaSquaredSerFcn(
             const Poly& poly, const Functor& fcn,
             const double* coeffs, const unsigned maxdeg)
             : poly_(poly), fcn_(fcn), coeffs_(coeffs), maxdeg_(maxdeg) {}
 
         inline virtual ~DeltaSquaredSerFcn() {}
 
-        inline double operator()(const double& x) const
+        inline virtual double operator()(const double& x) const
         {
             const double delta = poly_.series(coeffs_, maxdeg_, x) - fcn_(x);
             return delta*delta;
         }
 
     private:
         const Poly& poly_;
         const Functor& fcn_;
         const double* coeffs_;
         unsigned maxdeg_;
     };
 
     /**
     // Useful functor for calculating powers of another functor or
     // a cmath-style function. Does not make a copy of the argument
     // functor, using only a reference.
     */
     template <class Functor>
     class FunctorPowerFcnHelper : public Functor1<double, double>
     {
     public:
         inline FunctorPowerFcnHelper(const Functor& fcn, const double power)
             : fcn_(fcn), power_(power) {}
 
         inline virtual ~FunctorPowerFcnHelper() {}
         
-        inline double operator()(const double& x) const
+        inline virtual double operator()(const double& x) const
         {
             if (power_)
             {
                 const double tmp = fcn_(x);
                 if (power_ == 1.0)
                     return tmp;
                 else if (power_ == 2.0)
                     return tmp*tmp;
                 else if (power_ == -1.0)
                 {
                     assert(tmp);
                     return 1.0/tmp;
                 }
                 else
                 {
                     if (power_ < 0.0)
                         assert(tmp);
                     return std::pow(tmp, power_);
                 }
             }
             else
                 return 1.0;
         }
 
     private:
         const Functor& fcn_;
         double power_;
     };
 
     /** Utility function for making FunctorPowerFcnHelper objects */
     template <class Functor>
     inline FunctorPowerFcnHelper<Functor> FunctorPowerFcn(
         const Functor& fcn, const double power)
     {
         return FunctorPowerFcnHelper<Functor>(fcn, power);
     }
 
     /**
     // Useful functor for calculating powers of another functor or
     // a cmath-style function. Makes a copy of the argument functor.
     */
     template <class Functor>
     class FunctorPowerFcnCpHelper : public Functor1<double, double>
     {
     public:
         inline FunctorPowerFcnCpHelper(const Functor& fcn, const double power)
             : fcn_(fcn), power_(power) {}
 
         inline virtual ~FunctorPowerFcnCpHelper() {}
 
-        inline double operator()(const double& x) const
+        inline virtual double operator()(const double& x) const
             {return FunctorPowerFcnHelper<Functor>(fcn_, power_)(x);}
 
     private:
         Functor fcn_;
         double power_;
     };
 
     /** Utility function for making FunctorPowerFcnCpHelper objects */
     template <class Functor>
     inline FunctorPowerFcnCpHelper<Functor> FunctorPowerFcnCp(
         const Functor& fcn, const double power)
     {
         return FunctorPowerFcnCpHelper<Functor>(fcn, power);
     }
 
     /** Functor for multiplying the result of another functor by a constant */
     template <typename Result, typename Arg1>
     class MultiplyByConstHelper : public Functor1<Result, Arg1>
     {
     public:
         inline MultiplyByConstHelper(const Functor1<Result, Arg1>& fcn,
                                      const Result& factor)
             : fcn_(fcn), factor_(factor) {}
 
         inline virtual ~MultiplyByConstHelper() {}
 
-        inline double operator()(const Arg1& x) const
+        inline virtual double operator()(const Arg1& x) const
             {return factor_*fcn_(x);}
 
     private:
         const Functor1<Result, Arg1>& fcn_;
         Result factor_;
     };
 
     /** Utility function for making MultiplyByConstHelper objects */
     template <typename Result, typename Arg1, typename Numeric>
     inline MultiplyByConstHelper<Result, Arg1> MultiplyByConst(
         const Functor1<Result, Arg1>& fcn, const Numeric& factor)
     {
         return MultiplyByConstHelper<Result, Arg1>(fcn, static_cast<Result>(factor));
     }
 }
 
 #endif // NPSTAT_SIMPLEFUNCTORS_HH_