From 01642da17dc5c5061a96bb70ba95a1b67678ba24 Mon Sep 17 00:00:00 2001 From: MatejTomes Date: Fri, 17 Jan 2025 16:19:22 +0100 Subject: [PATCH 1/4] Add 2D Gauss-Legendre quadrature integrator --- cherab/core/math/integrators/__init__.pxd | 2 +- cherab/core/math/integrators/__init__.py | 2 +- .../core/math/integrators/integrators2d.pxd | 17 + .../core/math/integrators/integrators2d.pyx | 310 +++++++++++++++++- 4 files changed, 328 insertions(+), 3 deletions(-) diff --git a/cherab/core/math/integrators/__init__.pxd b/cherab/core/math/integrators/__init__.pxd index c4eff464..4b37ec97 100644 --- a/cherab/core/math/integrators/__init__.pxd +++ b/cherab/core/math/integrators/__init__.pxd @@ -17,5 +17,5 @@ # under the Licence. from cherab.core.math.integrators.integrators1d cimport Integrator1D, GaussianQuadrature -from cherab.core.math.integrators.integrators2d cimport Integrator2D +from cherab.core.math.integrators.integrators2d cimport Integrator2D, GaussianQuadrature2D diff --git a/cherab/core/math/integrators/__init__.py b/cherab/core/math/integrators/__init__.py index b1fa4516..1a0de26d 100644 --- a/cherab/core/math/integrators/__init__.py +++ b/cherab/core/math/integrators/__init__.py @@ -17,4 +17,4 @@ # under the Licence. from .integrators1d import Integrator1D, GaussianQuadrature -from .integrators2d import Integrator2D +from .integrators2d import Integrator2D, GaussianQuadrature2D diff --git a/cherab/core/math/integrators/integrators2d.pxd b/cherab/core/math/integrators/integrators2d.pxd index 62ba8667..5a6147ac 100644 --- a/cherab/core/math/integrators/integrators2d.pxd +++ b/cherab/core/math/integrators/integrators2d.pxd @@ -16,6 +16,8 @@ # See the Licence for the specific language governing permissions and limitations # under the Licence. +from numpy cimport ndarray + from raysect.core.math.function.float cimport Function1D, Function2D @@ -25,3 +27,18 @@ cdef class Integrator2D: Function2D function cdef double evaluate(self,double x_lower, double x_upper, Function1D y_lower, Function1D y_upper) except? -1e999 + + +cdef class GaussianQuadrature2D(Integrator2D): + + cdef: + int _x_min_order, _x_max_order, _y_min_order, _y_max_order + double _rtol + ndarray _x_roots, _x_weights, _y_roots, _y_weights + double[:] _x_roots_mv, _x_weights_mv, _y_roots_mv, _y_weights_mv + + cdef _build_cache(self) + + cpdef double profile_evaluate(self, int n, double x_lower, double x_upper, Function1D y_lower, Function1D y_upper) + + cpdef double evaluate_overhead(self, int n) \ No newline at end of file diff --git a/cherab/core/math/integrators/integrators2d.pyx b/cherab/core/math/integrators/integrators2d.pyx index 58627afe..05a78db3 100644 --- a/cherab/core/math/integrators/integrators2d.pyx +++ b/cherab/core/math/integrators/integrators2d.pyx @@ -18,7 +18,13 @@ # See the Licence for the specific language governing permissions and limitations # under the Licence. -from raysect.core.math.function.float cimport Function1D, autowrap_function2d +import numpy as np +from scipy.special import roots_legendre +cimport cython + +from raysect.core.math.function.float cimport Function1D, autowrap_function2d, Constant2D + +from libc.math cimport INFINITY cdef class Integrator2D: @@ -59,3 +65,305 @@ cdef class Integrator2D: """ return self.evaluate(x_lower, x_upper, y_lower, y_upper) + + +cdef class GaussianQuadrature2D(Integrator2D): + """ + Approximates an integral of a two-dimensional function over a finite interval. + + The integral is approximated with fixed-tolerance Gauss-Legendre quadrature. + The quadrature approximation is calculated as follows: + + .. math:: + \\int_{x_{\\text{lower}}}^{x_{\\text{upper}}} \\int_{y_{\\text{lower}}(x)}^{y_{\\text{upper}}(x)} f(x, y) \\, dy \\, dx + \\approx \\sum_{i=1}^{n_x} \\sum_{j=1}^{n_y} w_{i} w_{j} \\, + f\\left( B \\xi_i + A, \\, D(x) \\eta_j + C(x) \\right) \\, + B \\cdot D(x) + where: + - :math:`x_{\mathrm{lower}}`: Lower limit of integration for the x-dimension. + - :math:`x_{\mathrm{upper}}`: Upper limit of integration for the x-dimension. + - :math:`y_{\mathrm{lower}}(x)`: Lower limit of integration for the y-dimension, a function of :math:`x`. + - :math:`y_{\mathrm{upper}}(x)`: Upper limit of integration for the y-dimension, a function of :math:`x`. + - :math:`f(x, y)`: The function to be integrated over the specified region. + - :math:`\\xi`: The transformed variable for the x-dimension, ranging from -1 to 1. + - :math:`\\eta`: The transformed variable for the y-dimension, ranging from -1 to 1. + - :math:`w_i`: The weight corresponding to the :math:`i`-th root of the Legendre polynomial in the x-dimension. + - :math:`w_j`: The weight corresponding to the :math:`j`-th root of the Legendre polynomial in the y-dimension. + - :math:`\\xi_i`: The :math:`i`-th root of the Legendre polynomial in the x-dimension. + - :math:`\\eta_j`: The :math:`j`-th root of the Legendre polynomial in the y-dimension. + - :math:`n_x`: The number of roots (or nodes) in the x-dimension. + - :math:`n_y`: The number of roots (or nodes) in the y-dimension. + - :math:`A = \\frac{x_{\\mathrm{upper}} + x_{\\mathrm{lower}}}{2}`: Midpoint of the x-interval. + - :math:`B = \\frac{x_{\\mathrm{upper}} - x_{\\mathrm{lower}}}{2}`: Half-width of the x-interval. + - :math:`D(x) = \\frac{y_{\\mathrm{upper}}(x) - y_{\\mathrm{lower}}(x)}{2}`: Half-width of the y-interval. + - :math:`C(x) = \\frac{y_{\\mathrm{upper}}(x) + y_{\\mathrm{lower}}(x)}{2}`: Midpoint of the y-interval. + + :param Function2D integrand: A 2D function to integrate. Default is `Constant2D(0)`. + :param double relative_tolerance: Iteration stops when relative error between + last two iterates is less than this value. Default is `1.e-5`. + :param int x_max_order: Maximum order on Gaussian quadrature in the x dimension. Default is `50`. + :param int x_min_order: Minimum order on Gaussian quadrature in the x dimension. Default is `1`. + :param int y_max_order: Maximum order on Gaussian quadrature in the y dimension. Default is `50`. + :param int y_min_order: Minimum order on Gaussian quadrature in the y dimension. Default is `1`. + + :ivar Function1D integrand: A 1D function to integrate. + :ivar double relative_tolerance: Iteration stops when relative error between + last two iterates is less than this value. + :ivar int x_max_order: Maximum order on Gaussian quadrature in the x dimension. + :ivar int x_min_order: Minimum order on Gaussian quadrature in the x dimension. + :ivar int y_max_order: Maximum order on Gaussian quadrature in the y dimension. + :ivar int y_min_order: Minimum order on Gaussian quadrature in the y dimension. + """ + def __init__(self, object integrand=Constant2D(0), double relative_tolerance=1.e-5, int x_max_order=50, int x_min_order=1, + int y_max_order=50, int y_min_order=1): + + self._check_order(x_min_order, x_max_order, "x") + self._check_order(y_min_order, y_max_order, "y") + + self._x_min_order = x_min_order + self._x_max_order = x_max_order + + self._y_min_order = y_min_order + self._y_max_order = y_max_order + + self._build_cache() + + self.integrand = integrand + + self.relative_tolerance = relative_tolerance + + def _check_order(self, int min_order, int max_order, str dimension): + """ + Check the order of Gaussian quadrature. + + :param int min_order: Minimum order on Gaussian quadrature. + :param int max_order: Maximum order on Gaussian quadrature. + :param str dimension: Dimension of Gaussian quadrature. + + :raises ValueError: If the order of Gaussian quadrature is invalid. + """ + + if min_order < 1 or max_order < 1: + raise ValueError("Order of Gaussian quadrature in the {} dimension must be >= 1.".format(dimension)) + + if min_order > max_order: + raise ValueError("Minimum order of Gaussian quadrature in the {} dimension must be less than or equal to the maximum order.".format(dimension)) + + @property + def x_min_order(self): + """ + Minimum order on Gaussian quadrature in the x dimension. + + :rtype: int + """ + return self._x_min_order + + @x_min_order.setter + def x_min_order(self, int value): + + self._check_order(value, self._x_max_order, "x") + + self._x_min_order = value + + self._build_cache() + + @property + def x_max_order(self): + """ + Maximum order on Gaussian quadrature in the x dimension. + + :rtype: int + """ + return self._x_max_order + + @x_max_order.setter + def x_max_order(self, int value): + + self._check_order(self._x_min_order, value, "x") + + self._x_max_order = value + + self._build_cache() + + @property + def y_min_order(self): + """ + Minimum order on Gaussian quadrature in the y dimension. + + :rtype: int + """ + return self._y_min_order + + @y_min_order.setter + def y_min_order(self, int value): + + self._check_order(value, self._y_max_order, "y") + + self._y_min_order = value + + self._build_cache() + + @property + def y_max_order(self): + """ + Maximum order on Gaussian quadrature in the y dimension. + + :rtype: int + """ + return self._y_max_order + + @y_max_order.setter + def y_max_order(self, int value): + + self._check_order(self._y_min_order, value, "y") + + self._y_max_order = value + + self._build_cache() + + @property + def relative_tolerance(self): + """ + Iteration stops when relative error between last two iterates is less than this value. + + :rtype: double + """ + return self._rtol + + @relative_tolerance.setter + def relative_tolerance(self, double value): + + if value <= 0: + raise ValueError("Relative tolerance must be positive.") + + self._rtol = value + + cdef _build_cache(self): + """ + Caches the roots and weights of the Gauss-Legendre quadrature. + """ + + cdef: + int order, n, i + + # cache roots and weights for x dimension + n = (self._x_max_order + self._x_min_order) * (self._x_max_order - self._x_min_order + 1) // 2 + + self._x_roots = np.zeros(n, dtype=np.float64) + self._x_weights = np.zeros(n, dtype=np.float64) + + i = 0 + for order in range(self._x_min_order, self._x_max_order + 1): + self._x_roots[i:i + order], self._x_weights[i:i + order] = roots_legendre(order) + i += order + + self._x_roots_mv = self._x_roots + self._x_weights_mv = self._x_weights + + # cache roots and weights for y dimension + n = (self._y_max_order + self._y_min_order) * (self._y_max_order - self._y_min_order + 1) // 2 + + self._y_roots = np.zeros(n, dtype=np.float64) + self._y_weights = np.zeros(n, dtype=np.float64) + + i = 0 + for order in range(self._y_min_order, self._y_max_order + 1): + self._y_roots[i:i + order], self._y_weights[i:i + order] = roots_legendre(order) + i += order + + self._y_roots_mv = self._y_roots + self._y_weights_mv = self._y_weights + + cpdef double profile_evaluate(self, int n, double x_lower, double x_upper, Function1D y_lower, Function1D y_upper): + + cdef: + int i + + for i in range(n): + self.evaluate(x_lower, x_upper, y_lower, y_upper) + + return 0. + + cpdef double evaluate_overhead(self, int n): + cdef: + int i + + for i in range(n): + continue + + return 0. + + @cython.boundscheck(False) + @cython.wraparound(False) + @cython.cdivision(True) + @cython.initializedcheck(False) + cdef double evaluate(self, double x_lower, double x_upper, Function1D y_lower, Function1D y_upper) except? -1e999: + """ + Integrates a two-dimensional function over a finite interval. + + :param double x_lower: Lower limit of integration in the x dimension. + :param double x_upper: Upper limit of integration in the x dimension. + :param Function1D y_lower: Lower limit of integration in the y dimension as a function of x. + :param Function1D y_upper: Upper limit of integration in the y dimension as a function of x. + + :returns: Gaussian quadrature approximation to integral. + """ + + cdef: + int x_order, y_order + int x_ibegin, y_ibegin + int i, j + double current_integral, previous_integral, y_contribution, error + double x, y, x_offset, y_offset, x_slope, y_slope, scale, rtol, y_lower_val, y_upper_val + + # set initial iteration values + previous_integral = INFINITY + x_ibegin = 0 + y_ibegin = 0 + + # x variable transformations from [-1, 1] to [x_lower, x_upper] and [y_lower, y_upper] + x_offset = 0.5 * (x_lower + x_upper) + x_slope = 0.5 * (x_upper - x_lower) + + # create local caches + x_order = self._x_min_order + y_order = self._y_min_order + + while True: + current_integral = 0. + for i in range(x_ibegin, x_ibegin + x_order): + x = x_offset + x_slope * self._x_roots_mv[i] + y_contribution = 0. + + # calculate integration limits and transformation from [-1, 1] to [y_lower, y_upper] + y_lower_val = y_lower.evaluate(x) + y_upper_val = y_upper.evaluate(x) + y_offset = 0.5 * (y_lower_val + y_upper_val) + y_slope = 0.5 * (y_upper_val - y_lower_val) + + for j in range(y_ibegin, y_ibegin + y_order): + y = y_offset + y_slope * self._y_roots_mv[j] + y_contribution += self._y_weights_mv[j] * self.function.evaluate(x, y) + + current_integral += self._x_weights_mv[i] * y_contribution + current_integral *= x_slope * y_slope + + error = abs(current_integral - previous_integral) + + # terminate integration if relative error is less than tolerance + # or if both x and y maximum orders have been reached + if error < self._rtol * abs(current_integral) or (x_order == self._x_max_order and y_order == self._y_max_order): + break + + previous_integral = current_integral + + # increase the orders of x and y and shift indexes + if x_order < self._x_max_order: + x_ibegin += x_order + x_order += 1 + + if y_order < self._y_max_order: + y_ibegin += y_order + y_order += 1 + + return current_integral From 4ce86493f1d4bda5bc482dc0167fc5492f078146 Mon Sep 17 00:00:00 2001 From: MatejTomes Date: Fri, 17 Jan 2025 16:20:11 +0100 Subject: [PATCH 2/4] Add tests for GaussianQuadrature2D --- cherab/core/math/tests/test_integrators.py | 128 +++++++++++++++++++-- 1 file changed, 119 insertions(+), 9 deletions(-) diff --git a/cherab/core/math/tests/test_integrators.py b/cherab/core/math/tests/test_integrators.py index fa967165..9eb6d855 100644 --- a/cherab/core/math/tests/test_integrators.py +++ b/cherab/core/math/tests/test_integrators.py @@ -16,8 +16,8 @@ # See the Licence for the specific language governing permissions and limitations # under the Licence. -from raysect.core.math.function.float import Exp1D, Arg1D -from cherab.core.math.integrators import GaussianQuadrature +from raysect.core.math.function.float import Exp1D, Arg1D, Exp2D, Arg2D, Constant1D +from cherab.core.math.integrators import GaussianQuadrature, GaussianQuadrature2D from math import sqrt, pi from scipy.special import erf import unittest @@ -30,8 +30,13 @@ def test_properties(self): """Test property assignment.""" min_order = 3 max_order = 30 - reltol = 1.e-6 - quadrature = GaussianQuadrature(integrand=Arg1D, relative_tolerance=reltol, max_order=max_order, min_order=min_order) + reltol = 1.0e-6 + quadrature = GaussianQuadrature( + integrand=Arg1D, + relative_tolerance=reltol, + max_order=max_order, + min_order=min_order, + ) self.assertEqual(quadrature.relative_tolerance, reltol) self.assertEqual(quadrature.max_order, max_order) @@ -53,7 +58,7 @@ def test_properties(self): min_order = 1 max_order = 20 - reltol = 1.e-5 + reltol = 1.0e-5 quadrature.relative_tolerance = reltol quadrature.min_order = min_order @@ -67,14 +72,119 @@ def test_properties(self): def test_integrate(self): """Test integration.""" - quadrature = GaussianQuadrature(relative_tolerance=1.e-8) + quadrature = GaussianQuadrature(relative_tolerance=1.0e-8) a = -0.5 - b = 3. - quadrature.integrand = (2 / sqrt(pi)) * Exp1D(- Arg1D() * Arg1D()) + b = 3.0 + quadrature.integrand = (2 / sqrt(pi)) * Exp1D(-Arg1D() * Arg1D()) exact_integral = erf(b) - erf(a) self.assertAlmostEqual(quadrature(a, b), exact_integral, places=8) -if __name__ == '__main__': +class TestGaussianQuadrature2D(unittest.TestCase): + """Gaussian quadrature 2D integrator tests.""" + + def test_properties(self): + """Test property assignment.""" + x_min_order = 3 + y_min_order = 4 + x_max_order = 30 + y_max_order = 40 + reltol = 1.0e-6 + integrand = Exp2D(Arg2D("x") + Arg2D("y")) + + quadrature = GaussianQuadrature2D( + integrand=integrand, + relative_tolerance=reltol, + x_max_order=x_max_order, + x_min_order=x_min_order, + y_max_order=y_max_order, + y_min_order=y_min_order, + ) + + self.assertEqual(quadrature.relative_tolerance, reltol) + self.assertEqual(quadrature.x_max_order, x_max_order) + self.assertEqual(quadrature.y_max_order, y_max_order) + self.assertEqual(quadrature.x_min_order, x_min_order) + self.assertEqual(quadrature.y_min_order, y_min_order) + self.assertEqual(quadrature.integrand, integrand) + + x_min_order = 50 # > x_max_order + x_max_order = 2 # < x_min_order + y_min_order = 50 # > y_max_order + y_max_order = 1 # < y_min_order + reltol = -1 + + with self.assertRaises(ValueError): + quadrature.x_max_order = x_max_order + + with self.assertRaises(ValueError): + quadrature.x_min_order = x_min_order + + with self.assertRaises(ValueError): + quadrature.y_max_order = y_max_order + + with self.assertRaises(ValueError): + quadrature.y_min_order = y_min_order + + with self.assertRaises(ValueError): + quadrature.relative_tolerance = reltol + + x_min_order = 0 + y_min_order = 0 + + with self.assertRaises(ValueError): + quadrature.x_min_order = x_min_order + + with self.assertRaises(ValueError): + quadrature.y_min_order = y_min_order + + x_min_order = 1 + x_max_order = 20 + y_min_order = 2 + y_max_order = 30 + reltol = 1.0e-5 + + quadrature.relative_tolerance = reltol + quadrature.x_min_order = x_min_order + quadrature.x_max_order = x_max_order + quadrature.y_min_order = y_min_order + quadrature.y_max_order = y_max_order + quadrature.integrand = Arg2D("x") + + self.assertEqual(quadrature.relative_tolerance, reltol) + self.assertEqual(quadrature.x_min_order, x_min_order) + self.assertEqual(quadrature.x_max_order, x_max_order) + self.assertEqual(quadrature.y_min_order, y_min_order) + self.assertEqual(quadrature.y_max_order, y_max_order) + self.assertEqual(quadrature.integrand, Arg2D("x")) + + def test_integrate(self): + """Test 2D integration.""" + quadrature = GaussianQuadrature2D(relative_tolerance=1.0e-8) + + # Integration limits + a_x, b_x = -2.0, 2.0 + a_y, b_y = -3.0, 3.0 + + # Bivariate Normal distribution with std_dev=1, mean=0 and no correlation + quadrature.integrand = ( + 1 / (2 * pi) * Exp2D(-0.5 * (Arg2D("x") ** 2 + Arg2D("y") ** 2)) + ) + + # Exact integral of the bivariate normal distribution + exact_integral = ( + 1 / 4.0 + * (erf(b_x / sqrt(2)) - erf(a_x / sqrt(2))) + * (erf(b_y / sqrt(2)) - erf(a_y / sqrt(2))) + ) + + self.assertAlmostEqual( + quadrature(a_x, b_x, Constant1D(a_y), Constant1D(b_y)), + exact_integral, + places=8, + ) + + +if __name__ == "__main__": unittest.main() From 156509c8f05652009eee8b340242f562762f4f88 Mon Sep 17 00:00:00 2001 From: MatejTomes Date: Fri, 17 Jan 2025 16:20:50 +0100 Subject: [PATCH 3/4] Add Integrator2D documentation --- docs/source/math/integrators.rst | 11 +++++++++++ docs/source/math/math.rst | 1 + 2 files changed, 12 insertions(+) create mode 100644 docs/source/math/integrators.rst diff --git a/docs/source/math/integrators.rst b/docs/source/math/integrators.rst new file mode 100644 index 00000000..de9234fb --- /dev/null +++ b/docs/source/math/integrators.rst @@ -0,0 +1,11 @@ + +Integrators +------------- + +.. autoclass:: cherab.core.math.integrators.integrators1d.GaussianQuadrature + :members: + :show-inheritance: + +.. autoclass:: cherab.core.math.integrators.integrators2d.GaussianQuadrature2D + :members: + :show-inheritance: \ No newline at end of file diff --git a/docs/source/math/math.rst b/docs/source/math/math.rst index 836cbe84..e0a85775 100644 --- a/docs/source/math/math.rst +++ b/docs/source/math/math.rst @@ -22,3 +22,4 @@ utilities that Cherab provides for slicing, dicing and projecting these function mask samplers slice + integrators From b16bd80345e55faf72c2baba24be424d11620803 Mon Sep 17 00:00:00 2001 From: MatejTomes Date: Fri, 17 Jan 2025 16:45:53 +0100 Subject: [PATCH 4/4] Add CHANGELOG.md record --- CHANGELOG.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index e3647dad..f3ebb93a 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -5,7 +5,7 @@ Release 1.6.0 (TBD) ------------------- New: -* Add Integrator2D base class for integration of two-dimensional functions. (#472) +* Add Integrator2D base class for integration of two-dimensional functions. Add GaussianQuadrature2D integrator. (#472, #475) Release 1.5.0 (27 Aug 2024) -------------------