.. _program_listing_file_src_utils_GKIntegrator.hpp:

Program Listing for File GKIntegrator.hpp
=========================================

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

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

.. code-block:: cpp

   
   #ifndef GKINTEGRATOR_HPP
   #define GKINTEGRATOR_HPP
   
   #include "GKIntegratorData.hpp"
   #include <algorithm>
   #include <array>
   #include <cmath>
   #include <cstddef>
   #include <functional>
   #include <tuple>
   #include <type_traits>
   #include <utility>
   #include <vector>
   
   namespace utils
   {
   
       struct GKIntegrationInterval
       {
           double a_;                  
           double b_;                  
           std::vector<double> quad_;  
           std::vector<double> err_;   
       };
   
       template <class F, GKOrder Order>
       class GKIntegrator
       {
       public:
   
           GKIntegrator(
               F f,
               std::size_t nInt,
               std::size_t maxIter = 0,
               double absTol = 0.0,
               double relTol = 1e-6
           ) :
               f_(std::move(f)),
               nInt_(nInt),
               maxIter_(maxIter),
               absTol_(absTol),
               relTol_(relTol)
           { }
   
           template <class... Args>
           [[nodiscard]] auto quadSingle(const double a, const double b, Args&&... args) const
           -> GKIntegrationInterval
           {
               using Rule = GKRule<Order>;
               constexpr std::size_t ngk = Rule::ngk;
               constexpr std::size_t gknum = (2 * ngk) - 1;
   
               const double center = 0.5 * (a + b);
               const double halfLength = 0.5 * (b - a);
   
               std::array<double, gknum> xk{};
               for (std::size_t i = 0; i < gknum / 2; ++i)
               {
                   xk[i] = center - (halfLength * Rule::xgk[i]); // NOLINT(cppcoreguidelines-pro-bounds-constant-array-index) -- i < gknum / 2 == ngk - 1 < ngk == Rule::xgk.size() by construction
                   xk[gknum - i - 1] = center + (halfLength * Rule::xgk[i]); // NOLINT(cppcoreguidelines-pro-bounds-constant-array-index) -- gknum - i - 1 >= gknum - gknum/2 > 0 and < gknum, and i < ngk as above
               }
               xk[gknum / 2] = center; // NOLINT(cppcoreguidelines-pro-bounds-constant-array-index) -- gknum / 2 < gknum by construction
   
               // Evaluate f_ at every point first, rather than
               // interleaving evaluation with the summation loops below,
               // to mirror integrate_gk's own structure -- each point's
               // own nInt_ values are copied out of whatever contiguous
               // container f_ returns (see this class's own @tparam F)
               // into a plain std::vector<double>, since nInt_ itself is
               // only known at runtime. args is forwarded exactly once,
               // into this tuple, rather than at each of the gknum calls
               // below: repeatedly forwarding the same forwarding-reference
               // parameter would risk moving from an rvalue-bound argument
               // on its first use, leaving it in a moved-from state for
               // every later point -- std::apply instead hands every call
               // the tuple's own elements, exactly as deduced (not const,
               // so that invokeF()/invokeFMember() can still bind them to
               // a non-const reference parameter when F's own signature --
               // e.g. a member-function-pointer instantiation bound to a
               // specific, concrete Args... -- demands one).
               std::tuple<Args...> argsTuple(std::forward<Args>(args)...);
               std::vector<std::vector<double>> fVal(gknum);
               for (std::size_t i = 0; i < gknum; ++i)
               {
                   const auto val = std::apply(
                       [this, &xk, i](auto&... unpackedArgs) { return invokeF(xk[i], unpackedArgs...); }, // NOLINT(cppcoreguidelines-pro-bounds-constant-array-index) -- i < gknum == xk.size() by construction
                       argsTuple);
                   fVal[i].assign(std::begin(val), std::end(val));
               }
   
               std::vector<double> result(nInt_, 0.0);
               std::vector<double> gaussQuad(nInt_, 0.0);
   
               // Central point: always part of the Kronrod sum (weight
               // wgk[ngk - 1], the last stored weight -- see GKRule's own
               // comment for why); part of the Gauss sum too only when
               // ngk is even (see this function's own comment for why).
               const std::size_t centerIdx = gknum / 2;
               for (std::size_t k = 0; k < nInt_; ++k)
               {
                   result[k] = fVal[centerIdx][k] * Rule::wgk[ngk - 1]; // NOLINT(cppcoreguidelines-pro-bounds-constant-array-index) -- ngk - 1 < ngk == Rule::wgk.size() by construction
               }
               if constexpr (ngk % 2 == 0)
               {
                   for (std::size_t k = 0; k < nInt_; ++k)
                   {
                       gaussQuad[k] = fVal[centerIdx][k] * Rule::wg[(ngk / 2) - 1]; // NOLINT(cppcoreguidelines-pro-bounds-constant-array-index) -- ngk / 2 - 1 < ng == Rule::wg.size() whenever ngk is even, by the Gauss/Kronrod pairing GKRule's own comment describes
                   }
               }
   
               // Points shared between the Gauss and Kronrod sums
               for (std::size_t i = 0; i < (ngk - 1) / 2; ++i)
               {
                   const std::size_t p1 = (2 * i) + 1;
                   const std::size_t p2 = gknum - (2 * i) - 2;
                   for (std::size_t k = 0; k < nInt_; ++k)
                   {
                       const double sum = fVal[p1][k] + fVal[p2][k]; // NOLINT(cppcoreguidelines-pro-bounds-constant-array-index) -- p1, p2 < gknum == fVal.size() by construction
                       gaussQuad[k] += Rule::wg[i] * sum; // NOLINT(cppcoreguidelines-pro-bounds-constant-array-index) -- i < (ngk - 1) / 2 <= ng == Rule::wg.size(), by the Gauss/Kronrod pairing GKRule's own comment describes
                       result[k] += Rule::wgk[p1] * sum; // NOLINT(cppcoreguidelines-pro-bounds-constant-array-index) -- p1 < ngk == Rule::wgk.size() by construction
                   }
               }
   
               // Points that appear only in the Kronrod sum
               for (std::size_t i = 0; i < ngk / 2; ++i)
               {
                   const std::size_t p1 = 2 * i;
                   const std::size_t p2 = gknum - (2 * i) - 1;
                   for (std::size_t k = 0; k < nInt_; ++k)
                   {
                       result[k] += Rule::wgk[p1] * (fVal[p1][k] + fVal[p2][k]); // NOLINT(cppcoreguidelines-pro-bounds-constant-array-index) -- p1 < ngk == Rule::wgk.size(), and p1, p2 < gknum == fVal.size(), by construction
                   }
               }
   
               std::vector<double> err(nInt_);
               for (std::size_t k = 0; k < nInt_; ++k)
               {
                   result[k] *= halfLength;
                   gaussQuad[k] *= halfLength;
                   err[k] = std::abs(result[k] - gaussQuad[k]);
               }
               return { a, b, std::move(result), std::move(err) };
           }
   
           template <class... Args>
           [[nodiscard]] auto integrate(const double a, const double b, Args&&... args) const -> std::vector<double>
           {
               // A zero-width interval's own quad_/err_ both come out
               // exactly 0.0 regardless of f_ (every Gauss-Kronrod
               // abscissa collapses to the same point a == b, and quad_/
               // err_ are both scaled by the interval's own zero half-
               // width) -- but that makes the relative-error check below
               // divide 0.0 by 0.0 (a NaN that adaptive bisection would
               // otherwise need to somehow terminate on), so short-circuit
               // here instead of relying on that division's own behavior.
               if (a == b) { return std::vector<double>(nInt_, 0.0); } // NOLINT(modernize-return-braced-init-list) -- a braced return ({nInt_, 0.0}) would call vector's initializer_list constructor instead of its (count, value) one, silently returning the 2-element vector {double(nInt_), 0.0} instead of nInt_ copies of 0.0
   
               // args is forwarded exactly once, into this tuple, for the
               // same reason quadSingle() itself does this -- see its own
               // comment (including for why this is not const).
               std::tuple<Args...> argsTuple(std::forward<Args>(args)...);
               const auto callQuadSingle = [this, &argsTuple](const double lo, const double hi) -> GKIntegrationInterval
               {
                   return std::apply(
                       [this, lo, hi](auto&... unpackedArgs) { return quadSingle(lo, hi, unpackedArgs...); },
                       argsTuple);
               };
   
               std::vector<GKIntegrationInterval> intervals;
               intervals.push_back(callQuadSingle(a, b));
   
               const double absErr = *std::ranges::max_element(intervals[0].err_);
               if (absErr < absTol_) { return intervals[0].quad_; }
               {
                   double maxRelErr = 0.0;
                   for (std::size_t k = 0; k < nInt_; ++k)
                   {
                       maxRelErr = std::max(maxRelErr,
                           std::abs(intervals[0].err_[k] / intervals[0].quad_[k])); // NOLINT(cppcoreguidelines-pro-bounds-constant-array-index) -- k < nInt_ == intervals[0].err_.size() == intervals[0].quad_.size() by construction
                   }
                   if (maxRelErr < relTol_) { return intervals[0].quad_; }
               }
   
               std::vector<double> quadSum = intervals[0].quad_;
               std::vector<double> errSum = intervals[0].err_;
               std::size_t itCounter = 1;
               std::size_t intervalPtr = 0;
   
               while (true)
               {
                   const double xMid = 0.5 * (intervals[intervalPtr].a_ + intervals[intervalPtr].b_);
                   const double bOld = intervals[intervalPtr].b_; // intervals[intervalPtr] itself is about to be overwritten below, taking its own b_ with it -- the right half bisected off it still needs the original upper limit, not the new one
                   for (std::size_t k = 0; k < nInt_; ++k)
                   {
                       quadSum[k] -= intervals[intervalPtr].quad_[k]; // NOLINT(cppcoreguidelines-pro-bounds-constant-array-index) -- k < nInt_ == quadSum.size() == intervals[intervalPtr].quad_.size() by construction
                       errSum[k] -= intervals[intervalPtr].err_[k]; // NOLINT(cppcoreguidelines-pro-bounds-constant-array-index) -- see above, for errSum/err_
                   }
   
                   intervals[intervalPtr] = callQuadSingle(intervals[intervalPtr].a_, xMid);
                   intervals.push_back(callQuadSingle(xMid, bOld));
   
                   const auto& left = intervals[intervalPtr];
                   const auto& right = intervals.back();
                   for (std::size_t k = 0; k < nInt_; ++k)
                   {
                       quadSum[k] += left.quad_[k] + right.quad_[k]; // NOLINT(cppcoreguidelines-pro-bounds-constant-array-index) -- k < nInt_ == quadSum.size() == left.quad_.size() == right.quad_.size() by construction
                       errSum[k] += left.err_[k] + right.err_[k]; // NOLINT(cppcoreguidelines-pro-bounds-constant-array-index) -- see above, for errSum/err_
                   }
   
                   ++itCounter;
                   if (maxIter_ != 0 && itCounter > maxIter_) { return quadSum; }
   
                   double maxAbsErr = 0.0;
                   double maxRelErr = 0.0;
                   for (std::size_t k = 0; k < nInt_; ++k)
                   {
                       maxAbsErr = std::max(maxAbsErr, std::abs(errSum[k])); // NOLINT(cppcoreguidelines-pro-bounds-constant-array-index) -- k < nInt_ == errSum.size() by construction
                       maxRelErr = std::max(maxRelErr, std::abs(errSum[k] / quadSum[k])); // NOLINT(cppcoreguidelines-pro-bounds-constant-array-index) -- see above, for quadSum
                   }
                   if (maxAbsErr < absTol_ || maxRelErr < relTol_) { return quadSum; }
   
                   double bestRelErr = -1.0;
                   std::size_t best = 0;
                   for (std::size_t i = 0; i < intervals.size(); ++i)
                   {
                       double relErr = 0.0;
                       for (std::size_t k = 0; k < nInt_; ++k)
                       {
                           relErr = std::max(relErr,
                               std::abs(intervals[i].err_[k] / quadSum[k])); // NOLINT(cppcoreguidelines-pro-bounds-constant-array-index) -- k < nInt_ == intervals[i].err_.size() == quadSum.size() by construction
                       }
                       if (relErr > bestRelErr) { bestRelErr = relErr; best = i; }
                   }
                   intervalPtr = best;
               }
           }
   
       private:
   
           template <class... CallArgs>
           [[nodiscard]] auto invokeF(const double x, CallArgs&&... args) const
           {
               if constexpr (std::is_member_function_pointer_v<F>)
               {
                   return invokeFMember(x, std::forward<CallArgs>(args)...);
               }
               else
               {
                   return std::invoke(f_, x, std::forward<CallArgs>(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)...);
           }
   
           F f_;                  
           std::size_t nInt_;     
           std::size_t maxIter_;  
           double absTol_;        
           double relTol_;        
       };
   
   } // namespace utils
   
   #endif // GKINTEGRATOR_HPP