Template Class PDFIntegrator

• Defined in File PDFIntegrator.hpp

Class Documentation

template<class F, GKOrder Order = GKOrder::GK15>
class PDFIntegrator

Evaluate integrals of the form int_a^b p(x) f(x) dx, using GKIntegrator.

Evaluates int_a^b p(x) f(x) dx by adaptively bisecting [a, b] with GKIntegrator, weighting f by p at every point it visits – see integrand()’s own comment for exactly how the two are combined. Benchmarked (as PDFIntegratorGK, its own name before replacing an earlier, cubature-package-based implementation of this same class) head to head against that implementation on a real, full-scale run: consistently faster, in one representative case by roughly an order of magnitude – see this class’s own git history for the comparison.

integrate() is not safe to call concurrently on a single, shared PDFIntegrator instance: integrand() (see its own comment) reads the current call’s own real-space integration bounds back from two mutable members integrate() sets just before it, rather than from a per-call context object. Every call site in this codebase constructs a fresh PDFIntegrator immediately before using it, so this has not been a practical limitation; revisit if a future caller ever needs to share one PDFIntegrator across concurrent integrate() calls.

Template Parameters:

• F – Type of the callable integrand: either a standalone callable, invocable as f(x, …), where x is the (scalar) integration variable and … is any number of additional arguments forwarded from integrate(); or a pointer to a member function, in which case the first of those additional arguments must instead be an instance of the class f is a member of (x is still passed as that member function’s own first ordinary argument, immediately after the instance). Either way, f must return a contiguous container of nInt_ doubles (e.g. std::vector<double> or std::array<double, nInt_>) – the value of the integrand at x.

• Order – Which Gauss-Kronrod pair GKIntegrator integrates each subinterval with – see GKOrder’s own comment. Defaults to GKOrder::GK15.

Public Functions

inline PDFIntegrator(const pdfs::PDF &p, F f, std::size_t nInt, bool logTransform = false, std::size_t maxIter = 0, double absTol = 0.0, double relTol = 1e-6)

Construct a PDFIntegrator.

Parameters:

• p – The PDF p(x) weighting the integrand; must outlive this PDFIntegrator, since it is stored by reference

• f – The integrand f(x); see this class’s own

• nInt – The number of quantities f returns per point

• logTransform – Whether to integrate over ln(x) instead of x directly – a mathematically exact change of variables (int p(x) f(x) dx = int p(e^u) f(e^u) e^u du, u = ln(x)) that leaves the value of the integral, f’s own required signature, and the coordinate f/p are evaluated at (always the real, untransformed x, regardless of this parameter) all unchanged – it only changes which variable GKIntegrator itself subdivides. Defaults to false (no transform). integrate() throws std::runtime_error if this is true but its own integration bounds (after clamping to p’s own support) are not both strictly positive, since ln() of a non-positive bound is undefined.

• maxIter – See GKIntegrator’s own maxIter parameter

• absTol – See GKIntegrator’s own absTol parameter

• relTol – See GKIntegrator’s own relTol parameter

Template Parameters:

F –

PDFIntegrator(const PDFIntegrator&) = delete

PDFIntegrator(PDFIntegrator&&) = delete

auto operator=(const PDFIntegrator&) -> PDFIntegrator& = delete

auto operator=(PDFIntegrator&&) -> PDFIntegrator& = delete

~PDFIntegrator() = default

template<class ...Args>
inline auto integrand(const double x, Args&... args) const -> std::vector<double>

Evaluate p(x) f(x) at a single point, for GKIntegrator’s own benefit.

First undoes the log transform, if any (exp(x) -> xReal), then clamps xReal into [aReal_, bReal_] – the current integrate() call’s own real-space integration bounds (see their own comment for why these, not p_’s own full getMin()/getMax(), are the correct clamp target: an isochrone-segment-scale sub-range of p_’s own domain is the common case in this codebase, and clamping to p_’s own full domain instead could let a point land outside such a segment’s own valid range). This clamp matters because exp() and ln() are not exact inverses under floating point, so exp(x) can land fractionally outside [aReal_, bReal_] even when x itself was exactly ln(aReal_)/ln(bReal_)’s own (correctly clamped) value.

Evaluates f_ at xReal (via invokeF(), which transparently supports F being a pointer to member function – see its own comment), multiplies elementwise by p_(xReal), and, if logTransform_, by xReal again (the Jacobian).

Args&… here, not the more usual Args&&… forwarding reference: integrate() only ever names this function through &PDFIntegrator::template integrand<Args…>, with Args already fixed to its own deduced types – at that point Args&& is no longer a forwarding reference at all, just an ordinary (possibly rvalue) reference type, and for any Args element integrate() deduced as a plain value (e.g. a pointer instance passed as a prvalue, as PDFIntegrator’s own integrate() itself does when it is F for an outer GKIntegrator/PDFIntegrator – see specsyn::Specsyn:: specCtsHelper() for a real example, nesting one PDFIntegrator inside another this way), that becomes a genuine rvalue-reference parameter. But GKIntegrator’s own internal machinery (see its quadSingle()/integrate() own comments on argsTuple) always hands every call an lvalue, to stay safe across the many repeated calls one adaptive integral makes – so a plain lvalue reference is the only parameter type guaranteed to bind correctly here, regardless of what value category the original argument at integrate()’s own call site had.

Parameters:

• x – The integration variable – if logTransform_ is true, this is ln of the real coordinate (the domain integrate() itself hands GKIntegrator), not the real coordinate directly

• args – Any additional arguments f_ requires; if F is a pointer to member function, the first of these must be an instance of the class it is a member of (see this class’s own

Template Parameters:

F) –

Returns:

p(xReal) * f(xReal, args…), elementwise, further scaled by xReal (the change-of-variables Jacobian dx = x d(ln x)) if logTransform_ is true, where xReal is x itself (logTransform_ false) or exp(x) (logTransform_ true), clamped into [aReal_, bReal_]

template<class ...Args>
inline auto integrate(double a, double b, Args&&... args) const -> std::vector<double>

Evaluate the integral of p(x) f(x) from a to b.

Clamps [a, b] to p_’s own support, records the result in aReal_/bReal_ for integrand()’s own benefit (see its own comment for why), and returns all zeros if the clamped interval is empty or has zero width (as PDF::integral(a, b) does); otherwise log-transforms a/b themselves if logTransform_ is set, then constructs a GKIntegrator whose own integrand is this class’s integrand() (a pointer to member function, so this is passed as the first extra argument to GKIntegrator::integrate() below – see GKIntegrator’s own

Parameters:

• a – Lower limit of integration

• b – Upper limit of integration

• args – Any additional arguments f_ requires, exactly as they would be passed to integrand()

Throws:

std::runtime_error – if logTransform_ is true, but a/b (after clamping to p_’s own support) are not both strictly positive

Template Parameters:

F – comment), passing maxIter_/ absTol_/relTol_ through unchanged. Finally adds the contribution of each of p_’s own delta-function segments (see pdfs::PDF::deltas()) lying within [a, b], inclusive: p_’s density excludes them, but each contributes its own weight times f at its own location, since the integral of f(x) delta(x - x0) is f(x0).

Returns:

The integral, one value per quantity f_ returns (length nInt_)