Program Listing for File PDFIntegrator.hpp

↰ Return to documentation for file (src/utils/PDFIntegrator.hpp)

#ifndef PDFINTEGRATOR_HPP
#define PDFINTEGRATOR_HPP

#include "../pdfs/PDF.hpp"
#include "GKIntegrator.hpp"
#include "GKIntegratorData.hpp"
#include "MiscUtils.hpp"
#include <algorithm>
#include <cmath>
#include <cstddef>
#include <functional>
#include <ranges>
#include <stdexcept>
#include <string>
#include <type_traits>
#include <utility>
#include <vector>

namespace utils
{

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

        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
        ) :
            p_(p),
            f_(std::move(f)),
            nInt_(nInt),
            logTransform_(logTransform),
            maxIter_(maxIter),
            absTol_(absTol),
            relTol_(relTol)
        { }

        // Disallow copying and moving: p_ holds a reference, so a copy
        // would alias the same PDF as the original rather than being
        // an independent instance, and there is no real need to copy
        // or move a PDFIntegrator in practice, since it is meant to be
        // constructed and used immediately
        PDFIntegrator(const PDFIntegrator&) = delete;
        PDFIntegrator(PDFIntegrator&&) = delete;
        auto operator=(const PDFIntegrator&) -> PDFIntegrator& = delete;
        auto operator=(PDFIntegrator&&) -> PDFIntegrator& = delete;
        ~PDFIntegrator() = default;

        template <class... Args>
        [[nodiscard]] auto integrand(const double x, Args&... args) const -> std::vector<double>
        {
            const double xReal = logTransform_ ? std::exp(x) : x;
            const double xClamped = std::clamp(xReal, aReal_, bReal_);

            const auto val = invokeF(xClamped, args...);

            double weight = p_(xClamped);
            if (logTransform_) { weight *= xClamped; }

            std::vector<double> result(std::begin(val), std::end(val));
            for (double& v : result) { v *= weight; }
            return result;
        }

        template <class... Args>
        [[nodiscard]] auto integrate(double a, double b, Args&&... args) const -> std::vector<double>
        {
            a = std::max(a, p_.getMin());
            b = std::min(b, p_.getMax());
            aReal_ = a;
            bReal_ = b;

            // An empty (or zero-width) clamped interval integrates to 0,
            // even with a delta-function segment at its one point,
            // matching PDF::integral(a, b)
            if (a >= b)
            {
                std::vector<double> zero(nInt_, 0.0); // not a braced return: {nInt_, 0.0} would be a two-element list
                return zero;
            }

            if (logTransform_)
            {
                if (a <= 0.0 || b <= 0.0)
                {
                    throw std::runtime_error(
                        "PDFIntegrator::integrate: logTransform is set, but the "
                        "integration bounds (after clamping to the PDF's own support) "
                        "are not both strictly positive: a = " + std::to_string(a) +
                        ", b = " + std::to_string(b));
                }
                a = std::log(a);
                b = std::log(b);
            }

            // Delta-function contributions, evaluated before args are
            // forwarded (and possibly moved from) below
            std::vector<double> deltaSum(nInt_, 0.0);
            for (const auto& [x0, w] : p_.deltas())
            {
                if (x0 < aReal_ || x0 > bReal_) { continue; }
                const auto val = invokeF(x0, args...);
                for (auto&& [d, v] : std::views::zip(deltaSum, val)) { d += w * v; }
            }

            using IntegrandFn = decltype(&PDFIntegrator::template integrand<Args...>);
            const GKIntegrator<IntegrandFn, Order> integ(
                &PDFIntegrator::template integrand<Args...>, nInt_, maxIter_, absTol_, relTol_);
            auto result = integ.integrate(a, b, this, std::forward<Args>(args)...);
            for (auto&& [r, d] : std::views::zip(result, deltaSum)) { r += d; }
            return result;
        }

    private:

        template <class... Args>
        [[nodiscard]] auto invokeF(const double x, Args&&... args) const
        {
            if constexpr (std::is_member_function_pointer_v<F>)
            {
                return invokeFMember(x, std::forward<Args>(args)...);
            }
            else
            {
                return std::invoke(f_, x, std::forward<Args>(args)...);
            }
        }

        template <class Obj, class... Rest>
        [[nodiscard]] auto invokeFMember(const double x, Obj&& obj, Rest&&... rest) const
        {
            return std::invoke(f_, std::forward<Obj>(obj), x, std::forward<Rest>(rest)...);
        }

        const pdfs::PDF& p_;
        F f_;
        std::size_t nInt_;
        bool logTransform_;
        std::size_t maxIter_;
        double absTol_;
        double relTol_;
        // The current integrate() call's own real-space integration
        // bounds (after clamping to p_'s own support) -- set at the
        // top of integrate(), read by integrand() (see its own
        // comment for why) -- mutable so integrate() itself can stay
        // const. See this class's own @details for the concurrency
        // trade-off this makes.
        mutable double aReal_ = 0.0;
        mutable double bReal_ = 0.0;
    };

    template <class F, GKOrder Order = GKOrder::GK15>
    auto integrateScaled(const pdfs::PDF& p, const F& f, const std::size_t nInt,
        const double a, const double b, const std::vector<double>& scalePoints,
        const std::size_t maxIter, const double absTol, const double relTol,
        const bool logTransform = false) -> std::vector<double>
    {
        const double scale = integrandScale(f, scalePoints);
        const auto scaled = [&f, scale](const double x) -> std::vector<double>
        {
            auto v = f(x);
            for (double& e : v) { e /= scale; }
            return v;
        };
        const PDFIntegrator<decltype(scaled), Order> integrator(
            p, scaled, nInt, logTransform, maxIter, absTol * p.integral(a, b), relTol);
        auto result = integrator.integrate(a, b);
        for (double& e : result) { e *= scale; }
        return result;
    }

} // namespace utils

#endif // PDFINTEGRATOR_HPP