diff --git a/doc/quadrature.rst b/doc/quadrature.rst index 0ed5199..5691cfe 100644 --- a/doc/quadrature.rst +++ b/doc/quadrature.rst @@ -25,6 +25,120 @@ Clenshaw-Curtis and Fejér quadrature in one dimension :members: :show-inheritance: +.. _quadrature-transplanted-1d: + +Transplanted quadrature in one dimension +---------------------------------------- + +The transplanted maps implemented here include the Hale-Trefethen +conformal-map family and the Kosloff-Tal-Ezer map. + +Given a base rule :math:`(s_i, w_i^{(s)})` on :math:`[-1,1]`, transplanted quadrature +uses a map :math:`x=g(s)` to build + +.. math:: + + x_i = g(s_i), \qquad \tilde w_i = w_i^{(s)} g'(s_i), + +so that + +.. math:: + + \int_{-1}^1 f(x)\,dx = \int_{-1}^1 f(g(s))\,g'(s)\,ds + \approx \sum_i \tilde w_i f(x_i). + +Map functions +~~~~~~~~~~~~~ + +.. currentmodule:: modepy.quadrature.transplanted + +Identity map +^^^^^^^^^^^^ + +Use ``map_name="identity"`` for the unmodified base rule. + +.. autofunction:: map_identity + +Sausage polynomial maps +^^^^^^^^^^^^^^^^^^^^^^^ + +Use ``map_name="sausage"`` with odd ``sausage_degree`` (for example +``sausage_degree=5``, ``9``, ``17``) for odd-degree normalized polynomial +truncations of :math:`\arcsin`. + +.. autofunction:: map_sausage + +Kosloff-Tal-Ezer map +^^^^^^^^^^^^^^^^^^^^ + +Use ``map_name="kte"`` (or ``"kosloff_tal_ezer"``). + +* ``kte_rho`` (``>1``) sets the default parameterization + :math:`\alpha = 2/(\rho + \rho^{-1})`. +* ``kte_alpha`` explicitly sets :math:`\alpha` (must satisfy ``0 1``. + +.. note:: + + The strip map requires interior nodes (``abs(s)<1``), so endpoint rules + (for example Gauss-Lobatto or Clenshaw-Curtis) are not valid with + ``map_name="strip"``. + +.. autofunction:: map_strip + +Map dispatcher +^^^^^^^^^^^^^^ + +.. autofunction:: map_trefethen_transplant + +Quadrature wrappers +~~~~~~~~~~~~~~~~~~~ + +.. currentmodule:: modepy + +.. autofunction:: transplanted_1d_quadrature + +.. autofunction:: transplanted_legendre_gauss_quadrature + +Example +~~~~~~~ + +.. code-block:: python + + import modepy as mp + + q_kte = mp.transplanted_legendre_gauss_quadrature( + 20, + map_name="kte", + kte_rho=1.4, + force_dim_axis=True, + ) + + q_sausage = mp.transplanted_legendre_gauss_quadrature( + 20, + map_name="sausage", + sausage_degree=9, + force_dim_axis=True, + ) + +References +~~~~~~~~~~ + +* N. Hale and L. N. Trefethen, *New Quadrature Formulas from Conformal + Maps*, *SIAM Journal on Numerical Analysis* 46(2), 930-948 (2008), + `doi:10.1137/07068607X `__. +* D. Kosloff and H. Tal-Ezer, *A Modified Chebyshev Pseudospectral Method + with an :math:`O(N^{-1})` Time Step Restriction*, + *Journal of Computational Physics* 104(2), 457-469 (1993), + `doi:10.1006/jcph.1993.1044 `__. + Quadratures on the simplex -------------------------- diff --git a/modepy/__init__.py b/modepy/__init__.py index 6f22b14..de1ee5b 100644 --- a/modepy/__init__.py +++ b/modepy/__init__.py @@ -88,6 +88,10 @@ LegendreGaussQuadrature, ) from modepy.quadrature.jaskowiec_sukumar import JaskowiecSukumarQuadrature +from modepy.quadrature.transplanted import ( + transplanted_1d_quadrature, + transplanted_legendre_gauss_quadrature, +) from modepy.quadrature.vioreanu_rokhlin import VioreanuRokhlinSimplexQuadrature from modepy.quadrature.witherden_vincent import WitherdenVincentQuadrature from modepy.quadrature.xiao_gimbutas import XiaoGimbutasSimplexQuadrature @@ -177,6 +181,8 @@ "submesh_for_shape", "symbolicize_function", "tensor_product_nodes", + "transplanted_1d_quadrature", + "transplanted_legendre_gauss_quadrature", "unit_vertices_for_shape", "vandermonde", "warp_and_blend_nodes", diff --git a/modepy/quadrature/transplanted.py b/modepy/quadrature/transplanted.py new file mode 100644 index 0000000..6189e7e --- /dev/null +++ b/modepy/quadrature/transplanted.py @@ -0,0 +1,417 @@ +from __future__ import annotations + + +r""" +.. currentmodule:: modepy.quadrature.transplanted + +Transplanted quadrature applies a smooth map :math:`x=g(s)` to an existing +one-dimensional rule on :math:`[-1,1]`. + +Given base nodes/weights :math:`(s_i, w_i^{(s)})`, the transplanted rule is + +.. math:: + + x_i = g(s_i), \qquad \tilde w_i = w_i^{(s)} g'(s_i), + +so that + +.. math:: + + \int_{-1}^1 f(x)\,dx = \int_{-1}^1 f(g(s)) g'(s)\,ds + \approx \sum_i \tilde w_i f(x_i). + +For map names, parameters, examples, and references, see +:ref:`quadrature-transplanted-1d`. +""" + +from functools import lru_cache +from importlib import import_module +from math import asin, isfinite, log, sqrt +from typing import TYPE_CHECKING, Protocol, cast + +import numpy as np + +from modepy.quadrature import Quadrature + + +if TYPE_CHECKING: + from collections.abc import Callable + + from modepy.typing import ArrayF + + +class _RootScalarResult(Protocol): + converged: bool + root: float + + +class _RootScalarFn(Protocol): + def __call__( + self, + f: Callable[[float], float], + *, + bracket: tuple[float, float], + method: str, + ) -> _RootScalarResult: ... + + +class _EllipkFn(Protocol): + def __call__(self, x: float) -> float: ... + + +class _EllipjFn(Protocol): + def __call__( + self, u: ArrayF, m: float + ) -> tuple[ArrayF, ArrayF, ArrayF, ArrayF]: ... + + +def _scipy_attr(module_name: str, attr_name: str) -> object: + try: + module = import_module(module_name) + except ImportError as exc: + raise RuntimeError( + "The Trefethen strip map requires scipy. " + "Install modepy with scipy support to use map_name='strip'." + ) from exc + + try: + return cast("object", getattr(module, attr_name)) + except AttributeError as exc: + raise RuntimeError( + f"scipy module '{module_name}' is missing required attribute '{attr_name}'" + ) from exc + + +def map_identity(s: ArrayF) -> tuple[ArrayF, ArrayF]: + """Identity transplant map on :math:`[-1, 1]`. + + Returns ``(s, 1)``. + """ + return np.array(s, copy=True), np.ones_like(s) + + +def _arcsin_taylor_coefficients(max_odd_degree: int) -> tuple[float, ...]: + if max_odd_degree < 1 or max_odd_degree % 2 == 0: + raise ValueError( + f"sausage must use positive odd degree, got: {max_odd_degree}" + ) + + nterms = (max_odd_degree + 1) // 2 + coeffs = [1.0] + + for k in range(1, nterms): + km1 = k - 1 + coeffs.append( + coeffs[km1] * (2.0 * km1 + 1.0) ** 2 / (2.0 * k * (2.0 * km1 + 3.0)) + ) + + return tuple(coeffs) + + +def map_sausage(s: ArrayF, degree: int) -> tuple[ArrayF, ArrayF]: + r"""Odd-degree polynomial sausage map from Hale-Trefethen (2008). + + This is the normalized odd Taylor truncation of :math:`\arcsin(s)` + through the monomial of degree *degree*. + + :arg degree: positive odd degree in ``{1, 3, 5, ...}``. + """ + coeffs = _arcsin_taylor_coefficients(degree) + denom = np.asarray(sum(coeffs), dtype=s.dtype) + + g = np.zeros_like(s) + gp = np.zeros_like(s) + + for k, c_k in enumerate(coeffs): + power = 2 * k + 1 + g = g + c_k * s**power + gp = gp + c_k * power * s ** (power - 1) + + return g / denom, gp / denom + + +def map_kosloff_tal_ezer( + s: ArrayF, + *, + rho: float = 1.4, + alpha: float | None = None, +) -> tuple[ArrayF, ArrayF]: + r"""Kosloff-Tal-Ezer map. + + The map is + + .. math:: + + g(s) = \frac{\arcsin(\alpha s)}{\arcsin(\alpha)}, + + where :math:`0 < \alpha < 1`. + + If *alpha* is not provided, it is chosen from *rho* using + + .. math:: + + \alpha = \frac{2}{\rho + \rho^{-1}}, + + matching the parameter choice discussed by Hale-Trefethen for a + :math:`\rho`-ellipse analyticity model. + + Reference: + D. Kosloff and H. Tal-Ezer, "A Modified Chebyshev Pseudospectral + Method with an O(N^{-1}) Time Step Restriction," Journal of + Computational Physics 104(2), 457-469 (1993), + doi:10.1006/jcph.1993.1044. + """ + if alpha is None: + if rho <= 1.0: + raise ValueError(f"KTE parameter rho must be > 1: {rho}") + + alpha = float(2.0 / (rho + 1.0 / rho)) + + if not (0.0 < alpha < 1.0) or not isfinite(alpha): + raise ValueError(f"KTE parameter alpha must satisfy 0 < alpha < 1: {alpha}") + + alpha = float(alpha) + denom = asin(alpha) + + g = np.asarray(np.arcsin(alpha * s) / denom, dtype=np.float64) + gp = np.asarray( + alpha / (denom * np.sqrt(1.0 - alpha**2 * s**2)), + dtype=np.float64, + ) + + return g, gp + + +def _map_preserves_exact_to(map_name: str, *, sausage_degree: int) -> bool: + if map_name == "identity": + return True + + legacy_sausage_degree = _sausage_degree_from_map_name(map_name) + if legacy_sausage_degree is not None: + return legacy_sausage_degree == 1 + + return map_name == "sausage" and sausage_degree == 1 + + +def _sausage_degree_from_map_name(map_name: str) -> int | None: + if not map_name.startswith("sausage_d"): + return None + + degree_text = map_name[len("sausage_d") :] + if not degree_text.isdigit(): + raise ValueError( + f"unsupported sausage map '{map_name}'. Expected format: sausage_d{{odd}}" + ) + + return int(degree_text) + + +@lru_cache(maxsize=16) +def _strip_map_parameter_m(rho: float) -> float: + if rho <= 1.0: + raise ValueError(f"strip map parameter rho must be > 1: {rho}") + + root_scalar = cast("_RootScalarFn", _scipy_attr("scipy.optimize", "root_scalar")) + ellipk = cast("_EllipkFn", _scipy_attr("scipy.special", "ellipk")) + + target = 4.0 * log(rho) / np.pi + + def f(m: float) -> float: + return float(ellipk(1.0 - m) / ellipk(m) - target) + + upper = 1.0 - 1.0e-8 + while f(upper) > 0.0 and 1.0 - upper > 1.0e-16: + upper = 1.0 - (1.0 - upper) / 10.0 + + if f(upper) > 0.0: + raise RuntimeError(f"failed to bracket strip-map parameter for rho={rho}") + + result = root_scalar(f, bracket=(1.0e-14, upper), method="brentq") + if not result.converged: + raise RuntimeError("failed to solve strip-map parameter m") + + return float(result.root) + + +def map_strip(s: ArrayF, *, rho: float = 1.4) -> tuple[ArrayF, ArrayF]: + r"""Strip map from Hale-Trefethen transplanted quadrature. + + :arg rho: strip parameter, must satisfy ``rho > 1``. + + .. important:: + + This map requires interior nodes (``abs(s) < 1``), so it is intended + for base rules such as Legendre-Gauss that do not include endpoints. + """ + if np.any(np.abs(s) >= 1.0): + raise ValueError("strip map expects interior nodes, i.e. abs(s) < 1") + + ellipj = cast("_EllipjFn", _scipy_attr("scipy.special", "ellipj")) + ellipk = cast("_EllipkFn", _scipy_attr("scipy.special", "ellipk")) + + m = _strip_map_parameter_m(rho) + m4 = sqrt(sqrt(m)) + k = float(ellipk(m)) + + omega = 2.0 * k * np.arcsin(s) / np.pi + sn_jacobi, cn_jacobi, dn_jacobi, _ = ellipj(omega, m) + sn = np.asarray(sn_jacobi, dtype=np.float64) + cn = np.asarray(cn_jacobi, dtype=np.float64) + dn = np.asarray(dn_jacobi, dtype=np.float64) + + g = np.asarray(np.arctanh(m4 * sn) / np.arctanh(m4), dtype=np.float64) + gp = np.asarray( + ( + 2.0 + * k + * m4 + * cn + * dn + / ( + np.pi + * np.sqrt(1.0 - s**2) + * (1.0 - np.sqrt(m) * sn**2) + * np.arctanh(m4) + ) + ), + dtype=np.float64, + ) + + return g, gp + + +def map_trefethen_transplant( + s: ArrayF, + map_name: str, + *, + sausage_degree: int = 9, + strip_rho: float = 1.4, + kte_rho: float = 1.4, + kte_alpha: float | None = None, +) -> tuple[ArrayF, ArrayF]: + """Map 1D nodes to a Trefethen transplanted quadrature rule. + + :arg s: nodes on :math:`[-1, 1]`. + :arg map_name: one of ``identity``, ``sausage``, ``kte``, + ``kosloff_tal_ezer``, ``strip``. + :arg sausage_degree: odd polynomial degree for ``map_name="sausage"``. + :arg strip_rho: strip-map parameter for ``map_name="strip"``. + :arg kte_rho: KTE parameter for ``map_name in {"kte", "kosloff_tal_ezer"}`` + when ``kte_alpha`` is not supplied. + :arg kte_alpha: explicit KTE :math:`\alpha` override. + + :returns: ``(mapped_nodes, jacobian)``. + + The supported maps are: + + * ``identity``: :func:`map_identity` + * ``sausage``: :func:`map_sausage` + * ``sausage_d{odd}`` (legacy alias): :func:`map_sausage` + * ``kte`` / ``kosloff_tal_ezer``: :func:`map_kosloff_tal_ezer` + * ``strip``: :func:`map_strip` + + """ + if map_name == "identity": + return map_identity(s) + + if map_name == "sausage": + return map_sausage(s, sausage_degree) + + legacy_sausage_degree = _sausage_degree_from_map_name(map_name) + if legacy_sausage_degree is not None: + return map_sausage(s, legacy_sausage_degree) + + if map_name == "strip": + return map_strip(s, rho=strip_rho) + + if map_name in {"kte", "kosloff_tal_ezer"}: + return map_kosloff_tal_ezer(s, rho=kte_rho, alpha=kte_alpha) + + raise ValueError( + "unsupported map_name " + f"'{map_name}'. Expected one of: " + "identity, sausage, sausage_d{odd}, kte, kosloff_tal_ezer, strip" + ) + + +def transplanted_1d_quadrature( + quadrature: Quadrature, + map_name: str = "sausage", + *, + sausage_degree: int = 9, + strip_rho: float = 1.4, + kte_rho: float = 1.4, + kte_alpha: float | None = None, +) -> Quadrature: + r"""Map an existing 1D quadrature rule using a Trefethen transplant map. + + The transformed rule approximates + + .. math:: + + \int_{-1}^1 f(x)\,dx = \int_{-1}^1 f(g(s)) g'(s)\,ds, + + by mapping existing nodes :math:`s_i` and scaling existing weights :math:`w_i` + with :math:`g'(s_i)`. + """ + base_nodes = quadrature.nodes + if base_nodes.ndim == 1: + nodes_1d = np.asarray(base_nodes) + force_dim_axis = False + elif base_nodes.ndim == 2 and base_nodes.shape[0] == 1: + nodes_1d = cast("ArrayF", base_nodes[0]) + force_dim_axis = True + else: + raise ValueError( + "transplanted_1d_quadrature requires a one-dimensional base quadrature" + ) + + mapped_nodes, jacobian = map_trefethen_transplant( + nodes_1d, + map_name=map_name, + sausage_degree=sausage_degree, + strip_rho=strip_rho, + kte_rho=kte_rho, + kte_alpha=kte_alpha, + ) + mapped_weights = quadrature.weights * jacobian + + if force_dim_axis: + mapped_nodes = np.reshape(mapped_nodes, (1, mapped_nodes.shape[0])) + + exact_to = None + if _map_preserves_exact_to(map_name, sausage_degree=sausage_degree): + try: + exact_to = quadrature.exact_to + except AttributeError: + exact_to = None + + return Quadrature(mapped_nodes, mapped_weights, exact_to=exact_to) + + +def transplanted_legendre_gauss_quadrature( + n: int, + map_name: str = "sausage", + *, + sausage_degree: int = 9, + strip_rho: float = 1.4, + kte_rho: float = 1.4, + kte_alpha: float | None = None, + backend: str | None = None, + force_dim_axis: bool = False, +) -> Quadrature: + r"""Legendre-Gauss quadrature transplanted by a Trefethen map.""" + from modepy.quadrature.jacobi_gauss import LegendreGaussQuadrature + + return transplanted_1d_quadrature( + LegendreGaussQuadrature( + n, + backend=backend, + force_dim_axis=force_dim_axis, + ), + map_name=map_name, + sausage_degree=sausage_degree, + strip_rho=strip_rho, + kte_rho=kte_rho, + kte_alpha=kte_alpha, + ) diff --git a/modepy/test/test_quadrature.py b/modepy/test/test_quadrature.py index 955f6d2..ec7907f 100644 --- a/modepy/test/test_quadrature.py +++ b/modepy/test/test_quadrature.py @@ -25,6 +25,7 @@ import logging +from typing import TYPE_CHECKING, cast import numpy as np import numpy.linalg as la @@ -33,16 +34,21 @@ import modepy as mp +if TYPE_CHECKING: + from modepy.typing import ArrayF + + logger = logging.getLogger(__name__) def test_transformed_quadrature(): """Test 1D quadrature on arbitrary intervals""" - def gaussian_density(x, mu, sigma): - return ( - 1 / (sigma * np.sqrt(2*np.pi)) - * np.exp(-np.sum((x-mu)**2, axis=0) / (2 * sigma**2)) + def gaussian_density(x: ArrayF, mu: float, sigma: float) -> ArrayF: + sq_dist = cast("ArrayF", np.sum((x - mu) ** 2, axis=0)) + return cast( + "ArrayF", + 1 / (sigma * np.sqrt(2 * np.pi)) * np.exp(-sq_dist / (2 * sigma**2)), ) from modepy.quadrature import Transformed1DQuadrature @@ -51,8 +57,10 @@ def gaussian_density(x, mu, sigma): mu = 17 sigma = 12 tq = Transformed1DQuadrature( - LegendreGaussQuadrature(20, force_dim_axis=True), - left=mu - 6*sigma, right=mu + 6*sigma) + LegendreGaussQuadrature(20, force_dim_axis=True), + left=mu - 6 * sigma, + right=mu + 6 * sigma, + ) result = tq(lambda x: gaussian_density(x, mu, sigma)) assert abs(result - 1) < 1.0e-9 @@ -67,10 +75,13 @@ def gaussian_density(x, mu, sigma): @pytest.mark.parametrize("backend", BACKENDS) -@pytest.mark.parametrize("quad_type", [ - mp.LegendreGaussQuadrature, - mp.LegendreGaussLobattoQuadrature, - ]) +@pytest.mark.parametrize( + "quad_type", + [ + mp.LegendreGaussQuadrature, + mp.LegendreGaussLobattoQuadrature, + ], +) def test_gauss_quadrature(backend, quad_type): for s in range(9 + 1): if quad_type == mp.LegendreGaussLobattoQuadrature and s == 0: @@ -79,13 +90,14 @@ def test_gauss_quadrature(backend, quad_type): quad = quad_type(s, backend=backend, force_dim_axis=True) - assert quad.nodes.shape[1] == s+1 + assert quad.nodes.shape[1] == s + 1 for deg in range(quad.exact_to + 1): + def f(x): return np.sum(x**deg, axis=0) # noqa: B023 i_f = quad(f) - i_f_true = 1 / (deg+1) * (1 - (-1)**(deg + 1)) + i_f_true = 1 / (deg + 1) * (1 - (-1) ** (deg + 1)) err = abs(i_f - i_f_true) assert err < 3.0e-15, (s, deg, err, i_f, i_f_true) @@ -95,17 +107,152 @@ def test_clenshaw_curtis_quadrature() -> None: for s in range(1, 9 + 1): quad = ClenshawCurtisQuadrature(s, force_dim_axis=True) - assert quad.nodes.shape[1] == s+1 + assert quad.nodes.shape[1] == s + 1 for deg in range(quad.exact_to + 1): + def f(x): return x**deg # noqa: B023 i_f = quad(f) - i_f_true = 1 / (deg+1) * (1 - (-1)**(deg + 1)) + i_f_true = 1 / (deg + 1) * (1 - (-1) ** (deg + 1)) err = abs(i_f - i_f_true) assert err < 2.0e-15, (s, deg, err, i_f, i_f_true) +@pytest.mark.parametrize( + ("map_name", "sausage_degree"), + [ + ("identity", 9), + ("sausage", 5), + ("sausage", 9), + ("sausage", 17), + ("kte", 9), + ("kosloff_tal_ezer", 9), + ], +) +def test_transplanted_legendre_gauss_quadrature( + map_name: str, + sausage_degree: int, +) -> None: + base = mp.LegendreGaussQuadrature(15, force_dim_axis=True) + transplanted = mp.transplanted_legendre_gauss_quadrature( + 15, + map_name=map_name, + sausage_degree=sausage_degree, + force_dim_axis=True, + ) + + base_nodes = cast("ArrayF", np.asarray(base.nodes[0], dtype=np.float64)) + trans_nodes = cast("ArrayF", np.asarray(transplanted.nodes[0], dtype=np.float64)) + base_weights = cast("ArrayF", np.asarray(base.weights, dtype=np.float64)) + trans_weights = cast("ArrayF", np.asarray(transplanted.weights, dtype=np.float64)) + + assert transplanted.nodes.shape == base.nodes.shape + assert transplanted.weights.shape == base.weights.shape + + from modepy.quadrature.transplanted import map_trefethen_transplant + + mapped_nodes, mapped_jacobian = map_trefethen_transplant( + base_nodes, + map_name=map_name, + sausage_degree=sausage_degree, + ) + + assert la.norm(trans_nodes - mapped_nodes, np.inf) < 1.0e-14 + assert la.norm(trans_weights - base_weights * mapped_jacobian, np.inf) < 1.0e-14 + + # Sausage maps are polynomial maps, so integrating constants + # should remain exact. + if map_name == "sausage": + err = abs(float(np.sum(trans_weights)) - 2.0) + assert err < 1.0e-14 + + +def test_transplanted_strip_map_quadrature() -> None: + pytest.importorskip("scipy") + + transplanted = mp.transplanted_legendre_gauss_quadrature( + 32, + map_name="strip", + strip_rho=1.4, + force_dim_axis=True, + ) + + strip_nodes = cast("ArrayF", np.asarray(transplanted.nodes[0], dtype=np.float64)) + strip_weights = cast("ArrayF", np.asarray(transplanted.weights, dtype=np.float64)) + + assert np.all(np.diff(strip_nodes) > 0) + assert np.all(strip_weights > 0) + assert abs(float(np.sum(strip_weights)) - 2.0) < 1.0e-9 + + +def test_transplanted_strip_map_rejects_endpoint_rules() -> None: + pytest.importorskip("scipy") + + with pytest.raises(ValueError, match="interior nodes"): + mp.transplanted_1d_quadrature( + mp.ClenshawCurtisQuadrature(5, force_dim_axis=True), map_name="strip" + ) + + +def test_transplanted_sausage_map_rejects_even_degrees() -> None: + with pytest.raises(ValueError, match="positive odd degree"): + mp.transplanted_legendre_gauss_quadrature( + 8, + map_name="sausage", + sausage_degree=4, + force_dim_axis=True, + ) + + +def test_transplanted_sausage_map_name_matches_direct_map() -> None: + from modepy.quadrature.transplanted import map_sausage, map_trefethen_transplant + + s = np.linspace(-0.95, 0.95, 33) + + for degree in [5, 9, 17]: + g, gp = map_sausage(s, degree=degree) + g_ref, gp_ref = map_trefethen_transplant( + s, + map_name="sausage", + sausage_degree=degree, + ) + assert la.norm(g - g_ref, np.inf) < 1.0e-15 + assert la.norm(gp - gp_ref, np.inf) < 1.0e-15 + + +def test_transplanted_kte_map_name_matches_direct_map() -> None: + from modepy.quadrature.transplanted import ( + map_kosloff_tal_ezer, + map_trefethen_transplant, + ) + + s = np.linspace(-1.0, 1.0, 33) + rho = 1.4 + alpha = 2.0 / (rho + 1.0 / rho) + + g, gp = map_kosloff_tal_ezer(s, rho=rho) + g_ref, gp_ref = map_kosloff_tal_ezer(s, alpha=alpha) + assert la.norm(g - g_ref, np.inf) < 1.0e-15 + assert la.norm(gp - gp_ref, np.inf) < 1.0e-15 + + g_name, gp_name = map_trefethen_transplant(s, map_name="kte", kte_rho=rho) + assert la.norm(g - g_name, np.inf) < 1.0e-15 + assert la.norm(gp - gp_name, np.inf) < 1.0e-15 + + +def test_transplanted_kte_map_rejects_invalid_parameters() -> None: + from modepy.quadrature.transplanted import map_trefethen_transplant + + s = np.array([0.0]) + + with pytest.raises(ValueError, match="rho must be > 1"): + map_trefethen_transplant(s, map_name="kte", kte_rho=1.0) + + with pytest.raises(ValueError, match="0 < alpha < 1"): + map_trefethen_transplant(s, map_name="kte", kte_alpha=1.0) + + @pytest.mark.parametrize("kind", [1, 2]) def test_fejer_quadrature(kind: int) -> None: from modepy.quadrature.clenshaw_curtis import FejerQuadrature @@ -118,17 +265,20 @@ def f(x): return x**deg # noqa: B023 i_f = quad(f) - i_f_true = 1 / (deg+1) * (1 - (-1)**(deg + 1)) + i_f_true = 1 / (deg + 1) * (1 - (-1) ** (deg + 1)) err = abs(i_f - i_f_true) assert err < 2.0e-15, (s, deg, err, i_f, i_f_true) -@pytest.mark.parametrize(("quad_class", "highest_order"), [ - (mp.XiaoGimbutasSimplexQuadrature, None), - (mp.JaskowiecSukumarQuadrature, None), - (mp.VioreanuRokhlinSimplexQuadrature, None), - (mp.GrundmannMoellerSimplexQuadrature, 3), - ]) +@pytest.mark.parametrize( + ("quad_class", "highest_order"), + [ + (mp.XiaoGimbutasSimplexQuadrature, None), + (mp.JaskowiecSukumarQuadrature, None), + (mp.VioreanuRokhlinSimplexQuadrature, None), + (mp.GrundmannMoellerSimplexQuadrature, 3), + ], +) @pytest.mark.parametrize("dim", [2, 3]) def test_simplex_quadrature(quad_class, highest_order, dim) -> None: """Check that quadratures on simplices works as advertised""" @@ -154,6 +304,7 @@ def test_simplex_quadrature(quad_class, highest_order, dim) -> None: if 0: import matplotlib.pyplot as pt + pt.plot(quad.nodes[0], quad.nodes[1]) pt.show() @@ -172,10 +323,13 @@ def test_simplex_quadrature(quad_class, highest_order, dim) -> None: @pytest.mark.parametrize("dim", [2, 3]) -@pytest.mark.parametrize(("quad_class", "max_order"), [ - (mp.WitherdenVincentQuadrature, np.inf), - (mp.LegendreGaussTensorProductQuadrature, 6), - ]) +@pytest.mark.parametrize( + ("quad_class", "max_order"), + [ + (mp.WitherdenVincentQuadrature, np.inf), + (mp.LegendreGaussTensorProductQuadrature, 6), + ], +) def test_hypercube_quadrature(dim, quad_class, max_order): from pytools import ( generate_nonnegative_integer_tuples_summing_to_at_most as gnitstam, @@ -183,14 +337,19 @@ def test_hypercube_quadrature(dim, quad_class, max_order): from modepy.tools import Monomial - def _check_monomial(quad, comb): + def _check_monomial(quad: mp.Quadrature, comb: tuple[int, ...]) -> float: f = Monomial(comb) - int_approx = quad(f) - int_exact = 2**dim * f.hypercube_integral() + int_approx_obj = quad(f) + if isinstance(int_approx_obj, np.ndarray): + int_approx = float(int_approx_obj.item()) + else: + int_approx = float(int_approx_obj) + int_exact = cast("float", 2**dim * f.hypercube_integral()) error = abs(int_approx - int_exact) / abs(int_exact) - logger.info("%s: %.5e %.5e / rel error %.5e", - comb, int_approx, int_exact, error) + logger.info( + "%s: %.5e %.5e / rel error %.5e", comb, int_approx, int_exact, error + ) return error @@ -202,14 +361,15 @@ def _check_monomial(quad, comb): logger.info("UNAVAILABLE at order %d", order) break + quad_exact_to = cast("int", quad.exact_to) + assert np.all(quad.weights > 0) - logger.info("quadrature: %s %d %d", - quad_class.__name__.lower(), order, quad.exact_to) - for comb in gnitstam(quad.exact_to, dim): + logger.info("quadrature: order %d %d", order, quad_exact_to) + for comb in gnitstam(quad_exact_to, dim): assert _check_monomial(quad, comb) < 5.0e-15 - comb = (0,) * (dim - 1) + (quad.exact_to + 1,) + comb = (0,) * (dim - 1) + (quad_exact_to + 1,) assert _check_monomial(quad, comb) > 5.0e-15 order += 2 @@ -229,6 +389,7 @@ def test_mass_quadrature_is_newton_cotes(shape: mp.Shape, order: int) -> None: assert basis.orthonormality_weight() == 1 from math import factorial + shape_volume = 2**shape.dim / factorial(shape.dim) # integrals are orthogonal to the constant @@ -237,9 +398,13 @@ def test_mass_quadrature_is_newton_cotes(shape: mp.Shape, order: int) -> None: integrals[0] = shape_volume * basis.functions[0](np.zeros(shape.dim)) newton_cotes_weights = la.solve(vdm.T, integrals) + mass_weights = cast("ArrayF", mass_weights) + newton_cotes_weights = cast("ArrayF", newton_cotes_weights) - assert (la.norm(newton_cotes_weights - mass_weights, np.inf) - / la.norm(newton_cotes_weights, np.inf)) < 1e-13 + assert ( + la.norm(newton_cotes_weights - mass_weights, np.inf) + / la.norm(newton_cotes_weights, np.inf) + ) < 1e-13 # You can test individual routines by typing @@ -247,10 +412,12 @@ def test_mass_quadrature_is_newton_cotes(shape: mp.Shape, order: int) -> None: if __name__ == "__main__": import sys + if len(sys.argv) > 1: exec(sys.argv[1]) else: from pytest import main + main([__file__]) # vim: fdm=marker