From bdf8d2b9c34946032065ae9ebc634bf6979ba029 Mon Sep 17 00:00:00 2001 From: Alexandre Hoffmann Date: Fri, 24 Jul 2026 18:21:53 +0200 Subject: [PATCH 1/7] Refactor adaptive quadrature interval handling Updated the scaling factor and improved interval handling in the adaptive quadrature implementation. Added checks for finite integrals and adjusted the logic for adding intervals to capture the entire support. --- .../AdaptiveQuadratureBase_impl.hpp | 84 ++++++++++--------- 1 file changed, 46 insertions(+), 38 deletions(-) diff --git a/include/LNIT/AdaptiveQuadratures/AdaptiveQuadratureBase_impl.hpp b/include/LNIT/AdaptiveQuadratures/AdaptiveQuadratureBase_impl.hpp index f7caedf..a71220a 100644 --- a/include/LNIT/AdaptiveQuadratures/AdaptiveQuadratureBase_impl.hpp +++ b/include/LNIT/AdaptiveQuadratures/AdaptiveQuadratureBase_impl.hpp @@ -115,45 +115,23 @@ auto AdaptiveQuadratureBase::integrateWithHints(Function&& f, const std using Iterator = typename std::vector::iterator; - constexpr Size scal = 3; + const Scalar scal = 2.58; assert(sigma > Scalar{}); - - GaussLaguerreQuadrature gLaguerreQuad; m_intervals.clear(); m_subIntergrals.clear(); m_subIntergralsErr.clear(); - m_intervals.reserve(2*mu.size()); + m_intervals.reserve(2*mu.size() + 2); - // First pass, we compute intervals near the peaks on which the function is not null. + // First pass, we compute intervals near the peaks. for (const double& mu_i : mu) { Interval curr(mu_i - scal*sigma, mu_i + scal*sigma); if (curr.first > curr.second) { swap(curr.first, curr.second); } - LongScalar leftIntegral = gLaguerreQuad.integrateLeftInfinite(f, curr.first); - while (isfinite(leftIntegral) and abs(leftIntegral) >= NumTraits::epsilon) - { - curr.first -= scal*sigma; - leftIntegral = gLaguerreQuad.integrateLeftInfinite(f, curr.first); - assert(curr.first <= curr.second); - } - - if (not isfinite(leftIntegral)) { return NumTraits::NaN; } - - LongScalar rightIntegral = gLaguerreQuad.integrateRightInfinite(f, curr.second); - while (isfinite(rightIntegral) and abs(rightIntegral) >= NumTraits::epsilon) - { - curr.second += scal*sigma; - rightIntegral = gLaguerreQuad.integrateRightInfinite(f, curr.second); - assert(curr.first <= curr.second); - } - - if (not isfinite(rightIntegral)) { return NumTraits::NaN; } - Iterator firstInterval = m_intervals.begin(); while (firstInterval != m_intervals.end() and firstInterval->second < curr.first) { ++firstInterval; } @@ -168,6 +146,43 @@ auto AdaptiveQuadratureBase::integrateWithHints(Function&& f, const std m_intervals.insert(it, curr); } + // Second pass: add interval in between the previously computed intervals + for (size_t i=0; i+1!=m_intervals.size(); ++i) + { + assert(m_intervals[i].second != m_intervals[i+1].first); + ret.emplace(std::next(m_intervals.begin(), i + 1), m_intervals[i].second, m_intervals[i+1].first); + } + + // now add interval at the front and rear to ensure we capture the whole support + GaussLaguerreQuadrature gLaguerreQuad; + + Scalar xmin = m_intervals.front().first - sigma; + LongScalar leftIntegral = gLaguerreQuad.integrateLeftInfinite(f, xmin); + + while (isfinite(leftIntegral) and abs(leftIntegral) >= NumTraits::epsilon) + { + xmin -= sigma; + leftIntegral = gLaguerreQuad.integrateLeftInfinite(f, xmin); + } + + if (not isfinite(leftIntegral)) { return NumTraits::NaN; } + + m_intervals.emplace(m_intervals.begin(), xmin, m_intervals.front().first); + + Scalar xmax = m_intervals.back().second + sigma; + LongScalar rightIntegral = gLaguerreQuad.integrateRightInfinite(f, xmax); + while (isfinite(rightIntegral) and abs(rightIntegral) >= NumTraits::epsilon) + { + xmax += sigma; + rightIntegral = gLaguerreQuad.integrateRightInfinite(f, xmax); + } + + if (not isfinite(rightIntegral)) { return NumTraits::NaN; } + + m_intervals.emplace_back(m_intervals.back().second, xmax); + + // Finally setup the adaptive quadrature. + LongScalar res; LongScalar estimatedErr; @@ -201,11 +216,10 @@ auto AdaptiveQuadratureBase::integrateLeftInfinite(Function&& f, const xmin *= 2; leftIntegral = gLaguerreQuad.integrateLeftInfinite(f, xmin); } - const LongScalar ret = isfinite(leftIntegral) + + return isfinite(leftIntegral) ? integrate(std::forward(f), xmin, xmax) : NumTraits::NaN; - - return ret; } template template @@ -223,11 +237,9 @@ auto AdaptiveQuadratureBase::integrateRightInfinite(Function&& f, const xmax *= 2; rightIntegral = gLaguerreQuad.integrateRightInfinite(f, xmax); } - const LongScalar ret = isfinite(rightIntegral) + return isfinite(rightIntegral) ? integrate(std::forward(f), xmin, xmax) : NumTraits::NaN; - - return ret; } template template @@ -246,11 +258,9 @@ auto AdaptiveQuadratureBase::integrate(Function&& f) -> LongScalar leftIntegral = gLaguerreQuad.integrateLeftInfinite(f, xmin); } - const LongScalar ret = isfinite(leftIntegral) + return isfinite(leftIntegral) ? integrateRightInfinite(std::forward(f), xmin) : NumTraits::NaN; - - return ret; } template template @@ -262,11 +272,9 @@ auto AdaptiveQuadratureBase::remapAndIntegrate(Function&& f) -> LongSca { const LongScalar fx = f(t / (1 - t*t)); const Scalar dxdt = (1 + t*t) / ((1 - t*t)*(1 - t*t)); - const LongScalar ret = isnan(fx*dxdt) + return isnan(fx*dxdt) ? LongScalar{} - : fx*dxdt; - - return ret; + : fx*dxdt; }, -1, 1); } From 691c6506d0db6039c51b936ead8bd4d0daae5056 Mon Sep 17 00:00:00 2001 From: Alexandre Hoffmann Date: Fri, 24 Jul 2026 18:26:57 +0200 Subject: [PATCH 2/7] Fix interval handling in AdaptiveQuadratureBase_impl --- .../AdaptiveQuadratures/AdaptiveQuadratureBase_impl.hpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/include/LNIT/AdaptiveQuadratures/AdaptiveQuadratureBase_impl.hpp b/include/LNIT/AdaptiveQuadratures/AdaptiveQuadratureBase_impl.hpp index a71220a..6a771f0 100644 --- a/include/LNIT/AdaptiveQuadratures/AdaptiveQuadratureBase_impl.hpp +++ b/include/LNIT/AdaptiveQuadratures/AdaptiveQuadratureBase_impl.hpp @@ -150,7 +150,7 @@ auto AdaptiveQuadratureBase::integrateWithHints(Function&& f, const std for (size_t i=0; i+1!=m_intervals.size(); ++i) { assert(m_intervals[i].second != m_intervals[i+1].first); - ret.emplace(std::next(m_intervals.begin(), i + 1), m_intervals[i].second, m_intervals[i+1].first); + m_intervals.emplace(std::next(m_intervals.begin(), i + 1), m_intervals[i].second, m_intervals[i+1].first); } // now add interval at the front and rear to ensure we capture the whole support @@ -189,9 +189,9 @@ auto AdaptiveQuadratureBase::integrateWithHints(Function&& f, const std m_subIntergrals.reserve(m_intervals.size()); m_subIntergralsErr.reserve(m_intervals.size()); - for (const auto& [xmin, xmax] : m_intervals) + for (const auto& [a, b] : m_intervals) { - std::tie(res, estimatedErr) = estimateIntegral(f, xmin, xmax); + std::tie(res, estimatedErr) = estimateIntegral(f, a, b); m_subIntergrals.push_back(res); m_subIntergralsErr.push_back(estimatedErr); From 158a069ac976c79e7224cd3e4c8a7c27ccc582c9 Mon Sep 17 00:00:00 2001 From: Alexandre Hoffmann Date: Fri, 24 Jul 2026 18:35:20 +0200 Subject: [PATCH 3/7] Potential fix for pull request finding Co-authored-by: Copilot Autofix powered by AI <175728472+Copilot@users.noreply.github.com> --- .../LNIT/AdaptiveQuadratures/AdaptiveQuadratureBase_impl.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/include/LNIT/AdaptiveQuadratures/AdaptiveQuadratureBase_impl.hpp b/include/LNIT/AdaptiveQuadratures/AdaptiveQuadratureBase_impl.hpp index 6a771f0..a4a3c0d 100644 --- a/include/LNIT/AdaptiveQuadratures/AdaptiveQuadratureBase_impl.hpp +++ b/include/LNIT/AdaptiveQuadratures/AdaptiveQuadratureBase_impl.hpp @@ -126,7 +126,7 @@ auto AdaptiveQuadratureBase::integrateWithHints(Function&& f, const std m_intervals.reserve(2*mu.size() + 2); // First pass, we compute intervals near the peaks. - for (const double& mu_i : mu) + for (const Scalar& mu_i : mu) { Interval curr(mu_i - scal*sigma, mu_i + scal*sigma); From e4c634bc6196311098625712af87b6c8836c4314 Mon Sep 17 00:00:00 2001 From: Alexandre Hoffmann Date: Fri, 24 Jul 2026 18:40:56 +0200 Subject: [PATCH 4/7] Fix loop condition in interval insertion logic --- .../LNIT/AdaptiveQuadratures/AdaptiveQuadratureBase_impl.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/include/LNIT/AdaptiveQuadratures/AdaptiveQuadratureBase_impl.hpp b/include/LNIT/AdaptiveQuadratures/AdaptiveQuadratureBase_impl.hpp index a4a3c0d..739c6a7 100644 --- a/include/LNIT/AdaptiveQuadratures/AdaptiveQuadratureBase_impl.hpp +++ b/include/LNIT/AdaptiveQuadratures/AdaptiveQuadratureBase_impl.hpp @@ -147,7 +147,7 @@ auto AdaptiveQuadratureBase::integrateWithHints(Function&& f, const std } // Second pass: add interval in between the previously computed intervals - for (size_t i=0; i+1!=m_intervals.size(); ++i) + for (size_t i=0; i+1 Date: Fri, 24 Jul 2026 18:45:33 +0200 Subject: [PATCH 5/7] Fix loop condition to use m_intervals size --- .../LNIT/AdaptiveQuadratures/AdaptiveQuadratureBase_impl.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/include/LNIT/AdaptiveQuadratures/AdaptiveQuadratureBase_impl.hpp b/include/LNIT/AdaptiveQuadratures/AdaptiveQuadratureBase_impl.hpp index 739c6a7..3ad8afa 100644 --- a/include/LNIT/AdaptiveQuadratures/AdaptiveQuadratureBase_impl.hpp +++ b/include/LNIT/AdaptiveQuadratures/AdaptiveQuadratureBase_impl.hpp @@ -147,7 +147,7 @@ auto AdaptiveQuadratureBase::integrateWithHints(Function&& f, const std } // Second pass: add interval in between the previously computed intervals - for (size_t i=0; i+1 Date: Fri, 24 Jul 2026 18:54:47 +0200 Subject: [PATCH 6/7] Change scal to constexpr for better optimization --- .../LNIT/AdaptiveQuadratures/AdaptiveQuadratureBase_impl.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/include/LNIT/AdaptiveQuadratures/AdaptiveQuadratureBase_impl.hpp b/include/LNIT/AdaptiveQuadratures/AdaptiveQuadratureBase_impl.hpp index 3ad8afa..adea15c 100644 --- a/include/LNIT/AdaptiveQuadratures/AdaptiveQuadratureBase_impl.hpp +++ b/include/LNIT/AdaptiveQuadratures/AdaptiveQuadratureBase_impl.hpp @@ -115,7 +115,7 @@ auto AdaptiveQuadratureBase::integrateWithHints(Function&& f, const std using Iterator = typename std::vector::iterator; - const Scalar scal = 2.58; + constexpr Scalar scal = Scalar(2.58); assert(sigma > Scalar{}); From d0bc7c61f0ac10afe2394736a21b94453f26b9eb Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Fri, 24 Jul 2026 16:56:20 +0000 Subject: [PATCH 7/7] Add early guard in integrateWithHints for empty mu --- .../LNIT/AdaptiveQuadratures/AdaptiveQuadratureBase_impl.hpp | 2 ++ 1 file changed, 2 insertions(+) diff --git a/include/LNIT/AdaptiveQuadratures/AdaptiveQuadratureBase_impl.hpp b/include/LNIT/AdaptiveQuadratures/AdaptiveQuadratureBase_impl.hpp index adea15c..8f2069a 100644 --- a/include/LNIT/AdaptiveQuadratures/AdaptiveQuadratureBase_impl.hpp +++ b/include/LNIT/AdaptiveQuadratures/AdaptiveQuadratureBase_impl.hpp @@ -123,6 +123,8 @@ auto AdaptiveQuadratureBase::integrateWithHints(Function&& f, const std m_subIntergrals.clear(); m_subIntergralsErr.clear(); + if (mu.empty()) { return integrate(std::forward(f)); } + m_intervals.reserve(2*mu.size() + 2); // First pass, we compute intervals near the peaks.