From 22c4aed8ea3c3098dca74d3244752511616e2f0d Mon Sep 17 00:00:00 2001 From: Nick Thompson Date: Wed, 28 Oct 2026 22:59:05 -0700 Subject: [PATCH] Add nested Leja quadrature and a dimension-adaptive sparse grid Legendre, Hermite, Laguerre and Jacobi Leja rules reuse earlier evaluations; the sparse grid refines where contributions are largest. Co-Authored-By: Claude Opus 5.5 --- doc/math.qbk | 1 + doc/quadrature/leja.qbk | 184 +++++ example/Jamfile.v2 | 1 + example/leja_parallel_assembly.cpp | 111 +++ include/boost/math/quadrature/leja.hpp | 727 ++++++++++++++++++ include/boost/math/quadrature/sparse_grid.hpp | 353 +++++++++ nanobind/src/quadrature.cpp | 95 +++ nanobind/tests/test_smoke.py | 52 +- .../leja_quadrature_performance.cpp | 136 ++++ test/Jamfile.v2 | 1 + test/test_leja_quadrature.cpp | 708 +++++++++++++++++ 11 files changed, 2367 insertions(+), 2 deletions(-) create mode 100644 doc/quadrature/leja.qbk create mode 100644 example/leja_parallel_assembly.cpp create mode 100644 include/boost/math/quadrature/leja.hpp create mode 100644 include/boost/math/quadrature/sparse_grid.hpp create mode 100644 reporting/performance/leja_quadrature_performance.cpp create mode 100644 test/test_leja_quadrature.cpp diff --git a/doc/math.qbk b/doc/math.qbk index 5539bfb53a..5ed278b334 100644 --- a/doc/math.qbk +++ b/doc/math.qbk @@ -774,6 +774,7 @@ and as a CD ISBN 0-9504833-2-X 978-0-9504833-2-0, Classification 519.2-dc22. [include quadrature/gauss.qbk] [include quadrature/gauss_kronrod.qbk] [include quadrature/gauss_hermite.qbk] +[include quadrature/leja.qbk] [include quadrature/double_exponential.qbk] [include quadrature/ooura_fourier_integrals.qbk] [include quadrature/naive_monte_carlo.qbk] diff --git a/doc/quadrature/leja.qbk b/doc/quadrature/leja.qbk new file mode 100644 index 0000000000..0e7b79de0d --- /dev/null +++ b/doc/quadrature/leja.qbk @@ -0,0 +1,184 @@ +[/ +Copyright (c) 2026 Nick Thompson +Use, modification and distribution are subject to the +Boost Software License, Version 1.0. (See accompanying file +LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) +] + +[section:leja Leja quadrature] + +[heading Synopsis] + +`#include ` + + namespace boost{ namespace math{ namespace quadrature{ + + template > + class leja_legendre // the integral of f over [-1, 1] + { + public: + leja_legendre(); + + const std::vector& abscissa() const; + const std::vector& weights(unsigned n) const; + static constexpr unsigned max_points(); + + template + auto integrate(F f, Real tolerance = tools::root_epsilon(), + Real* error = nullptr, Real* L1 = nullptr) const -> decltype(f(Real(0))); + + // The integral over [a, b]; leja_legendre and leja_jacobi only. + template + auto integrate(F f, Real a, Real b, Real tolerance = tools::root_epsilon(), + Real* error = nullptr, Real* L1 = nullptr) const -> decltype(f(Real(0))); + }; + + template > + class leja_hermite; // the integral of f(x) exp(-x[super 2]) over the real line, with the same members except integrate(f, a, b) + + // The integral of f(x) x[super alpha] exp(-x) over \[0, [infin]), alpha > -1; as leja_hermite, with the constructor + // explicit leja_laguerre(Real alpha = 0); + template > + class leja_laguerre; + + // The integral of f(x) (1 - x)[super alpha] (1 + x)[super beta] over \[-1, 1\], alpha, beta > -1, with the constructor + // leja_jacobi(Real alpha, Real beta); + template > + class leja_jacobi; + + }}} // namespaces + +[heading Description] + +A weighted Leja sequence is built greedily: each new point maximizes [radic]/w/ times the product of its distances to the points +already chosen (Leja; Narayan and Jakeman). `leja_legendre` has the Legendre weight /w/(/x/) = 1 on \[-1, 1\], which gives +the classical Leja points, here in the order 0, -1, 1; `leja_hermite` has /w/(/x/) = exp(-/x/[super 2]), `leja_laguerre` +/w/(/x/) = /x/[super [alpha]] exp(-/x/) on \[0, [infin]), and `leja_jacobi` /w/(/x/) = (1 - /x/)[super [alpha]] (1 + /x/)[super [beta]] on \[-1, 1\]. +All start from the mean of /w/, so that the one-point rule is the one-point Gauss rule, exact to degree 1, the natural first level +of a sparse grid. Where an exponent is negative, [radic]/w/ has no maximum, and the points are chosen as if it were 0; the weights +are those of /w/. `leja_jacobi` with [alpha] = [beta] = 0 has the points of `leja_legendre`. +The constructors raise a __domain_error unless the exponents are > -1. +The /n/-point rule uses the first /n/ points and has the weights that make it exact for polynomials of degree /n/ - 1. + +The rules are nested: going from /n/ to /n/ + 1 points costs one more evaluation of /f/, where a Gauss rule must start again. +This makes them the usual building block of sparse grids, and lets the adaptive `integrate` add one point at a time. +The price is degree: a Leja rule with /n/ points is exact to degree /n/ - 1, a Gauss rule to 2/n/ - 1. + +The points and weights are computed when they are first needed, up to `max_points()`: 1023 for the built-in floating point +types and 127 for others. `weights(n)` returns the /n/ weights of the /n/-point rule, computing it and any smaller ones first, +and raises a __domain_error unless 1 [le] /n/ [le] `max_points()`. `abscissa()` returns the points computed so far, in +the order they were chosen. All members are `const` and safe to call from several threads, which share what any of them +computes, as `tanh_sinh` does; points that another thread is computing read as NaN until it is done. +`example/leja_parallel_assembly.cpp` assembles the Hamiltonian of the anharmonic oscillator /x/[super 4] in 60 Hermite functions +on all cores with one `leja_legendre`, in agreement with the exact entries to 2[times]10[super -15] times the largest: 0.21 s on one core +of an Apple M-series and 0.026 s on ten, once the rule has its points. +`leja_hermite` suits expectations over a normal distribution, such as the mean of the logistic function of /X/ ~ N(1, 1): + + leja_hermite rule; + double error; + double mean = rule.integrate([](double x) { return 1 / (1 + std::exp(-1 - std::sqrt(2.0) * x)); }, 1e-10, &error) + / std::sqrt(boost::math::constants::pi()); + +takes 105 points for 0.69673467014367874, within 4.6[times]10[super -15] of the true value; the estimate of its error, `error` divided by [radic][pi], is 5.6[times]10[super -12]. + +When several points maximize the product to within rounding, the leftmost is taken; in the tests float, double, long double and +`cpp_bin_float_50` choose the same points. Some weights are negative, but not by much: up to 100 points the sum of |/w[sub i]/| +is at most 1.32 times the total mass for the exponents tried, but 2.44 times for [alpha] = -0.9, [beta] = 1/2. + +`integrate(f, tolerance, error, L1)` evaluates /f/ once at each point in turn, until /n/ [ge] 4 and +|/Q[sub m]/ - /Q[sub m - max(1, m/4)]/| [le] `tolerance` * /L1/ for every /m/ from /n/ - /n//4 to /n/, where /Q[sub n]/ is the /n/-point rule +and /L1/ = [sum] |/w[sub i]/| |/f/(/x[sub i]/)|. If `error` or `L1` are not null they are set to the largest of those differences, +a heuristic estimate of the error, and to /L1/. +If `max_points()` is reached first, its result is returned with its error estimate, as `tanh_sinh` does; no error is raised. +/f/ may return a real or complex value, or a vector of them: any type with `operator[]` and `size()`, such as `std::array` or +`std::vector`, whose values all have the same size. A vector is integrated componentwise, and |/y/| above is its max norm. +`leja_legendre` evaluates /f/ at the endpoints [plusminus]1, and `leja_laguerre` and `leja_jacobi` at an end of the support where the exponent is [le] 0. + +`integrate(f, a, b, tolerance, error, L1)` integrates /f/ over \[/a/, /b/\] with `leja_legendre`, and /f/(/x/) (/b/ - /x/)[super [alpha]] (/x/ - /a/)[super [beta]] +with `leja_jacobi`, by the affine map from \[-1, 1\], which also scales the error and /L1/. If /a/ > /b/ the result is minus the +integral over \[/b/, /a/\], with [alpha] still at /b/. Infinite or NaN bounds raise a __domain_error. + +[heading Accuracy and cost] + +The rules up to /N/ points take O(/N/[super 2]) operations: about 0.3 ms for 64 points, 1 ms for 128, +16 ms for 512 and 60 ms for 1023 in double precision, and 0.4 s for 127 in `cpp_bin_float_50`. +Each rule follows from the last in Newton form, with the integral of the nodal polynomial from a recurrence carried to twice the +working precision: the rules integrate /x/[super /k/] to within 20 /n/ [epsilon] [sum] |/w[sub i]/ /x[sub i]/[super /k/]|, and in double +precision the weights are within an ulp of the largest weight of the same rules in 50 digits, up to 1023 points. +A bounded sum of |/w[sub i]/| does not help integrands that grow like exp(/x/[super 2]) against `leja_hermite`, such as products of Hermite polynomials +of high degree: /L1/ is 10[super 19] times the integral of /h/[sub 40][super 2]. Integrate such products with their Gaussian over +a finite interval instead, as the example does. + +[heading References] + +* F. Leja, ['Sur certaines suites liées aux ensembles plans et leur application à la représentation conforme], Annales Polonici Mathematici 4 (1957), 8-13. +* L. Reichel, ['Newton interpolation at Leja points], BIT 30 (1990), 332-346, [@https://doi.org/10.1007/BF02017352 doi:10.1007/BF02017352]. +* A. Narayan and J. D. Jakeman, ['Adaptive Leja sparse grid constructions for stochastic collocation and high-dimensional approximation], SIAM Journal on Scientific Computing 36 (2014), no. 6, A2952-A2983, [@https://doi.org/10.1137/140966368 doi:10.1137/140966368]. + +[endsect] [/section:leja Leja quadrature] + +[section:sparse_grid Sparse grids] + +[heading Synopsis] + +`#include ` + + namespace boost{ namespace math{ namespace quadrature{ + + // Rule: leja_legendre, leja_hermite, leja_laguerre or leja_jacobi, Real = typename Rule::value_type; f(const std::vector& x), or f(x) with x a std::array, Eigen::VectorXd or other container with operator[] and size(). + template + auto sparse_grid(F f, const Rule& rule, unsigned dimension, std::size_t max_calls = 100000, + Real tolerance = tools::root_epsilon(), Real* error = nullptr, Real* L1 = nullptr) + -> decltype(f(x)); + + // Over the box [a[0], b[0]] x ... x [a[d-1], b[d-1]], d = a.size(), for leja_legendre and leja_jacobi; Box has operator[] and size(). + template + auto sparse_grid(F f, const Rule& rule, const Box& a, const Box& b, std::size_t max_calls = 100000, + Real tolerance = tools::root_epsilon(), Real* error = nullptr, Real* L1 = nullptr) + -> decltype(f(x)); + + }}} // namespaces + +[heading Description] + +`sparse_grid` integrates /f/ over \[-1, 1\][super /d/] with `leja_legendre`, or /f/(/x/) exp(-|/x/|[super 2]) over the whole space +with `leja_hermite`, or against the product of the weights of `leja_laguerre` or `leja_jacobi`, by the dimension-adaptive sparse grid of Gerstner and Griebel. Coordinate /i/ at level /k[sub i]/ uses the +first 2/k[sub i]/ - 1 points of the rule, and the grid is the sum of the tensor products of differences of consecutive levels +over a set of multi-indices /k/. The set grows where these contributions are largest, so the grid spends its points on the +coordinates, and the interactions between them, that matter. Exact integrals of products of polynomials of degree 2/k[sub i]/ - 2 +need only the multi-index /k/ and those below it, and a function of a few coordinates costs little more than in those dimensions. + +The levels add two points at a time, not one: for an even /f/ such as /x/[super 2] the rules with one and two points agree, and a +zero difference would stop the refinement of that coordinate. For the same reason, when the error estimate (the sum of the +norms of the contributions not yet refined) falls below `tolerance` * /L1/, every such contribution is refined once more, +and the result is accepted only if the estimate is still below the tolerance. Nothing is accepted while /L1/ is 0, as it is while /f/ has been +zero at every point: /x/[sub 1][super 2] /x/[sub 2][super 2], with `leja_legendre`, is zero at the first point, the mean of the weight, +and on the axes through it. After /d/ + 1 such rounds of refinement with /f/ still zero everywhere, the result is 0 with an error estimate of 0. An interaction of three or more variables +that vanishes on the planes through the mean can still be missed if /f/ is not zero there, such as 1 + /x/[sub 1][super 2] /x/[sub 2][super 2] /x/[sub 3][super 2]. /L1/ is the grid's estimate of the integral of |/f/| times the +weight, the scale of the tolerance; before the grid has converged it can be far off, even negative. It is not a condition +number: the sparse grid's weights have mixed signs, and the sum of |/w[sub i] f/(/x[sub i]/)| can be much larger. +If `error` or `L1` are not null they are set to these. The grid computes the points of `rule` it needs, and keeps them there for the next call. + +`max_calls` bounds the number of evaluations of /f/: a refinement that would exceed it is not made, and the result has an +error estimate above the tolerance. A contribution that would need more points than that in a coordinate is not refined +further, and its norm stays in the error estimate; +if the calls run out during the final refinements the error estimate is NaN. No error is raised. A __domain_error is +raised if `dimension` or `max_calls` is zero, or if /f/ takes a fixed-size container whose size is not `dimension`. /f/ may return a real or complex value or a vector of them, as for the +one-dimensional rules; a vector is integrated componentwise, and its error estimate and /L1/ use the max norm. +With bounds `a` and `b`, each coordinate is mapped as `integrate(f, a[i], b[i])` maps it, and the result, error and /L1/ are scaled by the +product of the factors; a __domain_error is raised if a bound is infinite or NaN or `b.size()` is not `a.size()`. + +The bookkeeping costs under 2 microseconds per evaluation of /f/ in double precision on Genz's test integrands, more where +the levels are high, so a sparse grid pays off when /f/ costs more than that. On Genz's test integrands, with 10[super 4] calls or more, it beats +quasi-Monte Carlo (Boost.Random's Sobol' sequence) by orders of magnitude up to four dimensions; in eight it still wins for +entire functions such as cos(0.7 + [sum] /c[sub i] x[sub i]/), but not for peaked ones such as Genz's product peak and Gaussian, +and in sixteen quasi-Monte Carlo is better unless /f/ depends mostly on a few coordinates. For integrands with kinks, such as +exp(-[sum] /c[sub i]/|/x[sub i]/ - /u[sub i]/|), quasi-Monte Carlo is better in any dimension, and the error estimate, a heuristic +as for the one-dimensional rule, can be too small. + +[heading References] + +* T. Gerstner and M. Griebel, ['Dimension-adaptive tensor-product quadrature], Computing 71 (2003), 65-87, [@https://doi.org/10.1007/s00607-003-0015-5 doi:10.1007/s00607-003-0015-5]. +* A. Narayan and J. D. Jakeman, as above. + +[endsect] [/section:sparse_grid Sparse grids] diff --git a/example/Jamfile.v2 b/example/Jamfile.v2 index 79a87b2713..ea8e32a9a8 100644 --- a/example/Jamfile.v2 +++ b/example/Jamfile.v2 @@ -151,6 +151,7 @@ test-suite examples : [ run reciprocal_fibonacci_constant.cpp : : : [ check-target-builds ../config//has_mpfr "MPFR Support" : "-lmpfr -lgmp" : no ] ] [ run reverse_mode_linear_regression_example.cpp ] [ run autodiff_reverse_black_scholes.cpp ] + [ compile leja_parallel_assembly.cpp : [ requires cxx17_if_constexpr cxx11_hdr_thread ] [ check-target-builds ../../multiprecision/config//has_eigen : : no ] ] ; run root_elliptic_finding.cpp /boost/timer : : : release static [ requires cxx11_unified_initialization_syntax cxx11_defaulted_functions ] freebsd:"-lrt" linux:"-lrt -lpthread" ; diff --git a/example/leja_parallel_assembly.cpp b/example/leja_parallel_assembly.cpp new file mode 100644 index 0000000000..ece8a14bfb --- /dev/null +++ b/example/leja_parallel_assembly.cpp @@ -0,0 +1,111 @@ +// Copyright Nick Thompson 2026. +// Use, modification and distribution are subject to the +// Boost Software License, Version 1.0. (See accompanying file +// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) + +// Parallel assembly of a Hamiltonian with one shared Leja rule. +// +// The anharmonic oscillator H = -d^2/dx^2 + x^4, in the basis of Hermite functions psi_n(x / s), n < N. With y = x / s, +// H_ij = int psi_i(y) ((2j + 1 - y^2) / s^2 + s^4 y^4) psi_j(y) dy, since -psi_j'' = (2j + 1 - y^2) psi_j. The integrands +// are entire and negligible beyond |y| = L, and the threads integrate them over [-L, L] with one leja_legendre: whichever +// first needs more points computes them, and the others use them. Each entry takes as many points as it needs, and gets +// an error estimate from the points it has already used. +// +// The matrix elements are also known exactly, from x = (a + a^+) / sqrt(2), and the lowest eigenvalues agree with those +// of Hioe and Montroll (1975). +// +// A rule for exp(-y^2), leja_hermite, would be the wrong choice: h_i h_j grows like exp(y^2), the Leja points spread +// further than the Gauss points, and sum |w_k h_i h_j(x_k)| is 1e19 times H_ij for i = j = 40. + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +// psi_i(y) psi_j(y), from h_{n+1} = y sqrt(2 / (n + 1)) h_n - sqrt(n / (n + 1)) h_{n-1} times exp(-y^2 / 2). +double psi_product(unsigned i, unsigned j, double y) +{ + double previous = 0; + double psi = std::exp(-y * y / 2) / std::sqrt(std::sqrt(boost::math::constants::pi())); + double psi_i = i == 0 ? psi : 0; + double psi_j = j == 0 ? psi : 0; + for (unsigned n = 0; n < (std::max)(i, j); ++n) + { + const double next = y * std::sqrt(2.0 / (n + 1)) * psi - std::sqrt(double(n) / (n + 1)) * previous; + previous = psi; + psi = next; + if (n + 1 == i) + psi_i = psi; + if (n + 1 == j) + psi_j = psi; + } + return psi_i * psi_j; +} + +// H from the matrix X of y, tridiagonal with X_{n,n+1} = sqrt((n + 1) / 2): X^2 and X^4 are exact in the first N rows +// when X has N + 4. +Eigen::MatrixXd exact(unsigned N, double s) +{ + Eigen::MatrixXd X = Eigen::MatrixXd::Zero(N + 4, N + 4); + for (unsigned n = 0; n < N + 3; ++n) + X(n, n + 1) = X(n + 1, n) = std::sqrt((n + 1) / 2.0); + const Eigen::MatrixXd X2 = X * X; + const Eigen::MatrixXd D = Eigen::VectorXd::LinSpaced(N, 1, 2 * N - 1).asDiagonal(); + return (D - X2.topLeftCorner(N, N)) / (s * s) + s * s * s * s * (X2 * X2).topLeftCorner(N, N); +} + +int main() +{ + const unsigned N = 60; + const double s = 0.75; + const double L = 14; + const unsigned threads = (std::max)(1u, std::thread::hardware_concurrency()); + const boost::math::quadrature::leja_legendre rule; + + Eigen::MatrixXd H(N, N); + Eigen::MatrixXi calls(N, N); + // The entries i <= j, taken in turn by the threads. + auto assemble = [&] { + std::atomic next{0}; + auto work = [&] { + for (unsigned e = next++; e < N * N; e = next++) + { + const unsigned i = e / N, j = e % N; + if (i > j) + continue; + unsigned count = 0; + H(i, j) = H(j, i) = rule.integrate([&](double y) { + ++count; + return psi_product(i, j, y) * ((2 * j + 1 - y * y) / (s * s) + s * s * s * s * y * y * y * y); + }, -L, L, 1e-13); + calls(i, j) = calls(j, i) = count; + } + }; + const auto start = std::chrono::steady_clock::now(); + std::vector pool; + for (unsigned t = 0; t < threads; ++t) + pool.emplace_back(work); + for (std::thread& t : pool) + t.join(); + return std::chrono::duration(std::chrono::steady_clock::now() - start).count(); + }; + const double cold = assemble(); + const double warm = assemble(); + + const Eigen::MatrixXd H_exact = exact(N, s); + const double worst = (H - H_exact).cwiseAbs().maxCoeff(), largest = H_exact.cwiseAbs().maxCoeff(); + std::cout << std::setprecision(17) << N * (N + 1) / 2 << " entries on " << threads << " threads, with " + << calls.minCoeff() << " to " << calls.maxCoeff() + << " points: " << cold << " s computing the rule as it goes, " << warm << " s with it computed.\n" + << "max |H - H_exact| = " << worst << ", " << worst / largest << " of the largest entry.\n"; + const Eigen::VectorXd E = Eigen::SelfAdjointEigenSolver(H, Eigen::EigenvaluesOnly).eigenvalues(); + const double reference[] = {1.0603620904841829, 3.7996730298013941, 7.4556979379867383, 11.644745511378162}; + for (unsigned n = 0; n < 4; ++n) + std::cout << "E_" << n << " = " << E(n) << ", Hioe and Montroll " << reference[n] << '\n'; +} diff --git a/include/boost/math/quadrature/leja.hpp b/include/boost/math/quadrature/leja.hpp new file mode 100644 index 0000000000..b5963ccbdf --- /dev/null +++ b/include/boost/math/quadrature/leja.hpp @@ -0,0 +1,727 @@ +// Copyright Nick Thompson 2026. +// Use, modification and distribution are subject to the +// Boost Software License, Version 1.0. (See accompanying file +// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) + +#ifndef BOOST_MATH_QUADRATURE_LEJA_HPP +#define BOOST_MATH_QUADRATURE_LEJA_HPP + +#ifdef _MSC_VER +#pragma once +#endif + +#include + +#ifdef BOOST_MATH_NO_CXX17_IF_CONSTEXPR +#error "The header requires C++17 or later." +#endif + +#ifndef BOOST_MATH_BUILD_MODULE +#include +#include +#include +#include +#include +#include +#include +#ifdef BOOST_MATH_HAS_THREADS +#include +#include +#endif +#endif +#include +#include +#include +#include +#include +#include +#include +#include + +BOOST_MATH_NAMESPACE_BEGIN namespace quadrature { namespace detail { + +// The most points a rule computes. It keeps every rule, O(N^2) numbers, and takes O(N^2) operations to compute them: +// about 60 ms for 1023 in double precision, and 0.4 s for 127 in cpp_bin_float_50. +template +constexpr unsigned leja_max_points = std::is_floating_point::value ? 1023 : 127; + +// hi + lo, to twice the precision of Real, for the few quantities the weights cannot afford to round. +template +using leja_double = BOOST_MATH_NAMESPACE::detail::double_word; +template +using leja_arithmetic = BOOST_MATH_NAMESPACE::detail::double_word_arithmetic; + +template +leja_double leja_ldexp(const leja_double& a, int e) +{ + using std::ldexp; + return {ldexp(a.hi, e), ldexp(a.lo, e)}; +} + +// A weight w on [lower(), upper()] is described by the total mass mu0, the coefficients of the recurrence +// x pi_k = pi_{k+1} + alpha_k pi_k + beta_k^2 pi_{k-1} of its monic orthogonal polynomials, to twice the precision of Real, +// and log sqrt(u) and its derivatives, where u places the points. u = w, except that a negative exponent of w at an end of the support, where +// sqrt(w) has no maximum, is 0 in u. Every exponent must be finite and > -1; exponent() is one that is not, if there is one. +// The points start from the mean of w, the one-point Gauss rule, not from 0, where x and x^2 vanish and a sparse grid +// of x_1 x_2 would see only zeros. +template +struct leja_legendre_weight +{ + static constexpr const char* function = "boost::math::quadrature::leja_legendre<%1%>"; + static Real lower() { return -1; } + static Real upper() { return 1; } + static Real exponent() { return 0; } + static std::vector start() { return {Real(0), Real(-1), Real(1)}; } + static Real mu0() { return 2; } + static Real half_log(const Real&) { return 0; } + static Real half_log_derivative(const Real&) { return 0; } + static Real half_log_second_derivative(const Real&) { return 0; } + static leja_double alpha(unsigned, const leja_arithmetic&) { return {0, 0}; } + static leja_double beta_squared(unsigned k, const leja_arithmetic& op) + { + return op.divide(leja_double{Real(k) * k, 0}, leja_double{Real(4) * k * k - 1, 0}); + } + // The factor h |h|^(alpha + beta) by which x = c + h t takes the integral over [-1, 1] to one over [c - h, c + h]. + static Real jacobian(const Real& h) { return h; } +}; + +template +struct leja_hermite_weight +{ + static constexpr const char* function = "boost::math::quadrature::leja_hermite<%1%>"; + static Real lower() { return -std::numeric_limits::infinity(); } + static Real upper() { return std::numeric_limits::infinity(); } + static Real exponent() { return 0; } + static std::vector start() { return {Real(0)}; } + static Real mu0() { return BOOST_MATH_NAMESPACE::constants::root_pi(); } + static Real half_log(const Real& x) { return -x * x / 2; } + static Real half_log_derivative(const Real& x) { return -x; } + static Real half_log_second_derivative(const Real&) { return -1; } + static leja_double alpha(unsigned, const leja_arithmetic&) { return {0, 0}; } + static leja_double beta_squared(unsigned k, const leja_arithmetic&) { return {Real(k) / 2, 0}; } +}; + +// x^a exp(-x) on [0, infinity). +template +struct leja_laguerre_weight +{ + static constexpr const char* function = "boost::math::quadrature::leja_laguerre<%1%>"; + explicit leja_laguerre_weight(const Real& alpha_ = 0) : a(alpha_), u((std::max)(alpha_, Real(0))) {} + static Real lower() { return 0; } + static Real upper() { return std::numeric_limits::infinity(); } + Real exponent() const { return a; } + std::vector start() const { return {a + 1}; } + Real mu0() const { return BOOST_MATH_NAMESPACE::detail::laguerre_family(0, a).area(); } + Real half_log(const Real& x) const + { + using std::log; + return (u > 0 ? u * log(x) / 2 : Real(0)) - x / 2; + } + Real half_log_derivative(const Real& x) const { return (u > 0 ? u / (2 * x) : Real(0)) - Real(1) / 2; } + Real half_log_second_derivative(const Real& x) const { return u > 0 ? -u / (2 * x * x) : Real(0); } + leja_double alpha(unsigned k, const leja_arithmetic& op) const { return op.add(D{Real(2 * k + 1), 0}, a); } + leja_double beta_squared(unsigned k, const leja_arithmetic& op) const { return op.multiply(op.add(D{Real(k), 0}, a), Real(k)); } + using D = leja_double; + + Real a, u; +}; + +// (1 - x)^a (1 + x)^b on [-1, 1]. +template +struct leja_jacobi_weight +{ + static constexpr const char* function = "boost::math::quadrature::leja_jacobi<%1%>"; + leja_jacobi_weight(const Real& alpha_, const Real& beta_) + : a(alpha_), b(beta_), ua((std::max)(alpha_, Real(0))), ub((std::max)(beta_, Real(0))) {} + static Real lower() { return -1; } + static Real upper() { return 1; } + Real exponent() const { return a > -1 && (BOOST_MATH_NAMESPACE::isfinite)(a) ? b : a; } + Real jacobian(const Real& h) const + { + using std::abs; + using std::pow; + return h * pow(abs(h), a + b); + } + std::vector start() const { return {(b - a) / (a + b + 2)}; } + Real mu0() const { return BOOST_MATH_NAMESPACE::detail::jacobi_family(0, a, b).area(); } + Real half_log(const Real& x) const + { + using std::log; + return (ua > 0 ? ua * log(1 - x) / 2 : Real(0)) + (ub > 0 ? ub * log(1 + x) / 2 : Real(0)); + } + Real half_log_derivative(const Real& x) const + { + return (ub > 0 ? ub / (2 * (1 + x)) : Real(0)) - (ua > 0 ? ua / (2 * (1 - x)) : Real(0)); + } + Real half_log_second_derivative(const Real& x) const + { + return -(ua > 0 ? ua / (2 * (1 - x) * (1 - x)) : Real(0)) - (ub > 0 ? ub / (2 * (1 + x) * (1 + x)) : Real(0)); + } + // https://dlmf.nist.gov/18.9; k = 0 and 1 apart, where 2k + a + b or 2k + a + b - 1 can vanish. + leja_double alpha(unsigned k, const leja_arithmetic& op) const + { + const D b_minus_a = op.add(D{b, 0}, -a), s = op.add(op.add(D{Real(2 * k), 0}, a), b), s2 = op.add(s, Real(2)); + if (k == 0) + return op.divide(b_minus_a, s2); + return op.divide(op.multiply(b_minus_a, op.add(D{b, 0}, a)), op.multiply(s, s2)); + } + leja_double beta_squared(unsigned k, const leja_arithmetic& op) const + { + const D s = op.add(op.add(D{Real(2 * k), 0}, a), b), s_s1 = op.multiply(op.multiply(s, s), op.add(s, Real(1))); + if (k == 1) + return op.divide(op.multiply(op.multiply(op.add(D{1, 0}, a), Real(4)), op.add(D{1, 0}, b)), s_s1); + const D ka = op.add(D{Real(k), 0}, a), kb = op.add(D{Real(k), 0}, b); + const D numerator = op.multiply(op.multiply(op.multiply(ka, Real(4 * k)), kb), op.add(ka, b)); + return op.divide(numerator, op.multiply(s_s1, op.add(s, Real(-1)))); + } + using D = leja_double; + + Real a, b, ua, ub; +}; + +// A vector-valued integrand returns a random-access container of reals or complex numbers, such as std::array or std::vector. +template +struct leja_is_vector : std::false_type {}; +template +struct leja_is_vector()[0]), decltype(std::declval().size())>> : std::true_type {}; + +// |a - b|, or |a| without b, in the max norm for a vector; NaN if any component is NaN. +template +Real leja_distance(const K& a, const K* b = nullptr) +{ + using std::abs; + if constexpr (leja_is_vector::value) + { + Real m = 0; + for (std::size_t j = 0; j < a.size(); ++j) + { + const Real d = b ? Real(abs(a[j] - (*b)[j])) : Real(abs(a[j])); + if (d != d) + return d; + m = (std::max)(m, d); + } + return m; + } + else + return b ? Real(abs(a - *b)) : Real(abs(a)); +} + +// NaN, in every component of a default-constructed vector (none for std::vector). +template +K leja_nan() +{ + if constexpr (leja_is_vector::value) + { + K y{}; + for (std::size_t j = 0; j < y.size(); ++j) + y[j] = std::numeric_limits::quiet_NaN(); + return y; + } + else + return K(std::numeric_limits::quiet_NaN()); +} + +// y = a x if first, else y += a x; a vector y takes its size from x. +template +void leja_accumulate(K& y, const Real& a, const K& x, bool first) +{ + if constexpr (leja_is_vector::value) + { + if (first) + y = x; + for (std::size_t j = 0; j < y.size(); ++j) + { + if (first) + y[j] = a * x[j]; + else + y[j] += a * x[j]; + } + } + else if (first) + y = a * x; + else + y += a * x; +} + +template +class leja_detail +{ +public: + using value_type = Real; + using policy_type = Policy; + + // The arguments are the weight's: alpha for leja_laguerre, alpha and beta for leja_jacobi. The points and weights + // are computed as they are needed, by any thread, up to max_points(). + template ::value, bool>::type = true> + explicit leja_detail(const Args&... args) : m_weight(args...) + { + reserve(); + if (!(m_weight.exponent() > -1 && (BOOST_MATH_NAMESPACE::isfinite)(m_weight.exponent()))) + policies::raise_domain_error(Weight::function, "The exponents of the weight must be finite and > -1, but got %1%.", m_weight.exponent(), Policy()); + } + + leja_detail(const leja_detail& rule) : m_weight(rule.m_weight) + { +#ifdef BOOST_MATH_HAS_THREADS + std::lock_guard lock(rule.m_mutex); +#endif + reserve(); + const unsigned size = rule.m_size; + m_abscissa.assign(rule.m_abscissa.begin(), rule.m_abscissa.begin() + size); + m_weights.assign(rule.m_weights.begin(), rule.m_weights.begin() + size + 1); + m_state = rule.m_state; + m_size = size; + } + + leja_detail& operator=(const leja_detail& rule) + { + leja_detail copy(rule); + m_weight = copy.m_weight; + m_abscissa.swap(copy.m_abscissa); + m_weights.swap(copy.m_weights); + std::swap(m_state, copy.m_state); + m_size = static_cast(copy.m_size); + return *this; + } + + // The points computed so far, in Leja order; the n-point rule uses the first n, which weights(n) computes. Those + // another thread is computing are NaN until it is done. + const std::vector& abscissa() const + { + return m_abscissa; + } + + const std::vector& weights(unsigned n) const + { + if (n == 0 || n > max_points()) + { + policies::raise_domain_error(Weight::function, "The number of points must be in [1, max_points()], but got %1%.", Real(n), Policy()); + return m_weights[0]; + } + if (n > m_size) + extend(n); + // Empty for a weight with an invalid exponent, after a domain error the policy ignored. + return n <= m_size ? m_weights[n] : m_weights[0]; + } + + static constexpr unsigned max_points() + { + return leja_max_points; + } + + // Adds one point at a time until |Q_m - Q_{m - max(1, m/4)}| is small for every m in [n - n/4, n]: a single difference, + // and |Q_n - Q_{n-1}| even more so, often underestimates the error of slowly converging integrands. Like tanh_sinh, + // it returns its best result if max_points() is exhausted first. + template + auto integrate(F f, Real tolerance = tools::root_epsilon(), Real* error = nullptr, Real* pL1 = nullptr) const + ->decltype(std::declval()(std::declval())) + { + using K = decltype(f(Real(0))); + static_assert(!std::is_integral::value, + "The return type cannot be integral, it must be a real or complex floating point type, or a vector of them."); + std::vector y, Q; + // NaN for an empty rule. + K result = leja_nan(); + Real L1 = std::numeric_limits::quiet_NaN(); + Real err = L1; + std::vector e; // e[n - 1] = |Q_n - Q_{n - max(1, n/4)}| + for (unsigned n = 1; n <= max_points() && !weights(n).empty(); ++n) + { + y.push_back(f(m_abscissa[n - 1])); + result = sum(y, n, L1); + Q.push_back(result); + e.push_back(n > 1 ? leja_distance(result, &Q[n - 1 - (std::max)(1u, n / 4)]) : L1); + err = *std::max_element(e.end() - (n / 4 + 1), e.end()); + // A NaN from f: max_element need not return it. + if (L1 != L1) + { + err = L1; + break; + } + if (n >= 4 && err <= tolerance * L1) + break; + } + if (error) + *error = err; + if (pL1) + *pL1 = L1; + return result; + } + + // The integral over [a, b], through x = (a + b)/2 + h t, h = (b - a)/2: of f for leja_legendre, and of + // f(x) |b - x|^alpha |x - a|^beta for leja_jacobi. For a > b it is minus the integral over [b, a]. + template ().jacobian(Real(0)))> + auto integrate(F f, const Real& a, const Real& b, Real tolerance = tools::root_epsilon(), Real* error = nullptr, + Real* pL1 = nullptr) const->decltype(std::declval()(std::declval())) + { + using std::abs; + using K = decltype(f(Real(0))); + const Real h = (b - a) / 2; + if (!((BOOST_MATH_NAMESPACE::isfinite)(a) && (BOOST_MATH_NAMESPACE::isfinite)(b)) || h == 0) + { + const Real r = h == 0 ? Real(0) : policies::raise_domain_error(Weight::function, + "The bounds must be finite, but got %1%; leja_hermite and leja_laguerre integrate over infinite intervals.", + (BOOST_MATH_NAMESPACE::isfinite)(a) ? b : a, Policy()); + if (error) + *error = r; + if (pL1) + *pL1 = r; + // 0, of the shape of f(a), over an empty interval. + K result = leja_nan(); + if (h == 0) + leja_accumulate(result, Real(0), f(a), true); + return result; + } + const Real c = (a + b) / 2, J = jacobian(h); + K result; + leja_accumulate(result, J, integrate([&](const Real& t) { return f(c + h * t); }, tolerance, error, pL1), true); + if (error) + *error *= abs(J); + if (pL1) + *pL1 *= abs(J); + return result; + } + + // The factor by which the map from [-1, 1] to an interval of half-width h, of the sign of its orientation, scales an + // integral; leja_legendre and leja_jacobi only. + template + auto jacobian(const Real& h) const->decltype(std::declval().jacobian(h)) + { + return m_weight.jacobian(h); + } + +private: + void reserve() + { + m_abscissa.reserve(max_points()); + m_weights.reserve(max_points() + 1); + m_weights.resize(1); + } + + // Computes the points and weights of the rules up to n points, which the other threads see once m_size is n. The + // vectors never reallocate, so the references the others hold stay valid. + void extend(unsigned n) const + { +#ifdef BOOST_MATH_HAS_THREADS + std::lock_guard lock(m_mutex); +#endif + const unsigned size = m_size; + if (n <= size || !(m_weight.exponent() > -1 && (BOOST_MATH_NAMESPACE::isfinite)(m_weight.exponent()))) + return; + std::vector x(m_abscissa.begin(), m_abscissa.begin() + size); + m_abscissa.resize(n, std::numeric_limits::quiet_NaN()); + compute_abscissa(x, n); + std::copy(x.begin() + size, x.end(), m_abscissa.begin() + size); + compute_weights(size, n); + m_size = n; + } + + template + K sum(const std::vector& y, unsigned n, Real& L1) const + { + using std::abs; + K result{}; + L1 = 0; + for (unsigned i = 0; i < n; ++i) + { + leja_accumulate(result, m_weights[n][i], y[i], i == 0); + L1 += abs(m_weights[n][i]) * leja_distance(y[i]); + } + return result; + } + + // F(t) = log sqrt(u(t)) + sum_k log|t - x_k| is concave between consecutive points, so each gap has one maximum, + // found by Newton's method on F' safeguarded by bisection. + Real derivatives(const std::vector& points, const Real& t, Real& second) const + { + Real first = m_weight.half_log_derivative(t); + second = m_weight.half_log_second_derivative(t); + for (const Real& x : points) + { + const Real r = 1 / (t - x); + first += r; + second -= r * r; + } + return first; + } + + // The points between l and r, or an end of the support where end_l or end_r. lower = F(t) is at most the maximum of + // F there, and upper at least; they are equal, at the maximum t, after a search. + struct gap + { + Real l, r, t, lower, upper, scale; + bool end_l, end_r; + }; + + // Sets t to the maximum of F on the gap, starting from t0. + void search(const std::vector& points, gap& g, Real t) const + { + using std::abs; + using std::frexp; + using std::log; + const Real inf = std::numeric_limits::infinity(); + Real l = g.l, r = g.r, second; + // At a finite end of the support, u > 0 can make the end itself the maximum. + const bool at_l = g.end_l && l > -inf && derivatives(points, l, second) <= 0; + if (at_l || (g.end_r && r < inf && derivatives(points, r, second) >= 0)) + t = at_l ? l : r; + else + { + if (l == -inf || r == inf) + { // Double the ray's length until F' changes sign. + const Real end = l == -inf ? r : l; + const Real direction = l == -inf ? -1 : 1; + Real step = 1; + while (direction * derivatives(points, end + direction * step, second) > 0) + step *= 2; + l = end + direction * step; + r = end + direction * step / 2; + if (step == 1) + r = end; + if (direction > 0) + std::swap(l, r); + } + // Relative to the gap, not to t, which can be 0. + const Real tolerance = 2 * tools::epsilon() * (std::max)(abs(l), abs(r)); + if (!(t > l && t < r)) + t = (l + r) / 2; + for (unsigned i = 0; i < 4 * tools::digits(); ++i) + { + const Real first = derivatives(points, t, second); + (first > 0 ? l : r) = t; + Real next = t - first / second; + if (abs(next - t) <= tolerance) + { + t = next; + break; + } + if (!(next > l && next < r)) + next = (l + r) / 2; + t = next; + if (r - l <= 2 * tolerance) + break; + } + } + // F(t), with one logarithm: the product of the |t - x_k| is P 2^e, renormalized every four factors. Each factor + // adds at most eps to its rounding error, and log rounds |log P + e ln 2| and |log sqrt(u)| with relative error eps. + Real P = 1; + int e = 0; + for (std::size_t k = 0; k < points.size(); k += 4) + { + Real block = abs(t - points[k]); + for (std::size_t j = k + 1; j < (std::min)(k + 4, points.size()); ++j) + block *= abs(t - points[j]); + int exponent; + P = frexp(P * block, &exponent); + e += exponent; + } + const Real half_log = m_weight.half_log(t); + const Real product = log(P) + e * constants::ln_two(); + g.t = t; + g.lower = g.upper = half_log + product; + g.scale = abs(half_log) + abs(product) + 2 * points.size(); + } + + // Adds points to points until it has N. A new point x adds log|t - x| to F, so a gap's maximum grows by at most the + // largest log|t - x| over the gap, and by at least log|t - x| at the last maximum: only the gaps whose upper bound + // reaches the largest lower bound need a search, about 8 of them for each point up to 1023. + void compute_abscissa(std::vector& points, unsigned N) const + { + using std::abs; + using std::log; + const Real eps = tools::epsilon(); + const Real inf = std::numeric_limits::infinity(); + const std::vector start = m_weight.start(); + while (points.size() < (std::min)(static_cast(N), start.size())) + points.push_back(start[points.size()]); + std::vector& gaps = m_state.gaps; + if (points.size() < N && gaps.empty()) + { + std::vector sorted = points; + std::sort(sorted.begin(), sorted.end()); + for (std::size_t g = 0; g <= sorted.size(); ++g) + { + gap G{g == 0 ? m_weight.lower() : sorted[g - 1], g == sorted.size() ? m_weight.upper() : sorted[g], 0, 0, 0, 0, + g == 0, g == sorted.size()}; + if (G.l < G.r) + { + search(points, G, (G.l + G.r) / 2); + gaps.push_back(G); + } + } + } + while (points.size() < N) + { + // The bounds gain a rounding error with each point. + Real best = -inf, scale = 0; + for (const gap& g : gaps) + { + best = (std::max)(best, g.lower); + scale = (std::max)(scale, g.scale); + } + const Real slack = 4 * (points.size() + 8) * eps * (abs(best) + scale); + std::vector searched; + for (std::size_t i = 0; i < gaps.size(); ++i) + { + if (gaps[i].upper < best - slack) + continue; + if (gaps[i].upper > gaps[i].lower) + search(points, gaps[i], gaps[i].t); + searched.push_back(i); + } + // Maxima within 8 eps scale, a few times the rounding error of F, are tied, and the leftmost wins. Exact ties + // are common (+-1/sqrt(3) for Legendre, +-1 for Hermite), and this way every precision picks the same point. + std::size_t top = searched[0]; + for (std::size_t i : searched) + if (gaps[i].lower > gaps[top].lower) + top = i; + std::size_t c = 0; + while (std::find(searched.begin(), searched.end(), c) == searched.end() + || gaps[c].lower < gaps[top].lower - 8 * eps * gaps[top].scale) + ++c; + const Real x = gaps[c].t; + for (gap& g : gaps) + { + g.lower += log(abs(g.t - x)); + g.upper += g.l == -inf || g.r == inf ? inf : log((std::max)(abs(g.l - x), abs(g.r - x))); + } + points.push_back(x); + gap left{gaps[c].l, x, 0, 0, 0, 0, gaps[c].end_l, false}, right{x, gaps[c].r, 0, 0, 0, 0, false, gaps[c].end_r}; + gaps.erase(gaps.begin() + c); + for (gap* g : {&right, &left}) + { + if (g->l < g->r) + { + search(points, *g, (g->l + g->r) / 2); + gaps.insert(gaps.begin() + c, *g); + } + } + } + } + + // The rules in Newton form: Q_{n+1} f = Q_n f + I[omega_n] f[x_0, ..., x_n], where omega_n = prod_{k; + const leja_arithmetic op; + state& f = m_state; + const Real mu0 = m_weight.mu0(); + auto normalize = [](D& v, int& e) { + int k; + frexp(v.hi, &k); + v = leja_ldexp(v, -k); + e += k; + }; + for (unsigned n = old; n < N; ++n) + { + const Real& x = m_abscissa[n]; + if (n == 0) + { + f.c = {D{1, 0}}; + f.e = {0}; + f.d = {D{1, 0}}; + f.d_exponent = {0}; + f.w = {D{mu0, 0}}; + m_weights.push_back({mu0}); + continue; + } + // omega_n = (x - x_{n-1}) omega_{n-1}. + f.log2_norm += log(m_weight.beta_squared(n, op).hi) / (2 * constants::ln_two()); + f.e.push_back(static_cast(floor(f.log2_norm + Real(0.5)))); + std::vector c(n + 1, D{0, 0}); + Real largest = 0; + for (unsigned k = 0; k <= n; ++k) + { + if (k > 0) + c[k] = leja_ldexp(f.c[k - 1], f.e[k] - f.e[k - 1]); + if (k < n) + c[k] = op.add(c[k], op.multiply(op.add(m_weight.alpha(k, op), -m_abscissa[n - 1]), f.c[k])); + if (k + 1 < n) + c[k] = op.add(c[k], leja_ldexp(op.multiply(m_weight.beta_squared(k + 1, op), f.c[k + 1]), f.e[k] - f.e[k + 1])); + largest = (std::max)(largest, Real(abs(c[k].hi))); + } + int e; + frexp(largest, &e); + for (D& v : c) + v = leja_ldexp(v, -e); + f.c_exponent += e; + f.c.swap(c); + // d_i gains the factor x_i - x, and d_n is new. + D d{1, 0}; + int d_exponent = 0; + for (unsigned i = 0; i < n; ++i) + { + f.d[i] = op.multiply(f.d[i], op.two_sum(m_abscissa[i], -x)); + normalize(f.d[i], f.d_exponent[i]); + d = op.multiply(d, op.two_sum(x, -m_abscissa[i])); + normalize(d, d_exponent); + } + f.d.push_back(d); + f.d_exponent.push_back(d_exponent); + f.w.push_back(D{0, 0}); + const D I = op.multiply(f.c[0], mu0); + std::vector w(n + 1); + for (unsigned i = 0; i <= n; ++i) + { + f.w[i] = op.add(f.w[i], leja_ldexp(op.divide(I, f.d[i]), f.c_exponent - f.d_exponent[i])); + w[i] = f.w[i].hi + f.w[i].lo; + } + m_weights.push_back(std::move(w)); + } + } + + // The Newton form of the rules so far, and the gaps between the points, so that a rule computes the same points and + // weights however it grows. + struct state + { + std::vector> c, d, w; + std::vector e, d_exponent; + int c_exponent = 0; + Real log2_norm = 0; + std::vector gaps; + }; + + Weight m_weight; + mutable std::vector m_abscissa; + mutable std::vector> m_weights; // m_weights[n] is the n-point rule; m_weights[0] is empty. + mutable state m_state; +#ifdef BOOST_MATH_HAS_THREADS + mutable std::mutex m_mutex; + mutable std::atomic m_size{0}; +#else + mutable unsigned m_size = 0; +#endif +}; + +} // namespace detail + +// Nested Leja quadrature for the integral of f over [-1, 1], the Legendre weight w = 1. +template > +using leja_legendre = detail::leja_detail, Policy>; + +// Nested weighted Leja quadrature for the integral of f(x) exp(-x^2) over the real line. +template > +using leja_hermite = detail::leja_detail, Policy>; + +// Nested weighted Leja quadrature for the integral of f(x) x^alpha exp(-x) over [0, infinity), alpha > -1: +// leja_laguerre(alpha = 0). +template > +using leja_laguerre = detail::leja_detail, Policy>; + +// Nested weighted Leja quadrature for the integral of f(x) (1 - x)^alpha (1 + x)^beta over [-1, 1], alpha, beta > -1: +// leja_jacobi(alpha, beta). +template > +using leja_jacobi = detail::leja_detail, Policy>; + +} // namespace quadrature +BOOST_MATH_NAMESPACE_END + +#endif // BOOST_MATH_QUADRATURE_LEJA_HPP diff --git a/include/boost/math/quadrature/sparse_grid.hpp b/include/boost/math/quadrature/sparse_grid.hpp new file mode 100644 index 0000000000..ecb99f3ef5 --- /dev/null +++ b/include/boost/math/quadrature/sparse_grid.hpp @@ -0,0 +1,353 @@ +// Copyright Nick Thompson 2026. +// Use, modification and distribution are subject to the +// Boost Software License, Version 1.0. (See accompanying file +// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) + +#ifndef BOOST_MATH_QUADRATURE_SPARSE_GRID_HPP +#define BOOST_MATH_QUADRATURE_SPARSE_GRID_HPP + +#ifdef _MSC_VER +#pragma once +#endif + +#include + +#ifndef BOOST_MATH_BUILD_MODULE +#include +#include +#include +#include +#include +#include +#include +#endif + +BOOST_MATH_NAMESPACE_BEGIN namespace quadrature { + +namespace detail { +// The point f is called with: std::vector if f takes one, otherwise the argument of f's only call operator, such as +// std::array or Eigen::VectorXd. +template +std::decay_t leja_argument(std::function); +template +struct leja_point { using type = std::vector; }; +template +struct leja_point&>, + std::void_t()}))>>> +{ + using type = decltype(leja_argument(std::function{std::declval()})); +}; + +template +struct leja_resizable : std::false_type {}; +template +struct leja_resizable().resize(0u))>> : std::true_type {}; + +// A rule on [-1, 1], which integrate(f, a, b) maps to [a, b]. +template +struct leja_bounded : std::false_type {}; +template +struct leja_bounded().jacobian(typename Rule::value_type(1)))>> : std::true_type {}; +} + +// Dimension-adaptive sparse grid (Gerstner and Griebel) for the integral of f over the product of the rule's weight, +// as [-1, 1]^d for leja_legendre and exp(-|x|^2) over R^d for leja_hermite. Level k of a coordinate is the 2k - 1 point +// Leja rule U_k, and the multi-index k contributes Delta_k f = (U_{k_1} - U_{k_1 - 1}) x ... x (U_{k_d} - U_{k_d - 1}) f, +// with U_0 = 0. The levels grow by two points, not one: for even f, as x^2, the rules with 1 and 2 points can agree, +// and a zero contribution would stop the refinement in that coordinate. f takes a container of dimension Reals, resized if it can be. +template ::type> +auto sparse_grid(F f, const Rule& rule, unsigned dimension, std::size_t max_calls = 100000, + typename Rule::value_type tolerance = tools::root_epsilon(), + typename Rule::value_type* error = nullptr, typename Rule::value_type* pL1 = nullptr) + -> decltype(f(std::declval())) +{ + using Real = typename Rule::value_type; + using K = decltype(f(std::declval())); + using index = std::vector; + static_assert(!std::is_integral::value, + "The return type cannot be integral, it must be a real or complex floating point type, or a vector of them."); + constexpr const char* function = "boost::math::quadrature::sparse_grid<%1%>"; + Point x{}; + if constexpr (detail::leja_resizable::value) + { + if (x.size() == 0) + x.resize(dimension); + } + if (dimension == 0 || max_calls == 0 || rule.weights(1).empty() || static_cast(x.size()) != dimension) + { + Real nan = policies::raise_domain_error(function, "The dimension (%1%) and max_calls must be > 0, and the dimension the size of f's argument.", + Real(dimension), typename Rule::policy_type()); + if (error) + *error = nan; + if (pL1) + *pL1 = nan; + return detail::leja_nan(); + } + + // delta[k][j], j < 2k - 1, are the weights of U_k - U_{k-1}; the rule computes them as they are needed. + std::vector> delta(1); + auto possible = [&](unsigned level) { + return 2 * level - 1 <= rule.max_points(); + }; + auto grow = [&](unsigned level) { + for (unsigned k = static_cast(delta.size()); k <= level; ++k) + { + delta.push_back(rule.weights(2 * k - 1)); + for (unsigned j = 0; k > 1 && j < 2 * k - 3; ++j) + delta[k][j] -= rule.weights(2 * k - 3)[j]; + } + }; + grow(1); + + // The points, in a trie by their indices j_i into rule.abscissa(): kids[n][j] is the child of node n, at depth i, for + // j_i = j; at depth d - 1 it indexes points, which hold f and the weight of the point in the sparse grid. + constexpr std::size_t none = static_cast(-1); + std::vector> kids(1); + std::vector> points; + Real L1 = 0; + // Adds c times the product of delta[k_l][j_l], l >= i, times f(x_j) to sum over the box j_l < 2 k_l - 1. + auto walk = [&](auto& self, const index& k, std::size_t n, unsigned i, const Real& c, K& sum, bool& first) -> void { + const unsigned m = 2 * k[i] - 1; + if (kids[n].size() < m) + kids[n].resize(m, none); + for (unsigned j = 0; j < m; ++j) + { + const Real cj = c * delta[k[i]][j]; + x[i] = rule.abscissa()[j]; + std::size_t next = kids[n][j]; + if (next == none) + { + next = i + 1 < dimension ? kids.size() : points.size(); + kids[n][j] = next; + if (i + 1 < dimension) + kids.emplace_back(); + else + points.emplace_back(f(x), Real(0)); + } + if (i + 1 < dimension) + { + self(self, k, next, i + 1, cj, sum, first); + continue; + } + detail::leja_accumulate(sum, cj, points[next].first, first); + first = false; + points[next].second += cj; + L1 += cj * detail::leja_distance(points[next].first); + } + }; + // Adds Delta_k f to result, evaluating f at the new points, and returns its max norm. The sum of the Delta_k rounds far + // less than the sum over the points, whose weights are large and of both signs. + K result{}; + bool result_first = true; + auto add = [&](const index& k) { + K sum{}; + bool first = true; + walk(walk, k, 0, 0, Real(1), sum, first); + detail::leja_accumulate(result, Real(1), sum, result_first); + result_first = false; + return detail::leja_distance(sum); + }; + // The new points of k are those with j_i in {2 k_i - 3, 2 k_i - 2} where k_i > 1, and j_i = 0 elsewhere. + auto cost = [&](const index& k) { + std::size_t n = 1; + for (unsigned i = 0; i < dimension; ++i) + n *= k[i] > 1 ? 2 : 1; + return n; + }; + + // The active indices, with the norms of their contributions, and the old ones, whose forward neighbours are added. + std::map active; + std::set old; + // The active indices by decreasing norm, ties to the smallest, and the running sum of their norms; added, the sum of all + // the norms ever added, bounds its rounding error. + auto by_norm = [](const std::pair& a, const std::pair& b) { + return a.first > b.first || (a.first == b.first && a.second < b.second); + }; + std::set, decltype(by_norm)> queue(by_norm); + Real active_sum = 0, added = 0; + std::size_t updates = 0; + auto activate = [&](const index& m) { + const Real v = add(m); + active[m] = v; + if (v == v) + queue.emplace(v, m); + active_sum += v; + added += v; + ++updates; + }; + auto deactivate = [&](const index& k) { + auto a = active.find(k); + if (a->second == a->second) + queue.erase({a->second, k}); + active_sum -= a->second; + active.erase(a); + ++updates; + }; + // The contributions of old indices that the rule has too few points to refine stay in the error estimate. + Real unrefined = 0; + activate(index(dimension, 1)); + std::size_t calls = 1; + Real err = 0; + // When the estimate passes, every active index is refined once more, and the result is accepted only if the estimate + // still passes, since a contribution can vanish by symmetry or by chance. Nothing is accepted while L1 is 0, as it is + // while f has been 0 at every point: every active index is refined, until an interaction that vanishes there shows, or + // for dimension + 1 rounds, which reach the product of the levels 2 that all of them need, and then f is taken to be 0. + bool checked = false; + unsigned zero_rounds = 0; + for (;;) + { + err = unrefined + active_sum; + // Near the tolerance the running sum's rounding matters: sum again. + if (err <= tolerance * L1 + (updates + 1) * tools::epsilon() * (added + unrefined)) + { + active_sum = 0; + for (const auto& a : active) + active_sum += a.second; + added = active_sum; + updates = 0; + err = unrefined + active_sum; + } + std::vector batch; + // A NaN from f: no refinement can help. + if (active.empty() || (err <= tolerance * L1 && checked && (L1 > 0 || zero_rounds > dimension)) || err != err || L1 != L1) + break; + if (err <= tolerance * L1) + { + for (const auto& a : active) + batch.push_back(a.first); + checked = true; + zero_rounds += L1 == 0; + } + else + { + batch.push_back(queue.begin()->second); + checked = false; + } + // The admissible forward neighbours, whose backward neighbours are all old once the batch is. + std::set next; + std::size_t new_calls = 0; + Real capped = 0; + for (const index& k : batch) + old.insert(k); + for (const index& k : batch) + { + bool full = false; + for (unsigned i = 0; i < dimension; ++i) + { + index m = k; + if (!possible(++m[i])) + { + full = true; + continue; + } + bool admissible = !active.count(m) && !old.count(m); + for (unsigned l = 0; l < dimension && admissible; ++l) + { + if (m[l] > 1) + { + --m[l]; + admissible = old.count(m) > 0; + ++m[l]; + } + } + if (admissible && next.insert(m).second) + new_calls += cost(m); + } + if (full) + capped += active[k]; + } + if (calls + new_calls > max_calls) + { + // Out of calls before the check is done: the estimate is not to be trusted. + if (checked) + err = std::numeric_limits::quiet_NaN(); + break; + } + unrefined += capped; + for (const index& m : next) + grow(*std::max_element(m.begin(), m.end())); + for (const index& k : batch) + deactivate(k); + for (const index& m : next) + activate(m); + calls += new_calls; + } + + L1 = 0; + for (const auto& p : points) + L1 += p.second * detail::leja_distance(p.first); + if (error) + *error = err; + if (pL1) + *pL1 = L1; + return result; +} + +// The same over the box [a_0, b_0] x ... x [a_{d-1}, b_{d-1}], d = a.size(), for leja_legendre and leja_jacobi: each coordinate +// is mapped as integrate(f, a_i, b_i) maps it, and the result, error and L1 scale by the product of the factors. a and b +// are containers with operator[] and size(), such as std::vector, std::array or Eigen::VectorXd. +template ::type, + typename std::enable_if::value, bool>::type = true> +auto sparse_grid(F f, const Rule& rule, const Box& a, const Box& b, std::size_t max_calls = 100000, + typename Rule::value_type tolerance = tools::root_epsilon(), + typename Rule::value_type* error = nullptr, typename Rule::value_type* pL1 = nullptr) + -> decltype(f(std::declval())) +{ + using Real = typename Rule::value_type; + using K = decltype(f(std::declval())); + using std::abs; + static_assert(detail::leja_bounded::value, "Bounds need a rule on [-1, 1], leja_legendre or leja_jacobi."); + const std::size_t d = a.size(); + std::vector c(d), h(d); + Real J = 1; + const bool same_size = static_cast(b.size()) == d; + bool valid = same_size; + std::size_t i = 0; + for (; i < d && valid; ++i) + { + valid = (BOOST_MATH_NAMESPACE::isfinite)(Real(a[i])) && (BOOST_MATH_NAMESPACE::isfinite)(Real(b[i])); + c[i] = (Real(a[i]) + Real(b[i])) / 2; + h[i] = (Real(b[i]) - Real(a[i])) / 2; + J *= h[i] == 0 ? Real(0) : rule.jacobian(h[i]); + } + Point y{}; + if constexpr (detail::leja_resizable::value) + { + if (y.size() == 0) + y.resize(d); + } + if (!valid || (J == 0 && static_cast(y.size()) == d)) + { + Real r = valid ? Real(0) : policies::raise_domain_error("boost::math::quadrature::sparse_grid<%1%>", + same_size ? "The bounds must be finite, but got %1%." : "The bounds must be of the same size, but b has size %1%.", + same_size ? Real((BOOST_MATH_NAMESPACE::isfinite)(Real(a[i - 1])) ? b[i - 1] : a[i - 1]) : Real(b.size()), typename Rule::policy_type()); + if (error) + *error = r; + if (pL1) + *pL1 = r; + // A box of zero volume integrates to 0, of the shape of f(a). + K result = detail::leja_nan(); + for (std::size_t j = 0; valid && j < d; ++j) + y[j] = a[j]; + if (valid) + detail::leja_accumulate(result, Real(0), f(y), true); + return result; + } + auto g = [&](const Point& x) { + for (std::size_t i = 0; i < d; ++i) + y[i] = c[i] + h[i] * x[i]; + return f(y); + }; + K result; + detail::leja_accumulate(result, J, sparse_grid(g, rule, static_cast(d), max_calls, tolerance, error, pL1), true); + if (error) + *error *= abs(J); + if (pL1) + *pL1 *= abs(J); + return result; +} + +} // namespace quadrature +BOOST_MATH_NAMESPACE_END + +#endif // BOOST_MATH_QUADRATURE_SPARSE_GRID_HPP diff --git a/nanobind/src/quadrature.cpp b/nanobind/src/quadrature.cpp index 958c35dfce..602b7914c6 100644 --- a/nanobind/src/quadrature.cpp +++ b/nanobind/src/quadrature.cpp @@ -15,8 +15,10 @@ #include #include #include +#include #include #include +#include #include #include #include @@ -25,6 +27,7 @@ #include #include #include +#include #include #include @@ -107,6 +110,21 @@ struct py_gauss_jacobi : bq::gauss_jacobi py_gauss_jacobi(unsigned n_, double alpha_, double beta_) : bq::gauss_jacobi(n_, alpha_, beta_), n(n_), alpha(alpha_), beta(beta_) {} }; +// The Leja rules with parameters keep them to themselves; remember them for repr. +struct py_leja_laguerre : bq::leja_laguerre +{ + double alpha; + explicit py_leja_laguerre(double alpha_) : bq::leja_laguerre(alpha_), alpha(alpha_) {} +}; + +struct py_leja_jacobi : bq::leja_jacobi +{ + double alpha, beta; + py_leja_jacobi(double alpha_, double beta_) : bq::leja_jacobi(alpha_, beta_), alpha(alpha_), beta(beta_) {} +}; + +std::string py_repr(double x) { return nb::cast(nb::repr(nb::float_(x))); } + // The n-point Gauss rule of a family, computed afresh: (nodes ascending, weights). A symmetric family keeps only // its nonnegative nodes, so mirror them. template @@ -134,6 +152,59 @@ std::pair gauss_rule(unsigned n, Family const& family, bool sy return {to_array(x), to_array(w)}; } +// A nested Leja rule, whose integrate adds points until the tolerance; over [a, b] for the rules on [-1, 1]. +template +nb::class_ bind_leja(nb::module_& m, char const* name, char const* doc) +{ + nb::class_ c(m, name, doc); + if constexpr (bq::detail::leja_bounded::value) + c.def("integrate", [](Rule const& rule, nb::callable f, double a, double b, double tolerance) { + double error = 0, L1 = 0; + double Q = rule.integrate(scalar_function{f}, a, b, tolerance, &error, &L1); + return std::make_tuple(Q, error, L1); + }, nb::arg("f"), nb::arg("a"), nb::arg("b"), nb::kw_only(), nb::arg("tolerance") = root_eps, "(integral over [a, b], error estimate, L1 norm)"); + return c + .def("integrate", [](Rule const& rule, nb::callable f, double tolerance) { + double error = 0, L1 = 0; + double Q = rule.integrate(scalar_function{f}, tolerance, &error, &L1); + return std::make_tuple(Q, error, L1); + }, nb::arg("f"), nb::kw_only(), nb::arg("tolerance") = root_eps, "(integral, error estimate, L1 norm), adding one point at a time") + .def_prop_ro("max_points", [](Rule const&) { return Rule::max_points(); }) + .def_prop_ro("abscissa", [](Rule const& rule) { return to_array(rule.abscissa()); }, nb::rv_policy::move, "the points computed so far, in the order chosen") + .def("weights", [](Rule const& rule, unsigned n) { return to_array(rule.weights(n)); }, nb::kw_only(), nb::arg("n"), nb::rv_policy::move, + "the weights of the n-point rule, whose points are abscissa[:n]"); +} + +// sparse_grid(f, rule=..., ...), one overload per rule; over a box, the rule defaults to leja_legendre. +template +void bind_sparse_grid(nb::module_& m, nb::object default_rule = nb::object()) +{ + m.def("sparse_grid", [](nb::callable f, Rule const& rule, unsigned dimension, std::size_t max_calls, double tolerance) { + double error = 0, L1 = 0; + double Q = bq::sparse_grid([&f](std::vector const& x) { return nb::cast(f(to_array(x))); }, + rule, dimension, max_calls, tolerance, &error, &L1); + return std::make_tuple(Q, error, L1); + }, nb::arg("f"), nb::kw_only(), nb::arg("rule"), nb::arg("dimension"), nb::arg("max_calls") = 100000, nb::arg("tolerance") = root_eps, + "(integral, error estimate, L1) by the dimension-adaptive sparse grid of the rule; f takes a 1D array"); + if constexpr (bq::detail::leja_bounded::value) + { + auto box = [](nb::callable f, std::vector const& a, std::vector const& b, Rule const& rule, + std::size_t max_calls, double tolerance) { + double error = 0, L1 = 0; + double Q = bq::sparse_grid([&f](std::vector const& x) { return nb::cast(f(to_array(x))); }, + rule, a, b, max_calls, tolerance, &error, &L1); + return std::make_tuple(Q, error, L1); + }; + char const* doc = "(integral over the box [a, b], error estimate, L1); the dimension is len(a)"; + if (default_rule) + m.def("sparse_grid", box, nb::arg("f"), nb::arg("a"), nb::arg("b"), nb::kw_only(), nb::arg("rule") = default_rule, + nb::arg("max_calls") = 100000, nb::arg("tolerance") = root_eps, doc); + else + m.def("sparse_grid", box, nb::arg("f"), nb::arg("a"), nb::arg("b"), nb::kw_only(), nb::arg("rule"), + nb::arg("max_calls") = 100000, nb::arg("tolerance") = root_eps, doc); + } +} + void check_exponent(char const* name, double a) { if (!(a > -1) || !std::isfinite(a)) @@ -277,6 +348,30 @@ void bind_quadrature(nb::module_& m) + ", beta=" + nb::cast(nb::repr(nb::float_(rule.beta))) + ")"; }); + bind_leja>(m, "leja_legendre", "quadrature::leja_legendre(): nested Leja rules on [-1, 1]") + .def(nb::init<>()) + .def("__repr__", [](bq::leja_legendre const&) { return std::string("leja_legendre()"); }); + bind_leja>(m, "leja_hermite", "quadrature::leja_hermite(): nested Leja rules for exp(-x^2) on the real line") + .def(nb::init<>()) + .def("__repr__", [](bq::leja_hermite const&) { return std::string("leja_hermite()"); }); + bind_leja(m, "leja_laguerre", + "quadrature::leja_laguerre(alpha): nested Leja rules for x^alpha exp(-x) on [0, inf)") + .def(nb::init(), nb::kw_only(), nb::arg("alpha") = 0.0) + .def_ro("alpha", &py_leja_laguerre::alpha) + .def("__repr__", [](py_leja_laguerre const& rule) { return "leja_laguerre(alpha=" + py_repr(rule.alpha) + ")"; }); + bind_leja(m, "leja_jacobi", + "quadrature::leja_jacobi(alpha, beta): nested Leja rules for (1 - x)^alpha (1 + x)^beta on [-1, 1]") + .def(nb::init(), nb::kw_only(), nb::arg("alpha"), nb::arg("beta")) + .def_ro("alpha", &py_leja_jacobi::alpha) + .def_ro("beta", &py_leja_jacobi::beta) + .def("__repr__", [](py_leja_jacobi const& rule) { + return "leja_jacobi(alpha=" + py_repr(rule.alpha) + ", beta=" + py_repr(rule.beta) + ")"; + }); + bind_sparse_grid>(m, nb::cast(bq::leja_legendre())); + bind_sparse_grid>(m); + bind_sparse_grid(m); + bind_sparse_grid(m); + m.def("gauss_hermite_rule", [](unsigned n) { return gauss_rule(n, BOOST_MATH_NAMESPACE::detail::hermite_family(n), true); }, arg("n"), nb::rv_policy::move, "(nodes, weights) of the n-point rule for exp(-x^2), computed on each call; as numpy's hermgauss"); diff --git a/nanobind/tests/test_smoke.py b/nanobind/tests/test_smoke.py index b24298ab64..2016604dbc 100644 --- a/nanobind/tests/test_smoke.py +++ b/nanobind/tests/test_smoke.py @@ -278,7 +278,7 @@ def test_array_matches_scalar(case): SIEVE = {"prime_sieve", "prime_range", "prime_count", "polynomial_roots", "polynomial_roots_aberth", "polynomial_roots_companion_matrix"} QUADRATURE = {"tanh_sinh", "exp_sinh", "sinh_sinh", "trapezoidal", "gauss", "gauss_kronrod", "gauss_hermite", "gauss_laguerre", "gauss_hermite_averaged", "gauss_laguerre_averaged", "gauss_hermite_rule", "gauss_laguerre_rule", "gauss_jacobi_rule", - "ooura_fourier_sin", "ooura_fourier_cos", "van_den_bos_unit_square"} + "ooura_fourier_sin", "ooura_fourier_cos", "van_den_bos_unit_square", "sparse_grid"} OPTIMIZATION = {"differential_evolution", "jso", "random_search", "multilevel_coordinate_search", "cma_es"} @@ -325,6 +325,17 @@ def second_moment(rule): (second_moment(bm.gauss_hermite_rule(5)), math.sqrt(math.pi) / 2), (second_moment(bm.gauss_laguerre_rule(5, alpha=0.5)), math.gamma(3.5)), (second_moment(bm.gauss_jacobi_rule(5, 0.5, -0.5)), math.pi / 2), + (bm.leja_legendre().integrate(lambda x: x * x), 2 / 3), + (bm.leja_hermite().integrate(lambda x: x * x), math.sqrt(math.pi) / 2), + (bm.leja_legendre().integrate(math.exp), math.e - 1 / math.e), + (bm.leja_hermite().integrate(math.cos), math.sqrt(math.pi) * math.exp(-0.25)), + (bm.leja_laguerre().integrate(lambda x: x * x), 2), + (bm.leja_laguerre(alpha=0.5).integrate(lambda x: 1.0), math.sqrt(math.pi) / 2), + (bm.leja_jacobi(alpha=-0.5, beta=-0.5).integrate(lambda x: x * x), math.pi / 2), + (bm.leja_legendre().integrate(math.exp, -2, 5), math.exp(5) - math.exp(-2)), + (bm.leja_legendre().integrate(lambda x: x, 2, 0, tolerance=1e-12), -2), + # (4 - x)^alpha (x - 1)^beta over [1, 4]: 3^(alpha + beta + 1) B(alpha + 1, beta + 1). + (bm.leja_jacobi(alpha=0.5, beta=-0.25).integrate(lambda x: 1.0, 1, 4), 3 ** 1.25 * math.gamma(1.5) * math.gamma(0.75) / math.gamma(2.25)), (bm.ooura_fourier_sin(lambda t: 1 / t, 1.0), math.pi / 2), (bm.ooura_fourier_cos(lambda t: 1 / (1 + t * t), 1.0), math.pi / (2 * math.e)), (bm.van_den_bos_unit_square(lambda x, y: x * y), 0.25), @@ -339,6 +350,43 @@ def second_moment(rule): assert abs(Q - math.sqrt(math.pi) * math.exp(-0.25)) <= error < 1e-6 with pytest.raises(RuntimeError): bm.gauss_jacobi(10, -0.9, -0.9).integrate_averaged(math.cos) + rule = bm.leja_legendre() + assert len(rule.weights(n=5)) == 5 and len(rule.abscissa) == 5 and rule.abscissa[0] == 0.0 + with pytest.raises(ValueError): + rule.weights(n=rule.max_points + 1) + with pytest.raises(TypeError): + rule.integrate(math.cos, 5) + assert rule.integrate(math.cos, 1, 1) == (0, 0, 0) + with pytest.raises(ValueError): + rule.integrate(math.cos, 0, math.inf) + with pytest.raises(TypeError): + bm.leja_hermite().integrate(math.cos, 0, 1) + # The box [0, 1] x [-1, 2] x [3, 5], from lists or arrays. + box = (math.exp(0.5) - 1) * 2 * (math.exp(1) - math.exp(-2)) * 4 * (math.exp(1.25) - math.exp(0.75)) + for a, b in [([0, -1, 3], [1, 2, 5]), (np.array([0.0, -1, 3]), np.array([1.0, 2, 5]))]: + Q, error, L1 = bm.sparse_grid(lambda x: math.exp(x[0] / 2 - x[1] + x[2] / 4), a, b, rule=rule, max_calls=100000) + assert abs(Q - box) <= max(error, 1e-13 * L1) <= 1.5e-8 * L1 + with pytest.raises(ValueError): + bm.sparse_grid(lambda x: 1.0, [0, 0], [1], rule=rule, max_calls=100) + with pytest.raises(TypeError): + bm.sparse_grid(lambda x: 1.0, [0], [1], rule=bm.leja_hermite(), max_calls=100) + Q, error, L1 = bm.sparse_grid(lambda x: math.exp(x[0] + x[1] / 2), rule=bm.leja_legendre(), dimension=2, max_calls=1000) + assert abs(Q - (math.e - 1 / math.e) * 2 * (math.exp(0.5) - math.exp(-0.5))) <= max(error, 1e-14) <= 1e-7 * L1 + Q, error, L1 = bm.sparse_grid(lambda x: x[2] ** 2, rule=bm.leja_hermite(), dimension=3, max_calls=1000) + assert abs(Q - math.pi ** 1.5 / 2) < 1e-13 + Q, error, L1 = bm.sparse_grid(lambda x: x[0] * x[1], rule=bm.leja_laguerre(), dimension=2, max_calls=1000) + assert abs(Q - 1) < 1e-13 + # 0 at the first point and on the axes through it. + Q, error, L1 = bm.sparse_grid(lambda x: x[0] ** 2 * x[1] ** 2, [-1, -1], [1, 1], rule=rule) + assert abs(Q - 4 / 9) < 1e-13 + Q, error, L1 = bm.sparse_grid(lambda x: x[0] ** 2 * x[1] ** 2, rule=bm.leja_hermite(), dimension=2) + assert abs(Q - math.pi / 4) < 1e-13 + with pytest.raises(ValueError): + bm.leja_jacobi(alpha=-1, beta=0) + with pytest.raises(ValueError): + bm.leja_laguerre(alpha=math.inf) + with pytest.raises(ValueError): + bm.leja_jacobi(alpha=math.inf, beta=0) with pytest.raises(ZeroDivisionError): bm.tanh_sinh(lambda x: 1 / 0, 0, 1) @@ -486,7 +534,7 @@ def test_array_errors_warn_and_give_nan(): def test_every_distribution_has_a_case(): names = {n for n in bm.__all__ if inspect.isclass(getattr(bm, n))} - assert names - TRANSFORMS - {"gauss_jacobi"} == {c[0] for c in CONTINUOUS + DISCRETE} + assert names - TRANSFORMS - {"gauss_jacobi", "leja_legendre", "leja_hermite", "leja_laguerre", "leja_jacobi"} == {c[0] for c in CONTINUOUS + DISCRETE} assert TRANSFORMS <= names diff --git a/reporting/performance/leja_quadrature_performance.cpp b/reporting/performance/leja_quadrature_performance.cpp new file mode 100644 index 0000000000..cb13345514 --- /dev/null +++ b/reporting/performance/leja_quadrature_performance.cpp @@ -0,0 +1,136 @@ +// Copyright Nick Thompson 2026. +// Use, modification and distribution are subject to the +// Boost Software License, Version 1.0. (See accompanying file +// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) + +/* + * Google Benchmark suite for the Leja rules: building the points and weights, the adaptive rule, + * and the sparse grid; the calls counter is the number of evaluations of f. + * + * Example: + * c++ -O3 -DNDEBUG -std=c++17 -Iinclude \ + * reporting/performance/leja_quadrature_performance.cpp \ + * -lbenchmark -lpthread -o leja_bench + */ + +#include +#include +#include + +#include +#include +#include + +namespace { + +using boost::math::quadrature::leja_hermite; +using boost::math::quadrature::leja_jacobi; +using boost::math::quadrature::leja_laguerre; +using boost::math::quadrature::leja_legendre; + +template +void BM_construct(benchmark::State& state) +{ + for (auto _ : state) + { + Rule rule; + benchmark::DoNotOptimize(rule.weights(static_cast(state.range(0))).data()); + } +} +BENCHMARK_TEMPLATE(BM_construct, leja_legendre)->RangeMultiplier(2)->Range(16, 128)->Unit(benchmark::kMicrosecond); +BENCHMARK_TEMPLATE(BM_construct, leja_hermite)->RangeMultiplier(2)->Range(16, 128)->Unit(benchmark::kMicrosecond); +BENCHMARK_TEMPLATE(BM_construct, leja_laguerre)->RangeMultiplier(2)->Range(16, 128)->Unit(benchmark::kMicrosecond); +void BM_construct_jacobi(benchmark::State& state) +{ + for (auto _ : state) + { + leja_jacobi rule(-0.5, 2.0); + benchmark::DoNotOptimize(rule.weights(static_cast(state.range(0))).data()); + } +} +BENCHMARK(BM_construct_jacobi)->RangeMultiplier(2)->Range(16, 128)->Unit(benchmark::kMicrosecond); +BENCHMARK_TEMPLATE(BM_construct, leja_legendre)->Arg(64)->Unit(benchmark::kMicrosecond); +BENCHMARK_TEMPLATE(BM_construct, leja_legendre)->Arg(64)->Unit(benchmark::kMicrosecond); + +// A rule with the points the benchmarks need. +template +const Rule& prepared() +{ + static const Rule r; + r.weights(128); + return r; +} +const auto& legendre = prepared>; +const auto& hermite = prepared>; + +// s keeps the integrand from being hoisted out of the loop. +template +void integrate(benchmark::State& state, const Rule& rule, F f, double tolerance) +{ + double s = 0; + unsigned calls = 0; + for (auto _ : state) + { + calls = 0; + benchmark::DoNotOptimize(s); + auto Q = rule.integrate([&](double x) { ++calls; return f(x + s); }, tolerance); + benchmark::DoNotOptimize(Q); + } + state.counters["calls"] = calls; +} + +auto cos3 = [](double x) { return std::cos(3 * x); }; +auto runge = [](double x) { return 1 / (1 + 16 * x * x); }; +auto lorentz = [](double x) { return 1 / (1 + x * x); }; +auto three = [](double x) { return std::array{std::exp(x), std::cos(3 * x), 1 / (1 + 16 * x * x)}; }; + +void BM_adaptive_legendre_cos3x(benchmark::State& state) { integrate(state, legendre(), cos3, 1e-12); } +BENCHMARK(BM_adaptive_legendre_cos3x); +void BM_adaptive_legendre_runge(benchmark::State& state) { integrate(state, legendre(), runge, 1e-6); } +BENCHMARK(BM_adaptive_legendre_runge); +void BM_adaptive_hermite_cos3x(benchmark::State& state) { integrate(state, hermite(), cos3, 1e-12); } +BENCHMARK(BM_adaptive_hermite_cos3x); +void BM_adaptive_hermite_lorentz(benchmark::State& state) { integrate(state, hermite(), lorentz, 1e-6); } +BENCHMARK(BM_adaptive_hermite_lorentz); +void BM_adaptive_legendre_vector3(benchmark::State& state) { integrate(state, legendre(), three, 1e-6); } +BENCHMARK(BM_adaptive_legendre_vector3); + +// exp(sum_i x_i / 2^i) and cos(sum_i x_i / (i + 1)) in range(0) dimensions. +template +void sparse(benchmark::State& state, const Rule& rule, F f, double tolerance) +{ + const unsigned d = static_cast(state.range(0)); + unsigned calls = 0; + for (auto _ : state) + { + calls = 0; + auto Q = boost::math::quadrature::sparse_grid([&](const std::vector& x) { ++calls; return f(x); }, rule, d, 1000000, tolerance); + benchmark::DoNotOptimize(Q); + } + state.counters["calls"] = calls; +} + +void BM_sparse_legendre_exp(benchmark::State& state) +{ + sparse(state, legendre(), [](const std::vector& x) { + double s = 0; + for (std::size_t i = 0; i < x.size(); ++i) + s += x[i] / (1 << i); + return std::exp(s); + }, 1e-8); +} +BENCHMARK(BM_sparse_legendre_exp)->Arg(4)->Arg(8)->Unit(benchmark::kMillisecond); +void BM_sparse_hermite_cos(benchmark::State& state) +{ + sparse(state, hermite(), [](const std::vector& x) { + double s = 0; + for (std::size_t i = 0; i < x.size(); ++i) + s += x[i] / (i + 1); + return std::cos(s); + }, 1e-6); +} +BENCHMARK(BM_sparse_hermite_cos)->Arg(4)->Arg(8)->Unit(benchmark::kMillisecond); + +} // namespace + +BENCHMARK_MAIN(); diff --git a/test/Jamfile.v2 b/test/Jamfile.v2 index 789702c032..06ee533a3c 100644 --- a/test/Jamfile.v2 +++ b/test/Jamfile.v2 @@ -1295,6 +1295,7 @@ test-suite quadrature : [ run test_gauss_jacobi_quadrature.cpp : : : $(float128_type) [ requires cxx17_if_constexpr ] off msvc:/bigobj release ] [ run test_gauss_jacobi_quadrature_eigen.cpp : : : [ requires cxx17_if_constexpr ] [ check-target-builds ../../multiprecision/config//has_eigen : : no ] ] + [ run test_leja_quadrature.cpp : : : [ requires cxx17_if_constexpr cxx11_hdr_thread ] msvc:/bigobj gcc-mingw:-Wa,-mbig-obj linux:"-pthread" ] [ run test_naive_monte_carlo.cpp : : : msvc:/bigobj TEST=1 [ requires cxx11_auto_declarations cxx11_lambdas cxx11_unified_initialization_syntax cxx11_hdr_thread cxx11_hdr_atomic cxx11_decltype cxx11_hdr_future cxx11_hdr_chrono cxx11_hdr_random cxx11_allocator ] diff --git a/test/test_leja_quadrature.cpp b/test/test_leja_quadrature.cpp new file mode 100644 index 0000000000..d262d4ab5d --- /dev/null +++ b/test/test_leja_quadrature.cpp @@ -0,0 +1,708 @@ +// Copyright Nick Thompson, 2026 +// Use, modification and distribution are subject to the +// Boost Software License, Version 1.0. +// (See accompanying file LICENSE_1_0.txt +// or copy at http://www.boost.org/LICENSE_1_0.txt) + +#include "math_unit_test.hpp" +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#ifdef BOOST_HAS_FLOAT128 +#include +#endif + +using boost::math::quadrature::leja_hermite; +using boost::math::quadrature::leja_jacobi; +using boost::math::quadrature::leja_laguerre; +using boost::math::quadrature::leja_legendre; + +// The n-point rule, sum_i w_i f(x_i), and L1 = sum_i |w_i| |f(x_i)|, for a real or complex f. +template +auto n_point(const Rule& rule, F f, unsigned n, typename Rule::value_type* L1 = nullptr) +{ + using std::abs; + const std::vector& w = rule.weights(n); + decltype(f(w[0])) Q(0); + typename Rule::value_type l = 0; + for (unsigned i = 0; i < n; ++i) + { + const auto y = f(rule.abscissa()[i]); + Q += w[i] * y; + l += abs(w[i]) * abs(y); + } + if (L1) + *L1 = l; + return Q; +} + +// The integral of x^k over [-1, 1], or of x^k exp(-x^2) over the real line, Gamma((k + 1)/2) for even k. +template +Real moment(const leja_legendre&, unsigned k) +{ + return k % 2 ? Real(0) : Real(2) / (k + 1); +} + +template +Real moment(const leja_hermite&, unsigned k) +{ + Real m = boost::math::constants::root_pi(); + for (unsigned j = 2; j <= k; j += 2) + m *= Real(j - 1) / 2; + return k % 2 ? Real(0) : m; +} + + +// Rounding relative to L1, in units of eps: the small Laguerre weights far out, where x^k is large, are the least accurate. +template +unsigned slack(const Rule&, unsigned n) +{ + return 8 * n; +} + +template +unsigned slack(const leja_laguerre&, unsigned n) +{ + return 32 * n; +} + +// The n-point rule is exact for polynomials of degree n - 1, up to rounding relative to L1 = sum |w_i x_i^k|; moment(k) is +// the integral of x^k. The sum of |w_i| is at most mass times the integral of 1. +template +void test_exactness(const Rule& rule, Moment moment, unsigned N, double mass = 1.25) +{ + using Real = typename Rule::value_type; + using std::abs; + using std::pow; + const Real eps = std::numeric_limits::epsilon(); + for (unsigned n = 1; n <= N; ++n) + { + CHECK_EQUAL(rule.weights(n).size(), std::size_t(n)); + Real sum = 0; + for (const Real& w : rule.weights(n)) + sum += abs(w); + // The weights can be negative, but not by much. + CHECK_LE(sum, Real(mass) * moment(0)); + for (unsigned k = 0; k < n; ++k) + { + Real L1; + const Real Q = n_point(rule, [k](const Real& x) { return pow(x, Real(k)); }, n, &L1); + CHECK_LE(abs(Q - moment(k)), slack(rule, n) * eps * L1); + } + } +} + +// A rule computed to 40 points at once is the rule computed one point at a time, and its copy. +template +void check_nested(const Rule& rule) +{ + const Rule steps(rule); + rule.weights(40); + const Rule copy(rule); + CHECK_EQUAL(copy.abscissa().size(), std::size_t(40)); + // Assigned, a rule keeps the points and extends them as the original does. + Rule assigned(steps); + assigned = copy; + CHECK_EQUAL(assigned.weights(41)[40], rule.weights(41)[40]); + for (unsigned n = 1; n <= 40; ++n) + { + const std::vector& w = steps.weights(n); + CHECK_EQUAL(steps.abscissa()[n - 1], rule.abscissa()[n - 1]); + for (unsigned i = 0; i < n; ++i) + { + CHECK_EQUAL(w[i], rule.weights(n)[i]); + CHECK_EQUAL(copy.weights(n)[i], rule.weights(n)[i]); + } + } +} + +template +void test_exactness(unsigned n = std::is_same::value ? 16 : 32) +{ + using boost::math::tgamma; + const leja_legendre legendre; + const leja_hermite hermite; + test_exactness(legendre, [&](unsigned k) { return moment(legendre, k); }, n); + test_exactness(hermite, [&](unsigned k) { return moment(hermite, k); }, n); + for (Real a : {Real(0), Real(-0.5L), Real(2.5L)}) + test_exactness(leja_laguerre(a), [a](unsigned k) { return tgamma(k + a + 1); }, n); + // n-point Gauss-Jacobi integrates x^k, k < 2n, to rounding. + for (Real a : {Real(-0.5L), Real(0.5L)}) + { + for (Real b : {Real(-0.75L), Real(2)}) + { + boost::math::quadrature::gauss_jacobi gauss(n, a, b); + // A negative exponent places the points as if it were 0, and the weights are more negative. + test_exactness(leja_jacobi(a, b), [&](unsigned k) { return gauss.integrate([k](const Real& x) { using std::pow; return pow(x, Real(k)); }); }, n, 1.5); + } + } +} + +template +void test_nested() +{ + leja_legendre large; + leja_hermite h_large; + check_nested(large); + check_nested(h_large); + check_nested(leja_laguerre(Real(1.5))); + check_nested(leja_jacobi(Real(-0.5), Real(2))); + // The first points: the mean of the weight, then 0, -1, 1 for Legendre, where exact ties go to the left. + CHECK_EQUAL(large.abscissa()[0], Real(0)); + CHECK_ULP_CLOSE(-1 / boost::math::constants::root_three(), large.abscissa()[3], 2); + CHECK_EQUAL(h_large.abscissa()[0], Real(0)); + CHECK_ULP_CLOSE(Real(-1), h_large.abscissa()[1], 2); + leja_laguerre laguerre; + laguerre.weights(2); + CHECK_EQUAL(laguerre.abscissa()[1], Real(0)); + leja_jacobi chebyshev(Real(-0.5), Real(-0.5)); + chebyshev.weights(3); + CHECK_EQUAL(chebyshev.abscissa()[1], Real(-1)); + CHECK_EQUAL(chebyshev.abscissa()[2], Real(1)); +} + +// Every precision picks the same point at every step. +template class Rule, class... Args> +void test_same_points(unsigned points, const Args&... args) +{ + using std::abs; + Rule> rule{Real(args)...}; + Rule> reference{double(args)...}; + rule.weights(points); + reference.weights(points); + const Real tol = 8 * (std::max)(Real(std::numeric_limits::epsilon()), std::numeric_limits::epsilon()); + for (unsigned i = 0; i < points; ++i) + CHECK_ABSOLUTE_ERROR(Real(reference.abscissa()[i]), rule.abscissa()[i], tol * (std::max)(Real(1), abs(rule.abscissa()[i]))); +} + +template +void test_smooth(unsigned n) +{ + using std::cos; + using std::exp; + using std::sin; + const Real eps = std::numeric_limits::epsilon(); + auto f = [](const Real& x) { return cos(3 * x); }; + leja_legendre rule; + CHECK_ABSOLUTE_ERROR(2 * sin(Real(3)) / 3, n_point(rule, f, n), 16 * eps); + // Hermite needs twice as many points: Leja rules are exact to degree n - 1, not 2n - 1. + leja_hermite h_rule; + CHECK_ABSOLUTE_ERROR(boost::math::constants::root_pi() * exp(Real(-9) / 4), n_point(h_rule, f, 2 * n), 16 * eps); +} + +template +void test_adaptive() +{ + using std::abs; + using std::atan; + using std::cos; + using std::exp; + using std::sin; + using std::sqrt; + const Real eps = std::numeric_limits::epsilon(); + leja_legendre rule; + leja_hermite h_rule; + leja_laguerre l_rule; + leja_laguerre l_half(Real(-0.5L)); + leja_jacobi chebyshev(Real(-0.5L), Real(-0.5L)); + Real error, L1; + // On these integrands the error estimate bounds the error, or the error is at the level of rounding. + auto check = [&](const Real& expected, const Real& Q, const Real& tolerance) { + CHECK_LE(error, tolerance * L1); + CHECK_LE(abs(Q - expected), (std::max)(error, 64 * eps * L1)); + }; + Real Q = rule.integrate([](const Real& x) { return cos(3 * x); }, Real(1e-12L), &error, &L1); + check(2 * sin(Real(3)) / 3, Q, Real(1e-12L)); + Q = rule.integrate([](const Real& x) { return 1 / (1 + 16 * x * x); }, Real(1e-6L), &error, &L1); + check(atan(Real(4)) / 2, Q, Real(1e-6L)); + Q = rule.integrate([](const Real& x) { return exp(x); }, boost::math::tools::root_epsilon(), &error, &L1); + check(exp(Real(1)) - exp(Real(-1)), Q, boost::math::tools::root_epsilon()); + Q = h_rule.integrate([](const Real& x) { return cos(3 * x); }, Real(1e-12L), &error, &L1); + check(boost::math::constants::root_pi() * exp(Real(-9) / 4), Q, Real(1e-12L)); + Q = h_rule.integrate([](const Real& x) { return 1 / (1 + x * x / 4); }, Real(1e-8L), &error, &L1); + // 2 pi e^4 erfc(2) + check(Real(1.604698360911363269322784609663024807L), Q, Real(1e-8L)); + Q = l_rule.integrate([](const Real& x) { return cos(x); }, Real(1e-10L), &error, &L1); + check(Real(0.5L), Q, Real(1e-10L)); + Q = l_rule.integrate([](const Real& x) { return 1 / (1 + x); }, Real(1e-6L), &error, &L1); + // e E_1(1) + check(Real(0.5963473623231940743410784993692793760742L), Q, Real(1e-6L)); + // Gamma(1/2) cos(pi/8) / 2^(1/4) and pi J_0(1). + Q = l_half.integrate([](const Real& x) { return cos(x); }, Real(1e-10L), &error, &L1); + check(Real(1.376996331853153438664376624047155467474L), Q, Real(1e-10L)); + Q = chebyshev.integrate([](const Real& x) { return cos(x); }, Real(1e-12L), &error, &L1); + check(Real(2.403939430634412998273324892591511223785L), Q, Real(1e-12L)); +} + +// e^{ix} over [-1, 1] with n points, and with exp(-x^2) with h_n. +template > +void test_complex(unsigned n, unsigned h_n) +{ + using std::cos; + using std::exp; + using std::sin; + const Real eps = std::numeric_limits::epsilon(); + const Real tolerance = boost::math::tools::root_epsilon(); + auto f = [](const Real& x) { return Complex(cos(x), sin(x)); }; + leja_legendre rule; + leja_hermite h_rule; + Real error, L1; + Complex Q = n_point(rule, f, n); + CHECK_ABSOLUTE_ERROR(2 * sin(Real(1)), Real(Q.real()), 8 * eps); + CHECK_ABSOLUTE_ERROR(Real(0), Real(Q.imag()), 8 * eps); + Q = rule.integrate(f, tolerance, &error, &L1); + CHECK_ABSOLUTE_ERROR(2 * sin(Real(1)), Real(Q.real()), (std::max)(error, 8 * eps)); + CHECK_LE(error, tolerance * L1); + Q = n_point(h_rule, f, h_n); + const Real expected = boost::math::constants::root_pi() * exp(Real(-1) / 4); + CHECK_ABSOLUTE_ERROR(expected, Real(Q.real()), 8 * eps); + CHECK_ABSOLUTE_ERROR(Real(0), Real(Q.imag()), 8 * eps); + Q = h_rule.integrate(f, tolerance, &error, &L1); + CHECK_ABSOLUTE_ERROR(expected, Real(Q.real()), (std::max)(error, 8 * eps)); + CHECK_LE(error, tolerance * L1); +} + +// A vector-valued integrand is integrated componentwise, and the adaptive rule waits for the slowest component. +template > +void test_vector() +{ + using std::abs; + using std::cos; + using std::exp; + using std::sin; + const Real eps = std::numeric_limits::epsilon(); + const Real tolerance = Real(1e-6L); + leja_legendre rule; + unsigned calls = 0; + auto f = [&calls](const Real& x) { + ++calls; + return std::array{x * x, cos(3 * x), 1 / (1 + x * x)}; + }; + const std::array expected{Real(2) / 3, 2 * sin(Real(3)) / 3, boost::math::constants::half_pi()}; + Real error, L1; + // A constant converges at 4 points, and L1 is the max norm 2 times the sum of |w_i|. + std::array C = rule.integrate([](const Real&) { return std::array{Real(1), Real(-2)}; }, tolerance, &error, &L1); + CHECK_ULP_CLOSE(Real(-4), C[1], 4); + CHECK_ULP_CLOSE(2 * n_point(rule, [](const Real&) { return Real(1); }, 4), L1, 4); + std::array Q = rule.integrate(f, tolerance, &error, &L1); + CHECK_LE(error, tolerance * L1); + for (std::size_t j = 0; j < 3; ++j) + CHECK_LE(abs(Q[j] - expected[j]), (std::max)(error, 64 * eps * L1)); + // A NaN component is not hidden by the max norm. + rule.integrate([](const Real&) { return std::array{std::numeric_limits::quiet_NaN(), Real(1)}; }, tolerance, &error, &L1); + CHECK_NAN(L1); + unsigned slowest = 0; + rule.integrate([&slowest](const Real& x) { ++slowest; return 1 / (1 + x * x); }, tolerance); + CHECK_LE(slowest, calls); + + // A std::vector of complex values, with the Hermite weight. + leja_hermite h_rule; + auto g = [](const Real& x) { return std::vector{Complex(cos(x), sin(x)), Complex(x * x, Real(0))}; }; + std::vector P = h_rule.integrate(g, tolerance, &error, &L1); + CHECK_EQUAL(P.size(), std::size_t(2)); + CHECK_LE(error, tolerance * L1); + CHECK_ABSOLUTE_ERROR(boost::math::constants::root_pi() * exp(Real(-1) / 4), Real(P[0].real()), (std::max)(error, 8 * eps)); + CHECK_ABSOLUTE_ERROR(Real(0), Real(P[0].imag()), (std::max)(error, 8 * eps)); + CHECK_ABSOLUTE_ERROR(boost::math::constants::root_pi() / 2, Real(P[1].real()), 8 * eps); +} + +// The sparse grid, on integrands with known integrals. It stops when its error estimate is below tolerance * L1. +template > +void test_sparse_grid(const Real& tolerance) +{ + using boost::math::quadrature::sparse_grid; + using std::abs; + using std::cos; + using std::exp; + using std::pow; + using std::sin; + using point = std::vector; + const Real eps = std::numeric_limits::epsilon(); + const Real pi = boost::math::constants::pi(); + leja_legendre rule; + leja_hermite h_rule; + Real error, L1; + auto check = [&](const Real& expected, const Real& Q) { + CHECK_LE(error, tolerance * L1); + CHECK_LE(abs(Q - expected), (std::max)(error, 64 * eps * L1)); + }; + + // Constants, odd functions, and x^2, whose rules with 1 and 2 points agree. + check(Real(32), sparse_grid([](const point&) { return Real(1); }, rule, 5, 1000, tolerance, &error, &L1)); + check(pow(pi, Real(5) / 2), sparse_grid([](const point&) { return Real(1); }, h_rule, 5, 1000, tolerance, &error, &L1)); + check(Real(0), sparse_grid([](const point& x) { return x[0] * x[1] * x[1] * x[1]; }, rule, 3, 1000, tolerance, &error, &L1)); + check(Real(2) / 3, sparse_grid([](const point& x) { return x[0] * x[0]; }, rule, 1, 1000, tolerance, &error, &L1)); + check(pi / 2, sparse_grid([](const point& x) { return x[1] * x[1]; }, h_rule, 2, 1000, tolerance, &error, &L1)); + // With the Laguerre weight, exp(-x_0/2 - x_1/3 - x_2/4) integrates to 2/3 * 3/4 * 4/5, and x_0 x_1, which vanishes on + // the axes through 0, to 1. + check(Real(2) / 5, sparse_grid([](const point& x) { return exp(-x[0] / 2 - x[1] / 3 - x[2] / 4); }, leja_laguerre(), 3, + 100000, tolerance, &error, &L1)); + check(Real(1), sparse_grid([](const point& x) { return x[0] * x[1]; }, leja_laguerre(), 2, 1000, tolerance, &error, &L1)); + // f may take a std::vector, be generic, or take a std::array, as a lambda or a function pointer. + static const auto e = [](const auto& x) { return exp(x[0] / 2 - x[1] + x[2] / 4); }; + const Real Qv = sparse_grid([](const point& x) { return e(x); }, rule, 3, 10000, tolerance); + CHECK_EQUAL(Qv, sparse_grid(e, rule, 3, 10000, tolerance)); + CHECK_EQUAL(Qv, sparse_grid([](const std::array& x) { return e(x); }, rule, 3, 10000, tolerance)); + CHECK_EQUAL(Qv, sparse_grid(+[](std::array x) { return e(x); }, rule, 3, 10000, tolerance)); + // With the Chebyshev weight, cos(x_0 + x_1/2) integrates to pi^2 J_0(1) J_0(1/2), the same from a rule that grows + // as from one that has its points already. + auto c2 = [](const point& x) { return cos(x[0] + x[1] / 2); }; + const leja_jacobi chebyshev(Real(-0.5L), Real(-0.5L)); + const Real Qc = sparse_grid(c2, chebyshev, 2, 10000, tolerance, &error, &L1); + check(Real(7.087510228266830598984765096456477556883L), Qc); + chebyshev.weights(63); + CHECK_EQUAL(Qc, sparse_grid(c2, chebyshev, 2, 10000, tolerance, &error, &L1)); + // A NaN stops the grid after the refinement that found it. + unsigned nan_calls = 0; + CHECK_NAN(sparse_grid([&](const point& x) { ++nan_calls; return x[0] > 0 ? Real(1) : std::numeric_limits::quiet_NaN(); }, + rule, 3, 1000, tolerance, &error, &L1)); + CHECK_LE(nan_calls, 1u); + // Out of calls before the check is done: a NaN error estimate. + sparse_grid([](const point& x) { return x[0] * x[0]; }, rule, 1, 1, tolerance, &error, &L1); + CHECK_NAN(error); + // A product of polynomials of degree 4, exact once every level is 3. + auto p = [](const Real& t) { return 1 + t - 2 * t * t + t * t * t * t; }; + check(Real(16 * 16 * 16) / (15 * 15 * 15), + sparse_grid([&](const point& x) { return p(x[0]) * p(x[1]) * p(x[2]); }, rule, 3, 1000, tolerance, &error, &L1)); + + // Smooth and anisotropic: exp(sum_i x_i / 2^i) and, with the Hermite weight, cos(sum_i x_i / (i + 1)). + const unsigned d = 6; + Real expected = 1; + for (unsigned i = 0; i < d; ++i) + expected *= (exp(Real(1) / (1 << i)) - exp(-Real(1) / (1 << i))) * (1 << i); + std::set seen; + unsigned calls = 0; + auto f = [&](const point& x) { + ++calls; + CHECK_EQUAL(seen.insert(x).second, true); + Real s = 0; + for (unsigned i = 0; i < d; ++i) + s += x[i] / (1 << i); + return exp(s); + }; + check(expected, sparse_grid(f, rule, d, 100000, tolerance, &error, &L1)); + Real a2 = 0; + for (unsigned i = 0; i < 4; ++i) + a2 += 1 / Real((i + 1) * (i + 1)); + auto g = [](const point& x) { + Real s = 0; + for (unsigned i = 0; i < 4; ++i) + s += x[i] / (i + 1); + return Complex(cos(s), sin(s)); + }; + Complex Q = sparse_grid(g, h_rule, 4, 100000, tolerance, &error, &L1); + check(pi * pi * exp(-a2 / 4), Real(Q.real())); + CHECK_LE(Real(abs(Q.imag())), (std::max)(error, 64 * eps * L1)); + // A vector is integrated componentwise, with the max norm. + auto h = [&](const point& x) { return std::array{g(x), Complex(x[0] * x[0], Real(0))}; }; + std::array R = sparse_grid(h, h_rule, 4, 100000, tolerance, &error, &L1); + check(pi * pi * exp(-a2 / 4), Real(R[0].real())); + check(pi * pi / 2, Real(R[1].real())); + + // Out of calls: never more than max_calls, and an error estimate above the tolerance. + calls = 0; + seen.clear(); + sparse_grid(f, rule, d, 50, tolerance, &error, &L1); + CHECK_LE(calls, 50u); + CHECK_LE(tolerance * L1, error); +} + +// Functions that vanish at the first point, the mean of the weight, and on the axes through it. +void test_vanishing() +{ + using boost::math::quadrature::sparse_grid; + using point = std::vector; + const double eps = std::numeric_limits::epsilon(); + const double tolerance = boost::math::tools::root_epsilon(); + const double e = std::exp(1.0); + double error, L1; + auto check = [&](const double& expected, const double& Q) { + CHECK_LE(error, tolerance * L1); + CHECK_LE(std::abs(Q - expected), (std::max)(error, 64 * eps * L1)); + }; + leja_legendre rule; + check(4.0 / 9, sparse_grid([](const point& x) { return x[0] * x[0] * x[1] * x[1]; }, rule, 2, 100000, tolerance, &error, &L1)); + check((2 * std::sin(1.0) - 2) * (2 * std::sin(1.0) - 2), + sparse_grid([](const point& x) { return (std::cos(x[0]) - 1) * (std::cos(x[1]) - 1); }, rule, 2, 100000, tolerance, &error, &L1)); + check(4 / (e * e), sparse_grid([](const point& x) { return x[0] * x[1] * std::exp(x[0] + x[1]); }, rule, 2, 100000, tolerance, &error, &L1)); + check(8 / (e * e * e), sparse_grid([](const point& x) { return x[0] * x[1] * x[2] * std::exp(x[0] + x[1] + x[2]); }, rule, 3, 100000, + tolerance, &error, &L1)); + check(boost::math::constants::pi() / 4, + sparse_grid([](const point& x) { return x[0] * x[0] * x[1] * x[1]; }, leja_hermite(), 2, 100000, tolerance, &error, &L1)); + // 0 everywhere: 0 after dimension + 1 rounds, which reach the product x_0^2 ... x_3^2. + for (unsigned d = 1; d <= 4; ++d) + { + unsigned calls = 0; + CHECK_EQUAL(sparse_grid([&](const point&) { ++calls; return 0.0; }, rule, d, 100000, tolerance, &error, &L1), 0.0); + CHECK_LE(calls, 1000u); + CHECK_EQUAL(error, 0.0); + CHECK_EQUAL(L1, 0.0); + } + check(16.0 / 81, sparse_grid([](const point& x) { return x[0] * x[0] * x[1] * x[1] * x[2] * x[2] * x[3] * x[3]; }, rule, 4, 100000, + tolerance, &error, &L1)); +} + +// Whether integrate(f, a, b) is declared. +template +struct has_bounds : std::false_type {}; +template +struct has_bounds().integrate(std::declval(), 0.0, 2.0))>> : std::true_type {}; +static_assert(has_bounds>::value && has_bounds>::value, ""); +static_assert(!has_bounds>::value && !has_bounds>::value, ""); + +// integrate(f, a, b) and sparse_grid(f, rule, a, b) map [-1, 1] affinely onto [a, b]. +template > +void test_bounds(const Real& tolerance) +{ + using boost::math::quadrature::sparse_grid; + using std::abs; + using std::cos; + using std::exp; + using std::pow; + using std::sin; + const Real eps = std::numeric_limits::epsilon(); + const leja_legendre rule; + Real error, L1, error2, L2; + auto check = [&](const Real& expected, const Real& Q) { + CHECK_LE(error, tolerance * L1); + CHECK_LE(abs(Q - expected), (std::max)(error, 64 * eps * L1)); + }; + // Every form resolves to the intended overload. + auto e = [](const Real& x) { return exp(x); }; + const Real e1 = exp(Real(1)) - exp(Real(-1)), e02 = exp(Real(2)) - 1; + const Real bound = (std::max)(boost::math::tools::root_epsilon(), Real(1e-8)); + CHECK_LE(abs(rule.integrate(e) - e1), bound); + CHECK_LE(abs(rule.integrate(e, 1e-10) - e1), bound); + CHECK_LE(abs(rule.integrate(e, 1e-10, &error, &L1) - e1), bound); + CHECK_LE(abs(rule.integrate(e, 0.0, 2.0) - e02), bound * e02); + CHECK_LE(abs(rule.integrate(e, 0.0, 2.0, 1e-10) - e02), bound * e02); + CHECK_LE(abs(rule.integrate(e, 0.0, 2.0, 1e-10, &error, &L1) - e02), bound * e02); + CHECK_LE(abs(rule.integrate(e, 0, 2) - e02), bound * e02); + + // A polynomial of degree 4 on [0, 3], exact from 5 points, and exp on [-2, 5]. + auto p = [](const Real& x) { return 1 + x - 2 * x * x + x * x * x * x; }; + check(Real(381) / 10, rule.integrate(p, Real(0), Real(3), tolerance, &error, &L1)); + const Real Q = rule.integrate(e, Real(-2), Real(5), tolerance, &error, &L1); + check(exp(Real(5)) - exp(Real(-2)), Q); + // The error and L1 of the rule for f(x(t)) on [-1, 1], times |b - a|/2. + const Real Qt = rule.integrate([](const Real& t) { return exp(Real(1.5) + Real(3.5) * t); }, tolerance, &error2, &L2); + CHECK_ULP_CLOSE(Real(3.5) * Qt, Q, 4); + CHECK_ULP_CLOSE(Real(3.5) * error2, error, 4); + CHECK_ULP_CLOSE(Real(3.5) * L2, L1, 4); + // a > b: minus the integral over [b, a], with a positive error and L1. + const Real R = rule.integrate(e, Real(5), Real(-2), tolerance, &error2, &L2); + CHECK_LE(error2, tolerance * L2); + CHECK_LE(abs(R + exp(Real(5)) - exp(Real(-2))), (std::max)(error2, 64 * eps * L2)); + // a == b: 0, of the shape of f(a), from one call. + unsigned calls = 0; + CHECK_EQUAL(rule.integrate([&](const Real& x) { ++calls; return x; }, Real(1.5), Real(1.5), tolerance, &error, &L1), Real(0)); + CHECK_EQUAL(calls, 1u); + CHECK_EQUAL(error, Real(0)); + CHECK_EQUAL(L1, Real(0)); + const std::vector Z = rule.integrate([](const Real& x) { return std::vector{x, x * x}; }, Real(1), Real(1)); + CHECK_EQUAL(Z.size(), std::size_t(2)); + CHECK_EQUAL(Z[0], Real(0)); + CHECK_EQUAL(Z[1], Real(0)); + const Real inf = std::numeric_limits::infinity(); + CHECK_THROW(rule.integrate(e, Real(0), inf), std::domain_error); + CHECK_THROW(rule.integrate(e, -inf, Real(0)), std::domain_error); + CHECK_THROW(rule.integrate(e, std::numeric_limits::quiet_NaN(), Real(0)), std::domain_error); + + // Complex and vector f: e^{ix} over [0, pi] is 2i, and (x, x^2) over [0, 2] is (2, 8/3). + const Real pi = boost::math::constants::pi(); + const Complex C = rule.integrate([](const Real& x) { return Complex(cos(x), sin(x)); }, Real(0), pi, tolerance, &error, &L1); + CHECK_LE(Real(abs(C - Complex(0, 2))), (std::max)(error, 64 * eps * L1)); + const std::array V = rule.integrate([](const Real& x) { return std::array{x, x * x}; }, Real(0), Real(2), tolerance, &error, &L1); + check(Real(2), V[0]); + check(Real(8) / 3, V[1]); + + // The integral of x^k (4 - x)^alpha (x - 1)^beta over [1, 4] is 3^(alpha + beta + 1) B(alpha + 1, beta + 1) for k = 0, + // and that times 1 + 3 (beta + 1)/(alpha + beta + 2) for k = 1. Over [4, 1], alpha goes with 1, and the sign flips. + for (Real alpha : {Real(0.5), Real(2)}) + { + for (Real beta : {Real(-0.25), Real(1.5)}) + { + const leja_jacobi jacobi(alpha, beta); + const Real m = pow(Real(3), alpha + beta + 1) * boost::math::beta(alpha + 1, beta + 1); + const Real mean = 1 + 3 * (beta + 1) / (alpha + beta + 2); + check(m, jacobi.integrate([](const Real&) { return Real(1); }, Real(1), Real(4), tolerance, &error, &L1)); + check(m * mean, jacobi.integrate([](const Real& x) { return x; }, Real(1), Real(4), tolerance, &error, &L1)); + check(-m * (5 - mean), jacobi.integrate([](const Real& x) { return x; }, Real(4), Real(1), tolerance, &error, &L1)); + // The error and L1 scale by (3/2)^(alpha + beta + 1). + const Real Qj = jacobi.integrate(e, Real(1), Real(4), tolerance, &error, &L1); + const Real s = pow(Real(1.5), alpha + beta + 1); + CHECK_ULP_CLOSE(s * jacobi.integrate([](const Real& t) { return exp(Real(2.5) + Real(1.5) * t); }, tolerance, &error2, &L2), Qj, 8); + CHECK_ULP_CLOSE(s * error2, error, 8); + CHECK_ULP_CLOSE(s * L2, L1, 8); + // x_0 in a sparse grid on [1, 4]^2. + check(m * m * mean, sparse_grid([](const std::vector& x) { return x[0]; }, jacobi, std::vector{1, 1}, + std::vector{4, 4}, 10000, tolerance, &error, &L1)); + } + } + + // exp(x_0/2 - x_1 + x_2/4) on [0, 1] x [-1, 2] x [3, 5] with std::array bounds and points, and x_0 x_1 x_2 on [0, 1]^3. + auto box = [](const Real& c, const Real& a, const Real& b) { return (exp(c * b) - exp(c * a)) / c; }; + const std::array a{0, -1, 3}, b{1, 2, 5}; + auto g = [](const std::array& x) { return exp(x[0] / 2 - x[1] + x[2] / 4); }; + const Real G = sparse_grid(g, rule, a, b, 100000, tolerance, &error, &L1); + check(box(Real(0.5), 0, 1) * box(Real(-1), -1, 2) * box(Real(0.25), 3, 5), G); + // The result, error and L1 of the grid on [-1, 1]^3 for g(x(t)), times the volume / 8. + const Real Gt = sparse_grid([&](const std::array& t) { + return g(std::array{Real(0.5) + t[0] / 2, Real(0.5) + Real(1.5) * t[1], 4 + t[2]}); + }, rule, 3, 100000, tolerance, &error2, &L2); + CHECK_ULP_CLOSE(Real(0.75) * Gt, G, 8); + CHECK_ULP_CLOSE(Real(0.75) * error2, error, 8); + CHECK_ULP_CLOSE(Real(0.75) * L2, L1, 8); + check(Real(1) / 8, sparse_grid([](const std::vector& x) { return x[0] * x[1] * x[2]; }, rule, std::vector(3, 0), + std::vector(3, 1), 1000, tolerance, &error, &L1)); + // A complex f and a reversed coordinate: the integrals of e^{ix} over [0, pi] and over [0, -pi/2] are 2i and -1 + i. + const Complex S = sparse_grid([](const std::vector& x) { return Complex(cos(x[0] + x[1]), sin(x[0] + x[1])); }, rule, + std::vector{0, 0}, std::vector{pi, -pi / 2}, 10000, tolerance, &error, &L1); + CHECK_LE(Real(abs(S - Complex(0, 2) * Complex(-1, 1))), (std::max)(error, 64 * eps * L1)); + // A box of zero volume, even where the Jacobi exponents are negative. + CHECK_EQUAL(sparse_grid(g, leja_jacobi(Real(-0.25), Real(-0.25)), a, std::array{1, -1, 5}, 1000, tolerance, &error, &L1), Real(0)); + CHECK_EQUAL(error, Real(0)); + CHECK_EQUAL(L1, Real(0)); + const std::vector Zg = sparse_grid([](const std::vector& x) { return std::vector{x[0], x[1]}; }, rule, + std::vector{0, 1}, std::vector{1, 1}, 1000, tolerance, &error, &L1); + CHECK_EQUAL(Zg.size(), std::size_t(2)); + CHECK_EQUAL(Zg[0], Real(0)); + CHECK_EQUAL(Zg[1], Real(0)); + CHECK_EQUAL(error, Real(0)); + CHECK_EQUAL(L1, Real(0)); + CHECK_THROW(sparse_grid(g, rule, a, std::array{1, 2, inf}, 1000), std::domain_error); + CHECK_THROW(sparse_grid(g, rule, std::vector(3, 0), std::vector(2, 1), 1000), std::domain_error); + CHECK_THROW(sparse_grid(g, rule, std::vector(2, 0), std::vector(2, 1), 1000), std::domain_error); + CHECK_THROW(sparse_grid(g, rule, std::vector(2, 0), std::vector(2, 0), 1000), std::domain_error); +} + +// Threads that share a rule, before it has any points, get what a rule of their own gives them. +void test_threads() +{ + const leja_hermite shared; + std::vector Q(8); + std::vector threads; + for (unsigned t = 0; t < Q.size(); ++t) + threads.emplace_back([&, t] { Q[t] = shared.integrate([t](double x) { return std::cos((t + 1) * x / 2); }, 1e-12); }); + for (std::thread& thread : threads) + thread.join(); + // The compiler may contract the sums into fused multiply-adds differently where they are inlined. + const leja_hermite alone; + for (unsigned t = 0; t < Q.size(); ++t) + CHECK_ULP_CLOSE(Q[t], alone.integrate([t](double x) { return std::cos((t + 1) * x / 2); }, 1e-12), 2); + const unsigned n = static_cast(shared.abscissa().size()); + for (unsigned i = 0; i < n; ++i) + CHECK_EQUAL(shared.weights(n)[i], alone.weights(n)[i]); +} + +void test_errors() +{ + CHECK_THROW(leja_laguerre(-1.0), std::domain_error); + CHECK_THROW(leja_jacobi(0.5, -1.5), std::domain_error); + CHECK_THROW(leja_jacobi(0.0, std::numeric_limits::quiet_NaN()), std::domain_error); + const double inf = std::numeric_limits::infinity(); + CHECK_THROW(leja_laguerre(+inf), std::domain_error); + CHECK_THROW(leja_jacobi(inf, 0.0), std::domain_error); + CHECK_THROW(leja_jacobi(0.0, inf), std::domain_error); + leja_legendre rule; + CHECK_THROW(rule.weights(0), std::domain_error); + CHECK_THROW(rule.weights(rule.max_points() + 1), std::domain_error); + // A NaN at the fifth point, 0.66, reaches the error estimate. + double error; + rule.integrate([](double x) { return x > 0 && x < 1 ? std::numeric_limits::quiet_NaN() : std::exp(x); }, 1e-10, &error); + CHECK_NAN(error); + auto one = [](const std::vector&) { return 1.0; }; + CHECK_THROW(boost::math::quadrature::sparse_grid(one, rule, 0, 10), std::domain_error); + CHECK_THROW(boost::math::quadrature::sparse_grid(one, rule, 2, 0), std::domain_error); + CHECK_THROW(boost::math::quadrature::sparse_grid([](const std::array&) { return 1.0; }, rule, 2, 10), std::domain_error); +} + +// With domain errors ignored, an invalid n or exponent gives no weights and a NaN integral. +void test_ignore_error() +{ + using ignore = boost::math::policies::policy>; + auto x = [](double t) { return t; }; + double error, L1; + const leja_legendre rule; + CHECK_EQUAL(rule.weights(0).size(), std::size_t(0)); + CHECK_EQUAL(rule.weights(rule.max_points() + 1).size(), std::size_t(0)); + const leja_laguerre laguerre(-1.0); + const leja_jacobi jacobi(0.5, -1.5); + CHECK_EQUAL(laguerre.weights(1).size(), std::size_t(0)); + CHECK_EQUAL((leja_jacobi(std::numeric_limits::infinity(), 0.0).weights(1).size()), std::size_t(0)); + CHECK_NAN(laguerre.integrate(x, 1e-10, nullptr, &L1)); + CHECK_NAN(L1); + CHECK_NAN(jacobi.integrate(x, 1e-10, &error, &L1)); + CHECK_NAN(error); + CHECK_NAN(L1); + CHECK_EQUAL(jacobi.abscissa().size(), std::size_t(0)); +} + +int main() +{ + test_exactness(); + test_exactness(); + test_exactness(); + test_exactness(20); + + test_nested(); + test_nested(); + + test_same_points(128); + test_same_points(100); + test_same_points(128); + test_same_points(100); + test_same_points(40); + test_same_points(40); + test_same_points(100); + test_same_points(100); + test_same_points(40); + test_same_points(100, 2.5); + test_same_points(100, -0.5, 2.0); + test_same_points(40, 0.5, -0.75); + + test_smooth(30); + // long double is quad precision on some platforms, aarch64 Linux among them, and needs more points. + const bool quad = std::numeric_limits::digits > 64; + test_smooth(quad ? 40 : 30); + + test_adaptive(); + test_adaptive(); + + test_complex(20, 30); + test_complex(20, 30); + test_complex(quad ? 30 : 20, quad ? 50 : 30); + test_complex(36, 64); + + test_vector(); + test_vector(); + test_vector(); + + test_sparse_grid(1e-4f); + test_sparse_grid(1e-10); + test_sparse_grid(1e-12); + test_vanishing(); + + test_bounds(1e-4f); + test_bounds(1e-10); +#ifdef BOOST_HAS_FLOAT128 + test_bounds(1e-20); +#endif + test_bounds(1e-20); + + test_threads(); + test_errors(); + test_ignore_error(); + return boost::math::test::report_errors(); +}