Add GaussianQuadrature2D

CHANGELOG.md

@@ -9,11 +9,12 @@ API changes:
 * Add emission model attribute access to line and lineshape . (#294)
 
 New:
+* Add GaussianQuadrature2D integrator. (#475)
+* Support Raysect 0.9. (#486)
+* Test against Python 3.9, 3.10, 3.11, 3.12, 3.13 and latest released Numpy. Drop Python 3.7, 3.8 and older Numpy from tests. (#486)
 * Add Function6D framework. (#478)
 * Add e_field attribute to Plasma object for electric field vector. (#465)
 * Add Integrator2D base class for integration of two-dimensional functions. (#472)
-* Support Raysect 0.9. (#486)
-* Test against Python 3.9, 3.10, 3.11, 3.12, 3.13 and latest released Numpy. Drop Python 3.7, 3.8 and older Numpy from tests. (#486)
 
 Release 1.5.0 (27 Aug 2024)
 -------------------

cherab/core/math/integrators/__init__.pxd

@@ -16,6 +16,6 @@
 # See the Licence for the specific language governing permissions and limitations
 # 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.integrators1d cimport Integrator1D, GaussianQuadrature1D, GaussianQuadrature
+from cherab.core.math.integrators.integrators2d cimport Integrator2D, GaussianQuadrature2D
 

cherab/core/math/integrators/__init__.py

@@ -16,5 +16,5 @@
 # See the Licence for the specific language governing permissions and limitations
 # under the Licence.
 
-from .integrators1d import Integrator1D, GaussianQuadrature
-from .integrators2d import Integrator2D
+from .integrators1d import Integrator1D, GaussianQuadrature1D, GaussianQuadrature
+from .integrators2d import Integrator2D, GaussianQuadrature2D

cherab/core/math/integrators/integrators1d.pxd

@@ -28,7 +28,7 @@ cdef class Integrator1D:
     cdef double evaluate(self, double a, double b) except? -1e999
 
 
-cdef class GaussianQuadrature(Integrator1D):
+cdef class GaussianQuadrature1D(Integrator1D):
 
     cdef:
         int _min_order, _max_order
@@ -37,3 +37,7 @@ cdef class GaussianQuadrature(Integrator1D):
         double[:] _roots_mv, _weights_mv
 
     cdef _build_cache(self)
+
+
+cdef class GaussianQuadrature(GaussianQuadrature1D):
+    pass

cherab/core/math/integrators/integrators1d.pyx

@@ -18,6 +18,8 @@
 # See the Licence for the specific language governing permissions and limitations
 # under the Licence.
 
+import warnings
+
 import numpy as np
 from scipy.special import roots_legendre
 
@@ -65,23 +67,25 @@ cdef class Integrator1D:
         return self.evaluate(a, b)
 
 
-cdef class GaussianQuadrature(Integrator1D):
+cdef class GaussianQuadrature1D(Integrator1D):
     """
     Compute an integral of a one-dimensional function over a finite interval
     using fixed-tolerance Gaussian quadrature.
     (see Scipy `quadrature <https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.quadrature.html>`).
 
+    The integration is performed by iteratively increasing the order of the Gaussian quadrature until the relative tolerance is met or the maximum order is reached.
+
     :param object integrand: A 1D function to integrate. Default is Constant1D(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 max_order: Maximum order on Gaussian quadrature. Default is 50.
-    :param int min_order: Minimum order on Gaussian quadrature. Default is 1.
+    :param int max_order: Maximum order on Gaussian quadrature the integration stops at. Default is 50.
+    :param int min_order: Minimum order on Gaussian quadrature the integration starts from. 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 max_order: Maximum order on Gaussian quadrature.
-    :ivar int min_order: Minimum order on Gaussian quadrature.
+    :ivar int max_order: Maximum order on Gaussian quadrature the integration stops at.
+    :ivar int min_order: Minimum order on Gaussian quadrature the integration starts from.
     """
 
     def __init__(self, object integrand=Constant1D(0), double relative_tolerance=1.e-5, int max_order=50, int min_order=1):
@@ -169,6 +173,8 @@ cdef class GaussianQuadrature(Integrator1D):
         cdef:
             int order, n, i
 
+        # Store the variable-length roots and weights for each quadrature order
+        # consecutively in packed 1D arrays to avoid rectangular-array padding.
         n = (self._max_order + self._min_order) * (self._max_order - self._min_order + 1) // 2
 
         self._roots = np.zeros(n, dtype=np.float64)
@@ -222,3 +228,31 @@ cdef class GaussianQuadrature(Integrator1D):
                 break
 
         return newval
+
+
+cdef class GaussianQuadrature(GaussianQuadrature1D):
+    """
+    Compute an integral of a one-dimensional function over a finite interval
+    using fixed-tolerance Gaussian quadrature.
+    (see Scipy `quadrature <https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.quadrature.html>`).
+
+    .. warning::
+        This class is deprecated and will be removed in cherab 1.7. Use :class:`GaussianQuadrature1D` instead.
+
+    :param object integrand: A 1D function to integrate. Default is Constant1D(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 max_order: Maximum order on Gaussian quadrature. Default is 50.
+    :param int min_order: Minimum order on Gaussian quadrature. 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 max_order: Maximum order on Gaussian quadrature.
+    :ivar int min_order: Minimum order on Gaussian quadrature.
+    """
+
+    def __init__(self, object integrand=Constant1D(0), double relative_tolerance=1.e-5, int max_order=50, int min_order=1):
+
+        warnings.warn("The GaussianQuadrature class is deprecated and will be removed in cherab 1.7. Use GaussianQuadrature1D instead.", DeprecationWarning, stacklevel=2)
+        super().__init__(integrand, relative_tolerance, max_order, min_order)

cherab/core/math/integrators/integrators2d.pxd

@@ -25,3 +25,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
+        object _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)
+
+    cdef double _evaluate_orders(self, double x_lower, double x_upper, Function1D y_lower, Function1D y_upper, int x_order, int y_order) except? -1e999
+
+    cdef inline Py_ssize_t _packed_offset(self, int order, int min_order) noexcept