Page MenuHomeHEPForge
Files F9501737

No OneTemporary
Actions

View Options

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_

File Metadata

Mime Type
text/x-diff
Expires
Sun, Feb 23, 3:01 PM (21 m, 58 s)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4486794
Default Alt Text
(46 KB)

Event Timeline