diff --git a/include/LNIT/AdaptiveQuadratures/AdaptiveQuadratureBase.hpp b/include/LNIT/AdaptiveQuadratures/AdaptiveQuadratureBase.hpp index 04629dd..799e2ee 100644 --- a/include/LNIT/AdaptiveQuadratures/AdaptiveQuadratureBase.hpp +++ b/include/LNIT/AdaptiveQuadratures/AdaptiveQuadratureBase.hpp @@ -76,7 +76,9 @@ class AdaptiveQuadratureBase * @return Approximation of the integral. */ template LongScalar integrate(const Function& f, const Scalar& xmin, const Scalar& xmax); - + + template LongScalar integrateWithHints(const Function& f, const std::span mu, const Scalar& sigma); + /** * @brief Perform adaptive quadrature on (-inf, xmax]. * @@ -168,6 +170,8 @@ class AdaptiveQuadratureBase constexpr std::span getSubIntervals() const { return m_intervals; } private: + template LongScalar adaptQuadrature(const Function& func); + std::vector m_intervals; std::vector m_subIntergrals; std::vector m_subIntergralsErr; diff --git a/include/LNIT/AdaptiveQuadratures/AdaptiveQuadratureBase_impl.hpp b/include/LNIT/AdaptiveQuadratures/AdaptiveQuadratureBase_impl.hpp index 907469c..8297c1f 100644 --- a/include/LNIT/AdaptiveQuadratures/AdaptiveQuadratureBase_impl.hpp +++ b/include/LNIT/AdaptiveQuadratures/AdaptiveQuadratureBase_impl.hpp @@ -28,15 +28,51 @@ AdaptiveQuadratureBase::AdaptiveQuadratureBase(const Size& maxIt, const } template template -auto AdaptiveQuadratureBase::integrate(const Function& f, const Scalar& xmin, const Scalar& xmax) -> LongScalar -{ - using std::ceil; +auto AdaptiveQuadratureBase::adaptQuadrature(const Function& f) -> LongScalar +{ using std::abs; using std::isfinite; using const_Iterator = typename std::vector::const_iterator; - + m_hasConverged = false; + + if (m_out) { fmt::print(m_out, "#Iteration integral estimated_error relative_tol absolute_tol\n"); } + + for (m_it=0; m_it!=m_maxIt; ++m_it) + { + const LongScalar I = getEstimatedIntegral(); + const LongScalar err = getEstimatedError(); + + if (m_out) { fmt::print(m_out, "{} {:10.4e} {:10.4e} {:10.4e} {:10.4e}\n", m_it, Scalar(I), Scalar(err), Scalar(abs(I))*m_relativeTol, Scalar(m_absoluteTol)); } + + if (not isfinite(I)) { return I; } + if (err < abs(I)*LongScalar(m_relativeTol) or err < LongScalar(m_absoluteTol)) { m_hasConverged = true; return I; } + + // we find the interval over which the integral is the least accurate + const const_Iterator maxErrIt = std::ranges::max_element(m_subIntergralsErr); + const Size maxErrIdx = Size(std::ranges::distance(m_subIntergralsErr.begin(), maxErrIt)); + // we split it in two + const auto& [a, b] = m_intervals[maxErrIdx]; + + Scalar midPoint = std::midpoint(a, b); // non-const because I want to move it when I do not need it. + // first interval + m_intervals[maxErrIdx] = Interval(a, midPoint); + std::tie(m_subIntergrals[maxErrIdx], m_subIntergralsErr[maxErrIdx]) = estimateIntegral(f, a, midPoint); + // second interval + const auto [res, estimatedErr] = estimateIntegral(f, midPoint, b); + m_intervals.emplace_back(std::move(midPoint), b); + m_subIntergrals.push_back(res); + m_subIntergralsErr.push_back(estimatedErr); + } + return getEstimatedIntegral(); +} + +template template +auto AdaptiveQuadratureBase::integrate(const Function& f, const Scalar& xmin, const Scalar& xmax) -> LongScalar +{ + using std::ceil; + m_intervals.clear(); m_subIntergrals.clear(); m_subIntergralsErr.clear(); @@ -44,15 +80,12 @@ auto AdaptiveQuadratureBase::integrate(const Function& f, const Scalar& LongScalar res; LongScalar estimatedErr; - if (m_out) { fmt::print(m_out, "#NumericalIntegrator addapting quadrature over [{}, {}]\n", xmin, xmax); } - if (m_out) { fmt::print(m_out, "#Iteration integral estimated_error relative_tol absolute_tol\n"); } - const Size N = Size(ceil(getMaxDeltaX(xmin, xmax))); m_intervals.reserve(N); m_subIntergrals.reserve(N); m_subIntergralsErr.reserve(N); - + for (Size i=0; i!=N; ++i) { const Scalar x_i = xmin + Scalar(i)*(xmax - xmin) / Scalar(N); @@ -64,34 +97,89 @@ auto AdaptiveQuadratureBase::integrate(const Function& f, const Scalar& m_subIntergrals.push_back(res); m_subIntergralsErr.push_back(estimatedErr); } + + if (m_out) { fmt::print(m_out, "#NumericalIntegrator adapting quadrature over [{}, {}]\n", xmin, xmax); } + return adaptQuadrature(f); +} + +template template +auto AdaptiveQuadratureBase::integrateWithHints(const Function& f, const std::span mu, const Scalar& sigma) -> LongScalar +{ + using std::abs; + using std::isfinite; + using std::swap; + + using Iterator = typename std::vector::iterator; - for (m_it=0; m_it!=m_maxIt; ++m_it) + constexpr Size scal = 3; + + assert(sigma > Scalar{}); + + GaussLaguerreQuadrature gLaguerreQuad; + + m_intervals.clear(); + m_subIntergrals.clear(); + m_subIntergralsErr.clear(); + + m_intervals.reserve(2*mu.size()); + + // First pass, we compute intervals near the peaks on which the function is not null. + for (const double& mu_i : mu) { - const LongScalar I = getEstimatedIntegral(); - const LongScalar err = getEstimatedError(); + Interval curr(mu_i - scal*sigma, mu_i + scal*sigma); + + if (curr.first > curr.second) { swap(curr.first, curr.second); } - if (m_out) { fmt::print(m_out, "{} {:10.4e} {:10.4e} {:10.4e} {:10.4e}\n", m_it, Scalar(I), Scalar(err), Scalar(abs(I))*m_relativeTol, Scalar(m_absoluteTol)); } + 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(I)) { return I; } - if (err < abs(I)*LongScalar(m_relativeTol) or err < LongScalar(m_absoluteTol)) { m_hasConverged = true; return I; } + if (not isfinite(leftIntegral)) { return NumTraits::NaN; } - // we find the interval over which the integral is the least accurate - const const_Iterator maxErrIt = std::ranges::max_element(m_subIntergralsErr); - const Size maxErrIdx = Size(std::ranges::distance(m_subIntergralsErr.begin(), maxErrIt)); - // we split it in two - const auto& [a, b] = m_intervals[maxErrIdx]; + 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; } + + Iterator boundInterval = firstInterval; + for ( ;boundInterval != m_intervals.end() and boundInterval->first <= curr.second; ++boundInterval) + { + curr.first = std::min(curr.first, boundInterval->first); + curr.second = std::max(curr.second, boundInterval->second); + } + + const Iterator it = m_intervals.erase(firstInterval, boundInterval); + m_intervals.insert(it, curr); + } + + LongScalar res; + LongScalar estimatedErr; + + m_subIntergrals.reserve(m_intervals.size()); + m_subIntergralsErr.reserve(m_intervals.size()); + + for (const auto& [xmin, xmax] : m_intervals) + { + std::tie(res, estimatedErr) = estimateIntegral(f, xmin, xmax); - Scalar midPoint = std::midpoint(a, b); // non-const because I want to move it when I do not need it. - // first interval - m_intervals[maxErrIdx] = Interval(a, midPoint); - std::tie(m_subIntergrals[maxErrIdx], m_subIntergralsErr[maxErrIdx]) = estimateIntegral(f, a, midPoint); - // second interval - std::tie(res, estimatedErr) = estimateIntegral(f, midPoint, b); - m_intervals.emplace_back(std::move(midPoint), b); m_subIntergrals.push_back(res); m_subIntergralsErr.push_back(estimatedErr); } - return getEstimatedIntegral(); + + if (m_out) { fmt::print(m_out, "#NumericalIntegrator adapting quadrature over {}\n", m_intervals); } + return adaptQuadrature(f); } template template @@ -100,7 +188,6 @@ auto AdaptiveQuadratureBase::integrateLeftInfinite(const Function& f, c using std::isfinite; using std::abs; - GaussLaguerreQuadrature gLaguerreQuad; Scalar xmin = -1;