.. _program_listing_file_src_utils_PDFIntegrator.hpp:

Program Listing for File PDFIntegrator.hpp
==========================================

|exhale_lsh| :ref:`Return to documentation for file <file_src_utils_PDFIntegrator.hpp>` (``src/utils/PDFIntegrator.hpp``)

.. |exhale_lsh| unicode:: U+021B0 .. UPWARDS ARROW WITH TIP LEFTWARDS

.. code-block:: cpp

   
   #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