Program Listing for File GKIntegrator.hpp
↰ Return to documentation for file (src/utils/GKIntegrator.hpp)
#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