diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml new file mode 100644 index 0000000..cfab53f --- /dev/null +++ b/.github/workflows/ci.yml @@ -0,0 +1,86 @@ +name: C/C++ CI + +on: + push: + branches: [ "main" ] + pull_request: + branches: [ "main" ] + +jobs: + build: + strategy: + matrix: + os: [ubuntu-latest] + compiler: [g++-14, clang++-18] + build_type: [Debug, Release] + + runs-on: ${{ matrix.os }} + + steps: + - uses: actions/checkout@v4 + + - name: Install toolchain and dependencies + run: | + sudo apt-get update + sudo apt-get install -y g++-14 clang-18 libfmt-dev + + - name: Configure & build + run: | + cmake -B build \ + -DCMAKE_BUILD_TYPE=${{ matrix.build_type }} \ + -DCMAKE_CXX_COMPILER=${{ matrix.compiler }} \ + -DLNIT_BUILD_DEMO=ON \ + -DLNIT_BUILD_BENCHMARK=ON + cmake --build build --parallel + + benchmark: + # only track regressions from main; PRs/forks shouldn't push to gh-pages + if: github.event_name == 'push' && github.ref == 'refs/heads/main' + runs-on: ubuntu-latest + permissions: + contents: write + + steps: + - uses: actions/checkout@v4 + + - name: Install toolchain and dependencies + run: | + sudo apt-get update + sudo apt-get install -y g++-14 libfmt-dev + + - name: Configure & build benchmarks + run: | + cmake -B build \ + -DCMAKE_BUILD_TYPE=Release \ + -DCMAKE_CXX_COMPILER=g++-14 \ + -DLNIT_BUILD_BENCHMARK=ON + cmake --build build --parallel + + - name: Run benchmarks + run: | + ./build/benchmark/benchmark_addaptiveQuadratures --benchmark_format=json --benchmark_out=bench_addaptiveQuadratures_result.json + ./build/benchmark/benchmark_LevermoreLikePDF --benchmark_format=json --benchmark_out=bench_LevermoreLikePDF_result.json + + - name: Store adaptive quadratures benchmark & detect regressions + uses: benchmark-action/github-action-benchmark@v1 + with: + tool: 'googlecpp' + name: 'Adaptive Quadratures' + output-file-path: bench_addaptiveQuadratures_result.json + github-token: ${{ secrets.GITHUB_TOKEN }} + auto-push: true + alert-threshold: '150%' + comment-on-alert: true + fail-on-alert: true + + - name: Store Levermore-like PDF benchmark & detect regressions + uses: benchmark-action/github-action-benchmark@v1 + with: + tool: 'googlecpp' + name: 'Levermore-like PDF' + output-file-path: bench_LevermoreLikePDF_result.json + github-token: ${{ secrets.GITHUB_TOKEN }} + auto-push: true + alert-threshold: '150%' + comment-on-alert: true + fail-on-alert: true diff --git a/include/LNIT/AdaptiveQuadratures/AdaptiveQuadratureBase.hpp b/include/LNIT/AdaptiveQuadratures/AdaptiveQuadratureBase.hpp index 04629dd..5d95deb 100644 --- a/include/LNIT/AdaptiveQuadratures/AdaptiveQuadratureBase.hpp +++ b/include/LNIT/AdaptiveQuadratures/AdaptiveQuadratureBase.hpp @@ -65,7 +65,7 @@ class AdaptiveQuadratureBase * * @return Pair (integral, estimated error). */ - template constexpr std::pair estimateIntegral(const Function& f, const Scalar& xmin, const Scalar& xmax) { return derived().estimateIntegralImpl(f, xmin, xmax); } + template constexpr std::pair estimateIntegral(Function&& f, const Scalar& xmin, const Scalar& xmax) { return derived().estimateIntegralImpl(std::forward(f), xmin, xmax); } /** * @brief Perform adaptive quadrature on [xmin, xmax]. @@ -75,8 +75,10 @@ class AdaptiveQuadratureBase * @param xmax Upper bound of interval. * @return Approximation of the integral. */ - template LongScalar integrate(const Function& f, const Scalar& xmin, const Scalar& xmax); - + template LongScalar integrate(Function&& f, const Scalar& xmin, const Scalar& xmax); + + template LongScalar integrateWithHints(Function&& f, const std::span mu, const Scalar& sigma); + /** * @brief Perform adaptive quadrature on (-inf, xmax]. * @@ -92,7 +94,7 @@ class AdaptiveQuadratureBase * @param xmax Upper bound of interval. * @return Approximation of the integral. */ - template LongScalar integrateLeftInfinite(const Function& f, const Scalar& xmax); + template LongScalar integrateLeftInfinite(Function&& f, const Scalar& xmax); /** * @brief Perform adaptive quadrature on [xmin, inf). @@ -109,7 +111,7 @@ class AdaptiveQuadratureBase * @param xmax Upper bound of interval. * @return Approximation of the integral. */ - template LongScalar integrateRightInfinite(const Function& f, const Scalar& xmin); + template LongScalar integrateRightInfinite(Function&& f, const Scalar& xmin); /** * @brief Perform adaptive quadrature on (-inf, inf). @@ -124,7 +126,7 @@ class AdaptiveQuadratureBase * \f] * Finally addapt the quadrature over [xmin, xmax]. */ - template LongScalar integrate(const Function& f); + template LongScalar integrate(Function&& f); /** * @brief Perform adaptive quadrature on (-inf, inf) using a coordinate-remapping technique. @@ -139,9 +141,9 @@ class AdaptiveQuadratureBase * \int_{-1}^{1} f(x(t))\frac{1 + t^2}{(1 - t^2)^2} dt. * \f] */ - template LongScalar remapAndIntegrate(const Function& f); + template LongScalar remapAndIntegrate(Function&& f); - template std::invoke_result_t integrateWithoutAdaptation(const Function& f) const; + template std::invoke_result_t integrateWithoutAdaptation(Function&& f) const; constexpr Size getMaxIt() const { return m_maxIt; } ///< @brief Maximum iterations allowed. constexpr Size getNits() const { return m_it; } ///< @brief Number of iterations performed. @@ -168,6 +170,8 @@ class AdaptiveQuadratureBase constexpr std::span getSubIntervals() const { return m_intervals; } private: + template LongScalar adaptQuadrature(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..8f2069a 100644 --- a/include/LNIT/AdaptiveQuadratures/AdaptiveQuadratureBase_impl.hpp +++ b/include/LNIT/AdaptiveQuadratures/AdaptiveQuadratureBase_impl.hpp @@ -9,6 +9,7 @@ #include #include #include +#include #include @@ -28,15 +29,54 @@ 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(Function&& f) -> LongScalar +{ using std::abs; using std::isfinite; using const_Iterator = typename std::vector::const_iterator; + + LongScalar res; + LongScalar estimatedErr; 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 + 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(); +} + +template template +auto AdaptiveQuadratureBase::integrate(Function&& f, const Scalar& xmin, const Scalar& xmax) -> LongScalar +{ + using std::ceil; + m_intervals.clear(); m_subIntergrals.clear(); m_subIntergralsErr.clear(); @@ -44,15 +84,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,43 +101,114 @@ 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(std::forward(f)); +} + +template template +auto AdaptiveQuadratureBase::integrateWithHints(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 Scalar scal = Scalar(2.58); + + assert(sigma > Scalar{}); + + m_intervals.clear(); + 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. + for (const Scalar& 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)); } + Iterator firstInterval = m_intervals.begin(); + while (firstInterval != m_intervals.end() and firstInterval->second < curr.first) { ++firstInterval; } - if (not isfinite(I)) { return I; } - if (err < abs(I)*LongScalar(m_relativeTol) or err < LongScalar(m_absoluteTol)) { m_hasConverged = true; return I; } + 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); + } - // 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]; + const Iterator it = m_intervals.erase(firstInterval, boundInterval); + m_intervals.insert(it, curr); + } + + // Second pass: add interval in between the previously computed intervals + for (size_t i=0; i+1 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; + + m_subIntergrals.reserve(m_intervals.size()); + m_subIntergralsErr.reserve(m_intervals.size()); + + for (const auto& [a, b] : m_intervals) + { + std::tie(res, estimatedErr) = estimateIntegral(f, a, b); - 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(std::forward(f)); } template template -auto AdaptiveQuadratureBase::integrateLeftInfinite(const Function& f, const Scalar& xmax) -> LongScalar +auto AdaptiveQuadratureBase::integrateLeftInfinite(Function&& f, const Scalar& xmax) -> LongScalar { using std::isfinite; using std::abs; - GaussLaguerreQuadrature gLaguerreQuad; Scalar xmin = -1; @@ -110,11 +218,14 @@ auto AdaptiveQuadratureBase::integrateLeftInfinite(const Function& f, c xmin *= 2; leftIntegral = gLaguerreQuad.integrateLeftInfinite(f, xmin); } - return isfinite(leftIntegral) ? integrate(f, xmin, xmax) : NumTraits::NaN; + + return isfinite(leftIntegral) + ? integrate(std::forward(f), xmin, xmax) + : NumTraits::NaN; } template template -auto AdaptiveQuadratureBase::integrateRightInfinite(const Function& f, const Scalar& xmin) -> LongScalar +auto AdaptiveQuadratureBase::integrateRightInfinite(Function&& f, const Scalar& xmin) -> LongScalar { using std::isfinite; using std::abs; @@ -128,11 +239,13 @@ auto AdaptiveQuadratureBase::integrateRightInfinite(const Function& f, xmax *= 2; rightIntegral = gLaguerreQuad.integrateRightInfinite(f, xmax); } - return isfinite(rightIntegral) ? integrate(f, xmin, xmax) : NumTraits::NaN; + return isfinite(rightIntegral) + ? integrate(std::forward(f), xmin, xmax) + : NumTraits::NaN; } template template -auto AdaptiveQuadratureBase::integrate(const Function& f) -> LongScalar +auto AdaptiveQuadratureBase::integrate(Function&& f) -> LongScalar { using std::isfinite; using std::abs; @@ -147,31 +260,28 @@ auto AdaptiveQuadratureBase::integrate(const Function& f) -> LongScalar leftIntegral = gLaguerreQuad.integrateLeftInfinite(f, xmin); } - return isfinite(leftIntegral) ? integrateRightInfinite(f, xmin) : NumTraits::NaN; + return isfinite(leftIntegral) + ? integrateRightInfinite(std::forward(f), xmin) + : NumTraits::NaN; } template template -auto AdaptiveQuadratureBase::remapAndIntegrate(const Function& f) -> LongScalar +auto AdaptiveQuadratureBase::remapAndIntegrate(Function&& f) -> LongScalar { using std::isnan; - constexpr Scalar eps = {}; - - const auto fref = [&f](const Scalar t) -> LongScalar + return integrate([&f](const Scalar t) -> LongScalar { - const LongScalar fx = f(t / (1 - t*t)); - const Scalar dxdt = (1 + t*t) / ((1 - t*t)*(1 - t*t)); - + const LongScalar fx = f(t / (1 - t*t)); + const Scalar dxdt = (1 + t*t) / ((1 - t*t)*(1 - t*t)); return isnan(fx*dxdt) ? LongScalar{} - : fx*dxdt; - }; - - return integrate(fref, -1 + eps, 1 - eps); + : fx*dxdt; + }, -1, 1); } template template -auto AdaptiveQuadratureBase::integrateWithoutAdaptation(const Function& f) const -> std::invoke_result_t +auto AdaptiveQuadratureBase::integrateWithoutAdaptation(Function&& f) const -> std::invoke_result_t { const auto localIntegrals = m_intervals | std::views::transform([&self = derived(), &f](const Interval& interval) -> std::invoke_result_t { diff --git a/include/LNIT/AdaptiveQuadratures/ClenshawCurtisAdaptiveQuadrature.hpp b/include/LNIT/AdaptiveQuadratures/ClenshawCurtisAdaptiveQuadrature.hpp index 8597b59..e2832f3 100644 --- a/include/LNIT/AdaptiveQuadratures/ClenshawCurtisAdaptiveQuadrature.hpp +++ b/include/LNIT/AdaptiveQuadratures/ClenshawCurtisAdaptiveQuadrature.hpp @@ -33,9 +33,9 @@ class ClenshawCurtisAdaptiveQuadrature : public AdaptiveQuadratureBase< Clenshaw { using Base = AdaptiveQuadratureBase< ClenshawCurtisAdaptiveQuadrature >; public: - using Size = Base::Size; ///< @brief Type for iteration counters. - using Scalar = Base::Scalar; ///< @brief Floating point type for integration (e.g., double). - using LongScalar = Base::LongScalar; ///< @brief Higher precision type for accumulation (e.g., long double). + using Size = typename Base::Size; ///< @brief Type for iteration counters. + using Scalar = typename Base::Scalar; ///< @brief Floating point type for integration (e.g., double). + using LongScalar = typename Base::LongScalar; ///< @brief Higher precision type for accumulation (e.g., long double). /** * @brief Estimate integral and error on [xmin, xmax]. @@ -45,9 +45,9 @@ class ClenshawCurtisAdaptiveQuadrature : public AdaptiveQuadratureBase< Clenshaw * @param xmax Upper bound. * @return Pair (integral, estimated error). */ - template constexpr std::pair estimateIntegralImpl(const Function& f, const Scalar& xmin, const Scalar& xmax); + template constexpr std::pair estimateIntegralImpl(Function&& f, const Scalar& xmin, const Scalar& xmax); - template constexpr std::invoke_result_t integrateImpl(const Function& f, const Scalar& xmin, const Scalar& xmax) const; + template constexpr std::invoke_result_t integrateImpl(Function&& f, const Scalar& xmin, const Scalar& xmax) const; constexpr Scalar getMaxDeltaXImpl(const Scalar& xmin, const Scalar& xmax) const { return (xmax - xmin)*misc::maxDiff(std::span{s_xi}); } private: diff --git a/include/LNIT/AdaptiveQuadratures/ClenshawCurtisAdaptiveQuadrature_impl.hpp b/include/LNIT/AdaptiveQuadratures/ClenshawCurtisAdaptiveQuadrature_impl.hpp index 756c74b..4bb6e06 100644 --- a/include/LNIT/AdaptiveQuadratures/ClenshawCurtisAdaptiveQuadrature_impl.hpp +++ b/include/LNIT/AdaptiveQuadratures/ClenshawCurtisAdaptiveQuadrature_impl.hpp @@ -21,7 +21,7 @@ extern template class ClenshawCurtisAdaptiveQuadrature //// method implementations //// template template -constexpr auto ClenshawCurtisAdaptiveQuadrature::estimateIntegralImpl(const Function& f, const Scalar& xmin, const Scalar& xmax) -> std::pair +constexpr auto ClenshawCurtisAdaptiveQuadrature::estimateIntegralImpl(Function&& f, const Scalar& xmin, const Scalar& xmax) -> std::pair { using std::abs; @@ -52,7 +52,7 @@ constexpr auto ClenshawCurtisAdaptiveQuadrature::estimateIntegralImpl(cons } template template -constexpr auto ClenshawCurtisAdaptiveQuadrature::integrateImpl(const Function& f, const Scalar& xmin, const Scalar& xmax) const -> std::invoke_result_t +constexpr auto ClenshawCurtisAdaptiveQuadrature::integrateImpl(Function&& f, const Scalar& xmin, const Scalar& xmax) const -> std::invoke_result_t { const auto fx = s_xi | std::views::transform([&f, &xmin, &xmax](const Scalar& xi) -> std::invoke_result_t { diff --git a/include/LNIT/AdaptiveQuadratures/ClenshawCurtisHybridAdaptiveQuadrature.hpp b/include/LNIT/AdaptiveQuadratures/ClenshawCurtisHybridAdaptiveQuadrature.hpp index 2a88473..8ef1005 100644 --- a/include/LNIT/AdaptiveQuadratures/ClenshawCurtisHybridAdaptiveQuadrature.hpp +++ b/include/LNIT/AdaptiveQuadratures/ClenshawCurtisHybridAdaptiveQuadrature.hpp @@ -36,9 +36,9 @@ class ClenshawCurtisHybridAdaptiveQuadrature : public AdaptiveQuadratureBase< Cl { using Base = AdaptiveQuadratureBase< ClenshawCurtisHybridAdaptiveQuadrature >; public: - using Size = Base::Size; ///< @brief Type for iteration counters. - using Scalar = Base::Scalar; ///< @brief Floating point type for integration (e.g., double). - using LongScalar = Base::LongScalar; ///< @brief Higher precision type for accumulation (e.g., long double). + using Size = typename Base::Size; ///< @brief Type for iteration counters. + using Scalar = typename Base::Scalar; ///< @brief Floating point type for integration (e.g., double). + using LongScalar = typename Base::LongScalar; ///< @brief Higher precision type for accumulation (e.g., long double). /** * @brief Estimate integral and error on [xmin, xmax]. @@ -48,9 +48,9 @@ class ClenshawCurtisHybridAdaptiveQuadrature : public AdaptiveQuadratureBase< Cl * @param xmax Upper bound. * @return Pair (integral, estimated error). */ - template constexpr std::pair estimateIntegralImpl(const Function& f, const Scalar& xmin, const Scalar& xmax); + template constexpr std::pair estimateIntegralImpl(Function&& f, const Scalar& xmin, const Scalar& xmax); - template constexpr std::invoke_result_t integrateImpl(const Function& f, const Scalar& xmin, const Scalar& xmax) const; + template constexpr std::invoke_result_t integrateImpl(Function&& f, const Scalar& xmin, const Scalar& xmax) const; constexpr Scalar getMaxDeltaXImpl(const Scalar& xmin, const Scalar& xmax) const { return (xmax - xmin)*misc::maxDiff(std::span{s_xi}); } private: diff --git a/include/LNIT/AdaptiveQuadratures/ClenshawCurtisHybridAdaptiveQuadrature_impl.hpp b/include/LNIT/AdaptiveQuadratures/ClenshawCurtisHybridAdaptiveQuadrature_impl.hpp index 8247edf..34bec7a 100644 --- a/include/LNIT/AdaptiveQuadratures/ClenshawCurtisHybridAdaptiveQuadrature_impl.hpp +++ b/include/LNIT/AdaptiveQuadratures/ClenshawCurtisHybridAdaptiveQuadrature_impl.hpp @@ -20,7 +20,7 @@ extern template class ClenshawCurtisHybridAdaptiveQuadrature template -constexpr auto ClenshawCurtisHybridAdaptiveQuadrature::estimateIntegralImpl(const Function& f, const Scalar& xmin, const Scalar& xmax) -> std::pair +constexpr auto ClenshawCurtisHybridAdaptiveQuadrature::estimateIntegralImpl(Function&& f, const Scalar& xmin, const Scalar& xmax) -> std::pair { using std::abs; @@ -38,7 +38,7 @@ constexpr auto ClenshawCurtisHybridAdaptiveQuadrature::estimateIntegralImp } template template -constexpr auto ClenshawCurtisHybridAdaptiveQuadrature::integrateImpl(const Function& f, const Scalar& xmin, const Scalar& xmax) const -> std::invoke_result_t +constexpr auto ClenshawCurtisHybridAdaptiveQuadrature::integrateImpl(Function&& f, const Scalar& xmin, const Scalar& xmax) const -> std::invoke_result_t { const auto fx = s_xi | std::views::transform([&f, &xmin, &xmax](const Scalar& xi) -> std::invoke_result_t { diff --git a/include/LNIT/AdaptiveQuadratures/GLCCAdaptiveQuadrature.hpp b/include/LNIT/AdaptiveQuadratures/GLCCAdaptiveQuadrature.hpp index eb90c42..1c58758 100644 --- a/include/LNIT/AdaptiveQuadratures/GLCCAdaptiveQuadrature.hpp +++ b/include/LNIT/AdaptiveQuadratures/GLCCAdaptiveQuadrature.hpp @@ -37,9 +37,9 @@ class GLCCAdaptiveQuadrature : public AdaptiveQuadratureBase< GLCCAdaptiveQuadra { using Base = AdaptiveQuadratureBase< GLCCAdaptiveQuadrature >; public: - using Size = Base::Size; ///< @brief Type for iteration counters. - using Scalar = Base::Scalar; ///< @brief Floating point type for integration (e.g., double). - using LongScalar = Base::LongScalar; ///< @brief Higher precision type for accumulation (e.g., long double). + using Size = typename Base::Size; ///< @brief Type for iteration counters. + using Scalar = typename Base::Scalar; ///< @brief Floating point type for integration (e.g., double). + using LongScalar = typename Base::LongScalar; ///< @brief Higher precision type for accumulation (e.g., long double). /** * @brief Estimate integral and error on [xmin, xmax]. @@ -49,9 +49,9 @@ class GLCCAdaptiveQuadrature : public AdaptiveQuadratureBase< GLCCAdaptiveQuadra * @param xmax Upper bound. * @return Pair (integral, estimated error). */ - template constexpr std::pair estimateIntegralImpl(const Function& f, const Scalar& xmin, const Scalar& xmax) const; + template constexpr std::pair estimateIntegralImpl(Function&& f, const Scalar& xmin, const Scalar& xmax) const; - template constexpr std::invoke_result_t integrateImpl(const Function& f, const Scalar& xmin, const Scalar& xmax) const; + template constexpr std::invoke_result_t integrateImpl(Function&& f, const Scalar& xmin, const Scalar& xmax) const; constexpr Scalar getMaxDeltaXImpl(const Scalar xmin, const Scalar xmax) const { return (xmax - xmin)*misc::maxDiff(std::span{s_xi_cc}); } private: diff --git a/include/LNIT/AdaptiveQuadratures/GLCCAdaptiveQuadrature_impl.hpp b/include/LNIT/AdaptiveQuadratures/GLCCAdaptiveQuadrature_impl.hpp index 906dadc..7306575 100644 --- a/include/LNIT/AdaptiveQuadratures/GLCCAdaptiveQuadrature_impl.hpp +++ b/include/LNIT/AdaptiveQuadratures/GLCCAdaptiveQuadrature_impl.hpp @@ -20,7 +20,7 @@ extern template class GLCCAdaptiveQuadrature; //// method implementations //// template template -constexpr auto GLCCAdaptiveQuadrature::estimateIntegralImpl(const Function& f, const Scalar& xmin, const Scalar& xmax) const -> std::pair +constexpr auto GLCCAdaptiveQuadrature::estimateIntegralImpl(Function&& f, const Scalar& xmin, const Scalar& xmax) const -> std::pair { using std::abs; @@ -43,7 +43,7 @@ constexpr auto GLCCAdaptiveQuadrature::estimateIntegralImpl(const Function } template template -constexpr auto GLCCAdaptiveQuadrature::integrateImpl(const Function& f, const Scalar& xmin, const Scalar& xmax) const -> std::invoke_result_t +constexpr auto GLCCAdaptiveQuadrature::integrateImpl(Function&& f, const Scalar& xmin, const Scalar& xmax) const -> std::invoke_result_t { const auto fx = s_xi_gl | std::views::transform([&f, &xmin, &xmax](const Scalar& xi) -> std::invoke_result_t { diff --git a/include/LNIT/AdaptiveQuadratures/GaussLegendreAdaptiveQuadrature.hpp b/include/LNIT/AdaptiveQuadratures/GaussLegendreAdaptiveQuadrature.hpp index 4d861fc..4b15474 100644 --- a/include/LNIT/AdaptiveQuadratures/GaussLegendreAdaptiveQuadrature.hpp +++ b/include/LNIT/AdaptiveQuadratures/GaussLegendreAdaptiveQuadrature.hpp @@ -24,9 +24,9 @@ class GaussLegendreAdaptiveQuadrature : public AdaptiveQuadratureBase< GaussLege { using Base = AdaptiveQuadratureBase< GaussLegendreAdaptiveQuadrature >; public: - using Size = Base::Size; ///< @brief Type for iteration counters. - using Scalar = Base::Scalar; ///< @brief Floating point type for integration (e.g., double). - using LongScalar = Base::LongScalar; ///< @brief Higher precision type for accumulation (e.g., long double). + using Size = typename Base::Size; ///< @brief Type for iteration counters. + using Scalar = typename Base::Scalar; ///< @brief Floating point type for integration (e.g., double). + using LongScalar = typename Base::LongScalar; ///< @brief Higher precision type for accumulation (e.g., long double). /** * @brief Estimate integral and error on [xmin, xmax]. @@ -36,9 +36,9 @@ class GaussLegendreAdaptiveQuadrature : public AdaptiveQuadratureBase< GaussLege * @param xmax Upper bound. * @return Pair (integral, estimated error). */ - template constexpr std::pair estimateIntegralImpl(const Function& f, const Scalar& xmin, const Scalar& xmax); + template constexpr std::pair estimateIntegralImpl(Function&& f, const Scalar& xmin, const Scalar& xmax); - template constexpr std::invoke_result_t integrateImpl(const Function& f, const Scalar& xmin, const Scalar& xmax) const; + template constexpr std::invoke_result_t integrateImpl(Function&& f, const Scalar& xmin, const Scalar& xmax) const; inline constexpr Scalar getMaxDeltaXImpl(const Scalar xmin, const Scalar xmax) const { return (xmax - xmin)*misc::maxDiff(std::span{s_xi}); } private: diff --git a/include/LNIT/AdaptiveQuadratures/GaussLegendreAdaptiveQuadrature_impl.hpp b/include/LNIT/AdaptiveQuadratures/GaussLegendreAdaptiveQuadrature_impl.hpp index bf3fdff..a348c84 100644 --- a/include/LNIT/AdaptiveQuadratures/GaussLegendreAdaptiveQuadrature_impl.hpp +++ b/include/LNIT/AdaptiveQuadratures/GaussLegendreAdaptiveQuadrature_impl.hpp @@ -20,7 +20,7 @@ extern template class GaussLegendreAdaptiveQuadrature; //// method implementations //// template template -constexpr auto GaussLegendreAdaptiveQuadrature::estimateIntegralImpl(const Function& f, const Scalar& xmin, const Scalar& xmax) -> std::pair +constexpr auto GaussLegendreAdaptiveQuadrature::estimateIntegralImpl(Function&& f, const Scalar& xmin, const Scalar& xmax) -> std::pair { using std::abs; @@ -43,7 +43,7 @@ constexpr auto GaussLegendreAdaptiveQuadrature::estimateIntegralImpl(const } template template -constexpr auto GaussLegendreAdaptiveQuadrature::integrateImpl(const Function& f, const Scalar& xmin, const Scalar& xmax) const -> std::invoke_result_t +constexpr auto GaussLegendreAdaptiveQuadrature::integrateImpl(Function&& f, const Scalar& xmin, const Scalar& xmax) const -> std::invoke_result_t { const auto fx = s_xi | std::views::transform([&f, &xmin, &xmax](const Scalar& xi) -> std::invoke_result_t { diff --git a/include/LNIT/GaussHermiteQuadrature.hpp b/include/LNIT/GaussHermiteQuadrature.hpp index 9f9fafd..c0c0f6b 100644 --- a/include/LNIT/GaussHermiteQuadrature.hpp +++ b/include/LNIT/GaussHermiteQuadrature.hpp @@ -24,7 +24,7 @@ class GaussHermiteQuadrature constexpr GaussHermiteQuadrature() {} ///< @brief Default constructor. - template constexpr LongScalar integrate(const Function& f) const; + template constexpr LongScalar integrate(Function&& f) const; private: static constexpr std::array s_wi = { Scalar(0.72500244352094499799), Scalar(0.55908539234713480356), Scalar(0.48645731979545141956), Scalar(0.44287166385270444111), Scalar(0.41289257970879246181), diff --git a/include/LNIT/GaussHermiteQuadrature_impl.hpp b/include/LNIT/GaussHermiteQuadrature_impl.hpp index bd88aab..904c196 100644 --- a/include/LNIT/GaussHermiteQuadrature_impl.hpp +++ b/include/LNIT/GaussHermiteQuadrature_impl.hpp @@ -22,14 +22,14 @@ extern template class GaussHermiteQuadrature; //// method implementations //// template template -constexpr LongScalar GaussHermiteQuadrature::integrate(const Function& f) const +constexpr LongScalar GaussHermiteQuadrature::integrate(Function&& f) const { const auto fx = s_xi | std::views::transform([&f](const Scalar& x) -> LongScalar { return f(x); }); - return std::inner_product(s_wi.begin(), s_wi.end(), fx.end(), LongScalar{}); + return std::inner_product(s_wi.begin(), s_wi.end(), fx.begin(), LongScalar{}); } } // namespace LNIT diff --git a/include/LNIT/GaussLaguerreQuadrature.hpp b/include/LNIT/GaussLaguerreQuadrature.hpp index 1d0522e..ed3df28 100644 --- a/include/LNIT/GaussLaguerreQuadrature.hpp +++ b/include/LNIT/GaussLaguerreQuadrature.hpp @@ -38,7 +38,7 @@ class GaussLaguerreQuadrature * @param a Upper bound of integration (default = 0). * @return Approximated integral value. */ - template constexpr LongScalar integrateLeftInfinite (const Function& f, const Scalar& a = Scalar{}) const; + template constexpr LongScalar integrateLeftInfinite (Function&& f, const Scalar& a = Scalar{}) const; /** * @brief Approximate integral over the right semi-infinite interval. * @@ -53,7 +53,7 @@ class GaussLaguerreQuadrature * @param a Lower bound of integration (default = 0). * @return Approximated integral value. */ - template constexpr LongScalar integrateRightInfinite (const Function& f, const Scalar& a = Scalar{}) const; + template constexpr LongScalar integrateRightInfinite (Function&& f, const Scalar& a = Scalar{}) const; private: static constexpr std::array s_wi = { Scalar(0.11077730587320757274), Scalar(0.25810528128189475158), Scalar(0.40622176868437369247), Scalar(0.55526230959922306292), Scalar(0.70555738765958285661), diff --git a/include/LNIT/GaussLaguerreQuadrature_impl.hpp b/include/LNIT/GaussLaguerreQuadrature_impl.hpp index 9fa7486..526ce53 100644 --- a/include/LNIT/GaussLaguerreQuadrature_impl.hpp +++ b/include/LNIT/GaussLaguerreQuadrature_impl.hpp @@ -18,25 +18,25 @@ extern template class GaussLaguerreQuadrature; //// method implementations //// template template -constexpr LongScalar GaussLaguerreQuadrature::integrateLeftInfinite(const Function& f, const Scalar& a) const +constexpr LongScalar GaussLaguerreQuadrature::integrateLeftInfinite(Function&& f, const Scalar& a) const { const auto fx = s_xi | std::views::transform([&f, &a](const Scalar& x) -> LongScalar { return f(a - x); }); - return std::inner_product(s_wi.begin(), s_wi.end(), fx.end(), LongScalar{}); + return std::inner_product(s_wi.begin(), s_wi.end(), fx.begin(), LongScalar{}); } template template -constexpr LongScalar GaussLaguerreQuadrature::integrateRightInfinite(const Function& f, const Scalar& a) const +constexpr LongScalar GaussLaguerreQuadrature::integrateRightInfinite(Function&& f, const Scalar& a) const { const auto fx = s_xi | std::views::transform([&f, &a](const Scalar& x) -> LongScalar { return f(x + a); }); - return std::inner_product(s_wi.begin(), s_wi.end(), fx.end(), LongScalar{}); + return std::inner_product(s_wi.begin(), s_wi.end(), fx.begin(), LongScalar{}); } } // namespace LNIT