From f9581b5a0888ad453a44c42430bea63216a1cbeb Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Sat, 18 Apr 2026 16:24:08 +0300 Subject: [PATCH 01/10] feat: implement LU-based kernel rewriting --- sumpy/kernel_rewrite.py | 386 ++++++++++++++++++++++++++++++ sumpy/test/test_kernel_rewrite.py | 183 ++++++++++++++ 2 files changed, 569 insertions(+) create mode 100644 sumpy/kernel_rewrite.py create mode 100644 sumpy/test/test_kernel_rewrite.py diff --git a/sumpy/kernel_rewrite.py b/sumpy/kernel_rewrite.py new file mode 100644 index 00000000..5ee7841f --- /dev/null +++ b/sumpy/kernel_rewrite.py @@ -0,0 +1,386 @@ +from __future__ import annotations + + +__copyright__ = """ +Copyright (C) 2012 Andreas Kloeckner +Copyright (C) 2020 Isuru Fernando +Copyright (C) 2026 Alexandru Fikl +""" + +__license__ = """ +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in +all copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN +THE SOFTWARE. +""" + +import logging +from dataclasses import dataclass +from typing import TYPE_CHECKING, Any, NamedTuple + +import numpy as np + +from pytools import ( + generate_nonnegative_integer_tuples_summing_to_at_most as gnitstam, +) + +import sumpy.symbolic as sym + + +if TYPE_CHECKING: + from collections.abc import Sequence + + import optype.numpy as onp + + from pymbolic.typing import ArithmeticExpression + + from sumpy.expansions.diff_op import MultiIndex + from sumpy.kernel import ScalarKernel + +logger = logging.getLogger(__name__) + + +# {{{ rewrite_using_base_kernel + +@dataclass(frozen=True) +class LinearOperatorRepresentation: + r"""Expresses a target kernel as a linear operator acting on a base kernel. + + .. math:: + + G(\boldsymbol{r}) = C + \sum_{|\boldsymbol{\alpha}| < p} + c_{\boldsymbol{\alpha}} + \frac{\partial^{|\boldsymbol{\alpha}|} G_0} + {\partial \boldsymbol{r}^{\boldsymbol{\alpha}}} + + .. autoattribute:: target_kernel + .. autoattribute:: base_kernel + .. autoattribute:: mis + .. autoattribute:: coeffs + """ + + target_kernel: ScalarKernel + """The target kernel expressed in terms of :attr:`base_kernel`.""" + base_kernel: ScalarKernel + """A base kernel used to express the target kernel.""" + + mis: Sequence[MultiIndex] + """Multi-indices for each non-zero term in the linear combination.""" + coeffs: Sequence[ArithmeticExpression] + """Constant coefficients in the linear combination. Note that the coefficients + have length ``len(mis) + 1``, where the first element is always the constant term. + """ + + def pretty(self) -> str: + from sumpy.kernel import AxisTargetDerivative + + terms = [] + if self.coeffs[0] != 0: + terms.append(str(self.coeffs[0])) + + for mi, c in zip(self.mis, self.coeffs[1:], strict=True): + expr = self.base_kernel + for d, n in enumerate(mi): + for _ in range(n): + expr = AxisTargetDerivative(d, expr) + + terms.append(f"{c} * {expr}") + + return f"{self.target_kernel} = " + " + ".join(terms) + + +def rewrite_using_base_kernel( + target_kernel: ScalarKernel, + base_kernel: ScalarKernel, + *, + order: int | None = None, + atol: float = 1.0e-10, + ) -> LinearOperatorRepresentation: + pde = base_kernel.get_pde_as_diff_op() + if order is None: + order = pde.order + + # TODO: pick the best algorithm here + return rewrite_using_base_kernel_lu( + target_kernel, base_kernel, + min_order=order, + ) + +# }}} + + +# {{{ rewrite_using_base_kernel_lu + +INT_MAX = 10 ** 15 + + +class FactorizationFailedError(Exception): + pass + + +class RewriteFailedError(Exception): + pass + + +class _LUDecomposition(NamedTuple): + L: sym.Matrix + U: sym.Matrix + permutation: Sequence[tuple[int, int]] + + mis: Sequence[MultiIndex] + """The multi-indices for which the derivatives were computed. These correspond + to rows in the matrix and should be used to recover the expansion of the + target kernel in terms of the base kernel. + """ + points: onp.Array2D[Any] + """An array of shape ``(dim, npoints)`` of points where the base kernel was + evaluated to compute the current LU factorization. + """ + + +def evalf(expr: sym.Expr, prec: int = 100) -> sym.Expr: + """Evaluate an expression numerically using ``prec`` number of bits.""" + from sumpy.symbolic import USE_SYMENGINE + + if USE_SYMENGINE: + return expr.n(prec=prec) + else: + import sympy + dps = int(sympy.log(2**prec, 10)) + return expr.n(n=dps) + + +def round_expr( + expr: sym.Basic, atol: float = 1.0e-8, rtol: float = 1.0e-5 + ) -> sym.Basic: + """Round all numeric values in *expr* to the nearest integer. + + This function clips all numbers close to zero (effectively removing them + from the expression due to SymPy's automatic symplifications) and rounds all + numbers to the nearest integer using the given *atol* and *rtol*. + """ + nums = expr.atoms(sym.Float) + replace_dict: dict[Any, float] = {} + + for num in nums: + value = float(num) + nearest_int = round(value) + + if abs(value - nearest_int) < atol + rtol * abs(value): + replace_dict[num] = sym.Integer(nearest_int) + else: + replace_dict[num] = value + + return expr.xreplace(replace_dict) + + +def simplify(expr: sym.Basic) -> sym.Basic: + if sym.USE_SYMENGINE: + return expr.simplify() + else: + import sympy as sp + + return sp.simplify(expr) + + +def rewrite_using_base_kernel_lu( + target_kernel: ScalarKernel, + base_kernel: ScalarKernel, + *, + min_order: int | None = None, + retries: int = 5, + rng: np.random.Generator | None = None, + ) -> LinearOperatorRepresentation: + """Find a relation between the *target_kernel* and the *base_kernel* using + a numerical LU-based algorithm. + + The algorithm samples the *base_kernel* and its derivatives at random + points to get a matrix ``A``. It also samples the target kernel at the same + points to get a vector ``b`` and solving for the system ``A c = b`` using + an LU factorization of ``A``. The solution ``c`` is the vector of coefficients + in the linear combination :class:`LinearOperatorRepresentation`. + + :arg min_order: starting minimum derivative order to use when attempting the + decomposition. By default, this will be the order of the PDE solved by + *base_kernel*. + :arg retries: maximum number of retries for each order. If the LU decomposition + fails due to a poor choice of random points, it is retried several times. + """ + + pde = base_kernel.get_pde_as_diff_op() + if min_order is None: + min_order = pde.order + + if min_order > pde.order: + raise NotImplementedError( + "Rewriting when the base kernel's derivatives are linearly dependent " + "is not implemented") + + if rng is None: + rng = np.random.default_rng() + + coeffs: list[sym.Basic] = [] + mis: list[MultiIndex] = [] + + dim = base_kernel.dim + dvec = sym.make_sym_vector("d", dim) + target_expr = target_kernel.get_expression(dvec) + + target_scaling = target_kernel.get_global_scaling_const() + base_scaling = base_kernel.get_global_scaling_const() + + order = min_order + to_pymbolic = sym.SympyToPymbolicMapperWithSymbols() + while order <= pde.order: + try: + lu = _get_base_kernel_matrix_lu_factorization( + base_kernel, order, rng=rng, retries=retries + ) + except FactorizationFailedError as exc: + if order == pde.order: + raise RewriteFailedError( + f"failed to compute LU factorization for orders in " + f"[{min_order}, {pde.order}] for base kernel {base_kernel} " + f"after {retries} retries" + ) from exc + + order += 1 + continue + + # evaluate right-hand side + b = sym.Matrix([ + target_expr.xreplace(dict(zip(dvec, lu.points[:, i], strict=True))) + for i in range(lu.points.shape[1]) + ]) + + # solve + all_coeffs = sym.solve_lu(lu.L, lu.U, lu.permutation, b, + postprocess=lambda x: x.expand()) + + # gather all non-zero coefficients from the result + const = 0 + coeffs = [] + mis = [] + for i, coeff in enumerate(all_coeffs): + coeff = round_expr(evalf(coeff)) + if coeff == 0: + continue + + if i == 0: + const = to_pymbolic(simplify(coeff * target_scaling)) + logger.debug(" %s", coeff) + else: + mi = lu.mis[i - 1] + coeff = simplify(coeff * target_scaling / base_scaling) + + mis.append(mi) + coeffs.append(to_pymbolic(coeff)) + logger.debug(" + %s*%s.diff%s", coeff, base_kernel, mi) + + if coeffs: + coeffs.insert(0, const) + break + + order += 1 + + if not coeffs: + raise RewriteFailedError( + f"could not express {target_kernel} in terms of {base_kernel}" + ) + + return LinearOperatorRepresentation(target_kernel, base_kernel, mis, coeffs) + + +def _get_base_kernel_matrix_lu_factorization( + base_kernel: ScalarKernel, + order: int, + *, + rng: np.random.Generator, + retries: int, + ) -> _LUDecomposition: + pde = base_kernel.get_pde_as_diff_op() + if order > pde.order: + raise NotImplementedError( + "Rewriting when the base kernel's derivatives are linearly dependent " + "is not implemented") + + dim = base_kernel.dim + + mis = list(gnitstam(order, dim)) + if order == pde.order: + pde_mis = [ident.mi for eq in pde.eqs for ident in eq] + pde_mis = [mi for mi in pde_mis if sum(mi) == order] + mis.remove(pde_mis[-1]) + + logger.debug("Removing %s to avoid linear dependent mis", pde_mis[-1]) + + # get sympy expression for the base kernel + dvec = sym.make_sym_vector("d", dim) + base_expr = base_kernel.get_expression(dvec) + + # evaluate all the needed derivatives + mi_to_derivative: dict[MultiIndex, sym.Basic] = {} + for mi in mis: + expr = base_expr + for i, nderivs in enumerate(mi): + if nderivs == 0: + continue + expr = expr.diff(dvec[i], nderivs) + + mi_to_derivative[mi] = expr + + # try to LU factorize on random points + for _ in range(retries): + # TODO: is it faster to generate numbers and then sympify them? + points = np.empty((dim, len(mis) + 1), dtype=object) + for i in range(points.shape[0]): + for j in range(points.shape[1]): + points[i, j] = sym.Integer(rng.integers(1, INT_MAX)) / INT_MAX + + # evaluate derivatives at points and construct matrix + entries: list[list[sym.Basic]] = [] + for i in range(points.shape[1]): + row: list[sym.Basic] = [sym.Integer(1)] + + for mi in mis: + expr = mi_to_derivative[mi].replace( + dict(zip(dvec, points[:, i], strict=True)) + ) + row.append(evalf(expr)) + + entries.append(row) + mat = sym.Matrix(entries) + + # TODO: LUdecomposition in symengine is not implemented for non-square matrices + try: + L, U, perm = mat.LUdecomposition() # ruff: ignore[non-lowercase-variable-in-function] + except RuntimeError: + # NOTE: symengine seems to throw a SymEngineError -> RuntimeError when + # it fails to do the LU factorization due to rank-deficiency + continue + else: + # NOTE: and sympy seems to set the last row of U to 0 + if not sym.USE_SYMENGINE and all(expr == 0 for expr in U[-1, :]): + continue + + return _LUDecomposition(L, U, perm, mis, points) + + raise FactorizationFailedError( + f"failed to compute LU factorization to order {order} for {base_kernel} " + f"after {retries} retries" + ) + +# }}} diff --git a/sumpy/test/test_kernel_rewrite.py b/sumpy/test/test_kernel_rewrite.py new file mode 100644 index 00000000..9d3c42c1 --- /dev/null +++ b/sumpy/test/test_kernel_rewrite.py @@ -0,0 +1,183 @@ +from __future__ import annotations + + +__copyright__ = """ +Copyright (C) 2020 Isuru Fernando +Copyright (C) 2026 Alexandru Fikl +""" + +import logging + +import numpy as np +import pytest + +import sumpy.symbolic as sym +from sumpy.kernel import ( + BiharmonicKernel, + LaplaceKernel, + StokesletComponentKernel, + StressletComponentKernel, +) +from sumpy.kernel_rewrite import ( + LinearOperatorRepresentation, + rewrite_using_base_kernel_fourier, + rewrite_using_base_kernel_lu, +) + + +logger = logging.getLogger(__name__) + + +def mi_derivative(expr: sym.Expr, x: sym.Matrix, mi: tuple[int, ...]) -> sym.Expr: + result = expr + for i, n in enumerate(mi): + result = result.diff(x[i], n) + + return result + + +def check_kernel_rewrite(op: LinearOperatorRepresentation) -> None: + from sumpy.kernel_rewrite import evalf, simplify + from sumpy.symbolic import PymbolicToSympyMapperWithSymbols + + dim = op.target_kernel.dim + dvec = sym.make_sym_vector("d", dim) + to_sympy = PymbolicToSympyMapperWithSymbols() + + target_expr = ( + op.target_kernel.get_global_scaling_const() + * op.target_kernel.get_expression(dvec)) + base_expr = ( + op.base_kernel.get_global_scaling_const() + * op.base_kernel.get_expression(dvec)) + + expr = to_sympy(op.coeffs[0]) + sum(( + to_sympy(c) * mi_derivative(base_expr, dvec, mi) + for c, mi in zip(op.coeffs[1:], op.mis, strict=True) + ), sym.Integer(0)) + + result = evalf(simplify(target_expr - expr)) + assert abs(result) < 3.0 * 1.0e-16 + + +# {{{ test_rewrite_using_base_kernel_lu_laplace_biharmonic + +@pytest.mark.parametrize("dim", [2, 3]) +def test_rewrite_using_base_kernel_lu_laplace_biharmonic(dim: int) -> None: + rng = np.random.default_rng(seed=42) + + base_kernel = BiharmonicKernel(dim) + target_kernel = LaplaceKernel(dim) + result = rewrite_using_base_kernel_lu(target_kernel, base_kernel, rng=rng) + + print(result.pretty()) + check_kernel_rewrite(result) + +# }}} + + +# {{{ test_rewrite_using_base_kernel_lu_stokeslet_biharmonic + +@pytest.mark.parametrize("dim", [2, 3]) +def test_rewrite_using_base_kernel_lu_stokeslet_biharmonic(dim: int) -> None: + from pytools import generate_nonnegative_integer_tuples_below as gnitb + + rng = np.random.default_rng(seed=42) + + base_kernel = BiharmonicKernel(dim) + for i, j in gnitb(dim, 2): + target_kernel = StokesletComponentKernel(dim, i, j, viscosity_mu_name="mu") + result = rewrite_using_base_kernel_lu(target_kernel, base_kernel, rng=rng) + print(result.pretty()) + check_kernel_rewrite(result) + + +# }}} + + +# {{{ test_rewrite_using_base_kernel_lu_stresslet_biharmonic + +@pytest.mark.parametrize("dim", [2, 3]) +def test_rewrite_using_base_kernel_lu_stresslet_biharmonic(dim: int) -> None: + from pytools import generate_nonnegative_integer_tuples_below as gnitb + + rng = np.random.default_rng(seed=42) + + base_kernel = BiharmonicKernel(dim) + for i, j, k in gnitb(dim, 3): + target_kernel = StressletComponentKernel(dim, i, j, k, viscosity_mu_name="mu") + result = rewrite_using_base_kernel_lu(target_kernel, base_kernel, rng=rng) + print(result.pretty()) + check_kernel_rewrite(result) + + +# }}} + + +# {{{ test_rewrite_using_base_kernel_fourier_laplace_biharmonic + + +@pytest.mark.parametrize("dim", [2, 3]) +def test_rewrite_using_base_kernel_fourier_laplace_biharmonic(dim: int) -> None: + """Test that the Fourier-based algorithm recovers Laplace from biharmonic.""" + base_kernel = BiharmonicKernel(dim) + target_kernel = LaplaceKernel(dim) + result = rewrite_using_base_kernel_fourier(target_kernel, base_kernel) + + logger.info(result.pretty()) + check_kernel_rewrite(result) + + +# }}} + + +# {{{ test_rewrite_using_base_kernel_fourier_stokeslet_biharmonic + + +@pytest.mark.parametrize("dim", [2, 3]) +def test_rewrite_using_base_kernel_fourier_stokeslet_biharmonic(dim: int) -> None: + """Test that the Fourier-based algorithm recovers the Stokeslet from biharmonic.""" + from itertools import product + + base_kernel = BiharmonicKernel(dim) + + for i, j in product(range(dim), repeat=2): + target_kernel = StokesletComponentKernel(dim, i, j, viscosity_mu_name="mu") + result = rewrite_using_base_kernel_fourier(target_kernel, base_kernel) + + logger.info(result.pretty()) + check_kernel_rewrite(result) + +# }}} + + +# {{{ test_rewrite_using_base_kernel_fourier_indivisible + + +@pytest.mark.parametrize("dim", [2, 3]) +def test_rewrite_using_base_kernel_fourier_biharmonic_laplace(dim: int) -> None: + from sumpy.kernel_rewrite import ( + RewriteFailedError, + rewrite_using_base_kernel_fourier, + ) + + # Laplace Fourier symbol (-|k|^2) is not divisible by biharmonic (|k|^4) + base_kernel = LaplaceKernel(dim) + target_kernel = BiharmonicKernel(dim) + + with pytest.raises(RewriteFailedError, match="cannot rewrite"): + rewrite_using_base_kernel_fourier(target_kernel, base_kernel) + + +# }}} + + +if __name__ == "__main__": + import sys + + if len(sys.argv) > 1: + exec(sys.argv[1]) + else: + pytest.main([__file__]) + +# vim: fdm=marker From 02b8092f246c32494e17e094e30007eefa8dd4b1 Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Fri, 28 Aug 2026 16:41:52 +0300 Subject: [PATCH 02/10] feat: add Fourier based kernel rewrite --- sumpy/kernel_rewrite.py | 140 ++++++++++++++++++++++++++++++++++++++-- 1 file changed, 135 insertions(+), 5 deletions(-) diff --git a/sumpy/kernel_rewrite.py b/sumpy/kernel_rewrite.py index 5ee7841f..339595d7 100644 --- a/sumpy/kernel_rewrite.py +++ b/sumpy/kernel_rewrite.py @@ -55,6 +55,11 @@ # {{{ rewrite_using_base_kernel + +class RewriteFailedError(Exception): + pass + + @dataclass(frozen=True) class LinearOperatorRepresentation: r"""Expresses a target kernel as a linear operator acting on a base kernel. @@ -119,19 +124,144 @@ def rewrite_using_base_kernel( min_order=order, ) + # }}} -# {{{ rewrite_using_base_kernel_lu +# {{{ rewrite_using_base_kernel_fourier -INT_MAX = 10 ** 15 +def rewrite_using_base_kernel_fourier( + target_kernel: ScalarKernel, + base_kernel: ScalarKernel, +) -> LinearOperatorRepresentation: + r"""Find a relation between *target_kernel* and *base_kernel* using the + Fourier symbol of their respective PDEs. -class FactorizationFailedError(Exception): - pass + The algorithm works by computing the Fourier symbols :math:`P_{\text{base}}` + and :math:`P_{\text{target}}` of the scalar PDEs satisfied by the two + kernels. The target kernel can be expressed as a differential operator + applied to the base kernel if and only if the ratio + + .. math:: + D(i \boldsymbol{k}) = \frac{P_{\text{base}}(\boldsymbol{k})} + {P_{\text{target}}(\boldsymbol{k})} -class RewriteFailedError(Exception): + is a polynomial in :math:`\boldsymbol{k}`. When it is, each monomial + :math:`\prod_j (i k_j)^{\alpha_j}` corresponds to the derivative + :math:`\partial^{|\alpha|} / \partial \boldsymbol{r}^{\alpha}`. + """ + import sympy as sp + + # FIXME: want to add a general check for this? + try: + _ = base_kernel.get_pde_system_kernel() + except TypeError: + pass + else: + raise ValueError( + f"'base_kernel' cannot be part of a system: {type(base_kernel)}" + ) + + dim = base_kernel.dim + if target_kernel.dim != dim: + raise ValueError( + f"kernel dimension mismatch: {target_kernel.dim} (target_kernel) and " + f"{dim} (base_kernel)" + ) + + from sumpy.expansion.diff_op import to_fourier_matrix + + ks = sp.Matrix([sp.Symbol(f"_k{j}") for j in range(dim)]) + + pde_base = base_kernel.get_pde_as_diff_op() + fourier_base = to_fourier_matrix(pde_base, ks).inv() + + try: + target_system_kernel, idx = target_kernel.get_pde_system_kernel() + except TypeError: + target_system_kernel, idx = None, None + + if target_system_kernel is None: + pde_target = target_kernel.get_pde_as_diff_op() + fourier_target = to_fourier_matrix(pde_target, ks).inv() + else: + assert idx is not None + + pde_target_system = target_system_kernel.get_pde_as_diff_op() + fourier_target_system = to_fourier_matrix(pde_target_system, ks) + fourier_target_system_inv = sp.simplify(fourier_target_system.inv()) + + fourier_target = sp.Matrix([fourier_target_system_inv[idx]]) + + p_base: sp.Expr = sp.Integer(1) + for entry in fourier_base: + p_base = p_base * entry + + p_target: sp.Expr = sp.Integer(1) + for entry in fourier_target: + p_target = p_target * entry + + quotient = sp.simplify(p_target / p_base) + if not quotient.is_polynomial(*ks): + raise RewriteFailedError( + f"cannot rewrite {target_kernel} in terms of {base_kernel}" + ) + + mis: list[MultiIndex] = [] + coeffs: list[sp.Expr] = [] + + for exponent, coeff in sp.Poly(quotient, *ks).as_dict().items(): + coeff = sp.simplify(coeff) + if coeff == sp.Integer(0): + continue + + mis.append(exponent) + coeffs.append(coeff) + + if not mis: + raise RewriteFailedError( + f"cannot rewrite {target_kernel} in terms of {base_kernel}" + ) + + # compute constant term, if any + dvec = sym.make_sym_vector("d", dim) + base_scaled = ( + base_kernel.get_global_scaling_const() + * base_kernel.get_expression(dvec) + ) + + # build the differential operator part in real space + const: sp.Expr = ( + target_kernel.get_global_scaling_const() + * target_kernel.get_expression(dvec) + ) + + for mi, c in zip(mis, coeffs, strict=True): + term = c * base_scaled + for i, n in enumerate(mi): + term = term.diff(dvec[i], n) + const = const - sp.simplify(term) + const = sp.simplify(const) + + return LinearOperatorRepresentation( + target_kernel, + base_kernel, + mis, + [sym.to_pymbolic(expr) for expr in [const, *coeffs]], + ) + + +# }}} + + +# {{{ rewrite_using_base_kernel_lu + +INT_MAX = 10**15 + + +class FactorizationFailedError(Exception): pass From ce6644b593a4ceb0ca3d17a4bd806bcb56f6e54f Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Fri, 11 Sep 2026 14:31:51 +0300 Subject: [PATCH 03/10] feat: improve point generation --- sumpy/kernel_rewrite.py | 47 +++++++++++++++++++++++++++++++++++------ 1 file changed, 40 insertions(+), 7 deletions(-) diff --git a/sumpy/kernel_rewrite.py b/sumpy/kernel_rewrite.py index 339595d7..174f619b 100644 --- a/sumpy/kernel_rewrite.py +++ b/sumpy/kernel_rewrite.py @@ -258,8 +258,6 @@ def rewrite_using_base_kernel_fourier( # {{{ rewrite_using_base_kernel_lu -INT_MAX = 10**15 - class FactorizationFailedError(Exception): pass @@ -434,6 +432,45 @@ def rewrite_using_base_kernel_lu( return LinearOperatorRepresentation(target_kernel, base_kernel, mis, coeffs) +def generate_points( + dim: int, + npoints: int, *, + nshells: int = 2, + rmin: float = 0.25, + rmax: float = 2.0, + rng: np.random.Generator | None = None + ) -> onp.Array2D[Any]: + if dim < 1: + raise ValueError(f"'dim' must be >= 1: {dim!r}") + + if npoints < 1: + raise ValueError(f"'npoints' must be >= 1: {npoints!r}") + + if nshells < 1: + raise ValueError(f"'nshells' must be >= 1: {nshells!r}") + + if rmin >= rmax: + raise ValueError(f"'rmin' must be smaller than 'rmax': {rmin} >= {rmax}") + + if rng is None: + rng = np.random.default_rng() + + # make log spaced shell radii + radii = np.geomspace(rmin, rmax, nshells) + shell = np.arange(npoints) % nshells + + # jitter the radius within each shell so that the points do not all share + # the same (few) radii, which would alias radial modes of the kernels + ratio = (rmax / rmin)**(1 / (nshells - 1)) if nshells > 1 else 1.0 + jitter = rng.uniform(1.0, ratio, npoints) if ratio > 1.0 else 1.0 + + # generate points on the unit sphere + p = rng.standard_normal((dim, npoints)) + p /= np.linalg.norm(p, axis=0) + + return radii[shell] * jitter * p + + def _get_base_kernel_matrix_lu_factorization( base_kernel: ScalarKernel, order: int, @@ -474,11 +511,7 @@ def _get_base_kernel_matrix_lu_factorization( # try to LU factorize on random points for _ in range(retries): - # TODO: is it faster to generate numbers and then sympify them? - points = np.empty((dim, len(mis) + 1), dtype=object) - for i in range(points.shape[0]): - for j in range(points.shape[1]): - points[i, j] = sym.Integer(rng.integers(1, INT_MAX)) / INT_MAX + points = generate_points(dim, len(mis) + 1, rng=rng) # evaluate derivatives at points and construct matrix entries: list[list[sym.Basic]] = [] From 4bddd8bc383c560c85896e875fa46fed97cffef2 Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Fri, 11 Sep 2026 14:34:01 +0300 Subject: [PATCH 04/10] feat: keep fractions if possible --- sumpy/kernel_rewrite.py | 12 +++++++++--- 1 file changed, 9 insertions(+), 3 deletions(-) diff --git a/sumpy/kernel_rewrite.py b/sumpy/kernel_rewrite.py index 174f619b..5970fb59 100644 --- a/sumpy/kernel_rewrite.py +++ b/sumpy/kernel_rewrite.py @@ -294,12 +294,14 @@ def evalf(expr: sym.Expr, prec: int = 100) -> sym.Expr: def round_expr( expr: sym.Basic, atol: float = 1.0e-8, rtol: float = 1.0e-5 ) -> sym.Basic: - """Round all numeric values in *expr* to the nearest integer. + """Round all numeric values in *expr* to the nearest integer or fraction. This function clips all numbers close to zero (effectively removing them - from the expression due to SymPy's automatic symplifications) and rounds all + from the expression due to SymPy's automatic simplifications) and rounds all numbers to the nearest integer using the given *atol* and *rtol*. """ + from fractions import Fraction + nums = expr.atoms(sym.Float) replace_dict: dict[Any, float] = {} @@ -310,7 +312,11 @@ def round_expr( if abs(value - nearest_int) < atol + rtol * abs(value): replace_dict[num] = sym.Integer(nearest_int) else: - replace_dict[num] = value + frac = Fraction(value).limit_denominator(1000) + if abs(float(frac) - value) < atol + rtol * abs(value): + replace_dict[num] = sym.Rational(frac.numerator, frac.denominator) + else: + replace_dict[num] = value return expr.xreplace(replace_dict) From fd667a700bd904f969ec7f093dc8cf3e3d9bc925 Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Fri, 11 Sep 2026 15:15:38 +0300 Subject: [PATCH 05/10] feat: reformat lu matrix generation --- sumpy/kernel_rewrite.py | 142 +++++++++++++++++------------- sumpy/test/test_kernel_rewrite.py | 51 +++++++++++ 2 files changed, 130 insertions(+), 63 deletions(-) diff --git a/sumpy/kernel_rewrite.py b/sumpy/kernel_rewrite.py index 5970fb59..86ba955f 100644 --- a/sumpy/kernel_rewrite.py +++ b/sumpy/kernel_rewrite.py @@ -330,6 +330,82 @@ def simplify(expr: sym.Basic) -> sym.Basic: return sp.simplify(expr) +def _generate_points_shells( + dim: int, + npoints: int, *, + nshells: int = 2, + rmin: float = 0.25, + rmax: float = 2.0, + rng: np.random.Generator | None = None + ) -> onp.Array2D[Any]: + if dim < 1: + raise ValueError(f"'dim' must be >= 1: {dim!r}") + + if npoints < 1: + raise ValueError(f"'npoints' must be >= 1: {npoints!r}") + + if nshells < 1: + raise ValueError(f"'nshells' must be >= 1: {nshells!r}") + + if rmin >= rmax: + raise ValueError(f"'rmin' must be smaller than 'rmax': {rmin} >= {rmax}") + + if rng is None: + rng = np.random.default_rng() + + # make log spaced shell radii + # NOTE: we uniformly sample in each shell to ensure that: + # 1. we're not too close to the origin + # 2. the radii don't repeat to avoid issues with radially symmetric kernels + edges = np.geomspace(rmin, rmax, nshells + 1) + shell = np.arange(npoints) % nshells + log_r = rng.uniform(np.log(edges[shell]), np.log(edges[shell + 1])) + + # generate points on the unit sphere + p = rng.standard_normal((dim, npoints)) + p /= np.linalg.norm(p, axis=0) + + return np.exp(log_r) * p + + +def _make_expr_derivatives( + base_expr: sym.Expr, dvec: sym.Matrix, mis: Sequence[MultiIndex] + ) -> dict[MultiIndex, sym.Basic]: + mi_to_derivative: dict[MultiIndex, sym.Basic] = {} + for mi in mis: + expr = base_expr + for i, nderivs in enumerate(mi): + if nderivs == 0: + continue + expr = expr.diff(dvec[i], nderivs) + + mi_to_derivative[mi] = expr + + return mi_to_derivative + + +def _make_derivative_matrix( + points: onp.Array2D[Any], + dvec: sym.Matrix, + mis: Sequence[MultiIndex], + mi_to_derivative: dict[MultiIndex, sym.Basic] + ) -> sym.Matrix: + # evaluate derivatives at points and construct matrix + entries: list[list[sym.Basic]] = [] + for i in range(points.shape[1]): + row: list[sym.Basic] = [sym.Integer(1)] + + for mi in mis: + expr = mi_to_derivative[mi].xreplace( + dict(zip(dvec, points[:, i], strict=True)) + ) + row.append(evalf(expr)) + + entries.append(row) + + return sym.Matrix(entries) + + def rewrite_using_base_kernel_lu( target_kernel: ScalarKernel, base_kernel: ScalarKernel, @@ -438,45 +514,6 @@ def rewrite_using_base_kernel_lu( return LinearOperatorRepresentation(target_kernel, base_kernel, mis, coeffs) -def generate_points( - dim: int, - npoints: int, *, - nshells: int = 2, - rmin: float = 0.25, - rmax: float = 2.0, - rng: np.random.Generator | None = None - ) -> onp.Array2D[Any]: - if dim < 1: - raise ValueError(f"'dim' must be >= 1: {dim!r}") - - if npoints < 1: - raise ValueError(f"'npoints' must be >= 1: {npoints!r}") - - if nshells < 1: - raise ValueError(f"'nshells' must be >= 1: {nshells!r}") - - if rmin >= rmax: - raise ValueError(f"'rmin' must be smaller than 'rmax': {rmin} >= {rmax}") - - if rng is None: - rng = np.random.default_rng() - - # make log spaced shell radii - radii = np.geomspace(rmin, rmax, nshells) - shell = np.arange(npoints) % nshells - - # jitter the radius within each shell so that the points do not all share - # the same (few) radii, which would alias radial modes of the kernels - ratio = (rmax / rmin)**(1 / (nshells - 1)) if nshells > 1 else 1.0 - jitter = rng.uniform(1.0, ratio, npoints) if ratio > 1.0 else 1.0 - - # generate points on the unit sphere - p = rng.standard_normal((dim, npoints)) - p /= np.linalg.norm(p, axis=0) - - return radii[shell] * jitter * p - - def _get_base_kernel_matrix_lu_factorization( base_kernel: ScalarKernel, order: int, @@ -505,33 +542,12 @@ def _get_base_kernel_matrix_lu_factorization( base_expr = base_kernel.get_expression(dvec) # evaluate all the needed derivatives - mi_to_derivative: dict[MultiIndex, sym.Basic] = {} - for mi in mis: - expr = base_expr - for i, nderivs in enumerate(mi): - if nderivs == 0: - continue - expr = expr.diff(dvec[i], nderivs) - - mi_to_derivative[mi] = expr + mi_to_derivative = _make_expr_derivatives(base_expr, dvec, mis) # try to LU factorize on random points for _ in range(retries): - points = generate_points(dim, len(mis) + 1, rng=rng) - - # evaluate derivatives at points and construct matrix - entries: list[list[sym.Basic]] = [] - for i in range(points.shape[1]): - row: list[sym.Basic] = [sym.Integer(1)] - - for mi in mis: - expr = mi_to_derivative[mi].replace( - dict(zip(dvec, points[:, i], strict=True)) - ) - row.append(evalf(expr)) - - entries.append(row) - mat = sym.Matrix(entries) + points = _generate_points_shells(dim, len(mis) + 1, rng=rng) + mat = _make_derivative_matrix(points, dvec, mis, mi_to_derivative) # TODO: LUdecomposition in symengine is not implemented for non-square matrices try: diff --git a/sumpy/test/test_kernel_rewrite.py b/sumpy/test/test_kernel_rewrite.py index 9d3c42c1..d23547f8 100644 --- a/sumpy/test/test_kernel_rewrite.py +++ b/sumpy/test/test_kernel_rewrite.py @@ -26,6 +26,7 @@ logger = logging.getLogger(__name__) +logging.basicConfig(level=logging.INFO) def mi_derivative(expr: sym.Expr, x: sym.Matrix, mi: tuple[int, ...]) -> sym.Expr: @@ -60,6 +61,56 @@ def check_kernel_rewrite(op: LinearOperatorRepresentation) -> None: assert abs(result) < 3.0 * 1.0e-16 +# {{{ test_rewrite_using_base_kernel_lu_conditioning + + +@pytest.mark.parametrize("dim", [2, 3]) +def test_rewrite_using_base_kernel_lu_conditioning(dim: int) -> None: + from pytools import ( + generate_nonnegative_integer_tuples_summing_to_at_most as gnitstam, + ) + + from sumpy.kernel_rewrite import ( + _generate_points_shells, + _make_derivative_matrix, + _make_expr_derivatives, + ) + + rng = np.random.default_rng(42) + base_kernel = BiharmonicKernel(dim) + + dvec = sym.make_sym_vector("d", dim) + base_expr = base_kernel.get_expression(dvec) + + pde = base_kernel.get_pde_as_diff_op() + mis = list(gnitstam(pde.order, dim)) + pde_mis = [ident.mi for eq in pde.eqs for ident in eq] + pde_mis = [mi for mi in pde_mis if sum(mi) == pde.order] + mis.remove(pde_mis[-1]) + + mi_to_derivative = _make_expr_derivatives(base_expr, dvec, mis) + + nruns = 16 + kappa = np.empty(nruns) + + for i in range(nruns): + points = _generate_points_shells(dim, len(mis) + 1, rng=rng) + mat = _make_derivative_matrix(points, dvec, mis, mi_to_derivative) + + mat = np.array([ + [float(mat[i, j]) for j in range(mat.shape[1])] + for i in range(mat.shape[0]) + ]) + + kappa[i] = np.linalg.cond(mat) + logger.info("kappa = %.8e", kappa[i]) + + logger.info("median: %.8e max %.8e", np.median(kappa), np.max(kappa)) + assert np.max(kappa) < 2.0e+5 + +# }}} + + # {{{ test_rewrite_using_base_kernel_lu_laplace_biharmonic @pytest.mark.parametrize("dim", [2, 3]) From 5ed934e871921a5f6af9ad5d138d1aca4bf02346 Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Fri, 11 Sep 2026 15:33:15 +0300 Subject: [PATCH 06/10] feat: simplify generating points --- sumpy/kernel_rewrite.py | 13 ++----------- sumpy/test/test_kernel_rewrite.py | 2 +- 2 files changed, 3 insertions(+), 12 deletions(-) diff --git a/sumpy/kernel_rewrite.py b/sumpy/kernel_rewrite.py index 86ba955f..7e4e4cef 100644 --- a/sumpy/kernel_rewrite.py +++ b/sumpy/kernel_rewrite.py @@ -333,7 +333,6 @@ def simplify(expr: sym.Basic) -> sym.Basic: def _generate_points_shells( dim: int, npoints: int, *, - nshells: int = 2, rmin: float = 0.25, rmax: float = 2.0, rng: np.random.Generator | None = None @@ -344,9 +343,6 @@ def _generate_points_shells( if npoints < 1: raise ValueError(f"'npoints' must be >= 1: {npoints!r}") - if nshells < 1: - raise ValueError(f"'nshells' must be >= 1: {nshells!r}") - if rmin >= rmax: raise ValueError(f"'rmin' must be smaller than 'rmax': {rmin} >= {rmax}") @@ -354,18 +350,13 @@ def _generate_points_shells( rng = np.random.default_rng() # make log spaced shell radii - # NOTE: we uniformly sample in each shell to ensure that: - # 1. we're not too close to the origin - # 2. the radii don't repeat to avoid issues with radially symmetric kernels - edges = np.geomspace(rmin, rmax, nshells + 1) - shell = np.arange(npoints) % nshells - log_r = rng.uniform(np.log(edges[shell]), np.log(edges[shell + 1])) + radii = np.exp(rng.uniform(np.log(rmin), np.log(rmax), npoints)) # generate points on the unit sphere p = rng.standard_normal((dim, npoints)) p /= np.linalg.norm(p, axis=0) - return np.exp(log_r) * p + return radii * p def _make_expr_derivatives( diff --git a/sumpy/test/test_kernel_rewrite.py b/sumpy/test/test_kernel_rewrite.py index d23547f8..ae251e44 100644 --- a/sumpy/test/test_kernel_rewrite.py +++ b/sumpy/test/test_kernel_rewrite.py @@ -106,7 +106,7 @@ def test_rewrite_using_base_kernel_lu_conditioning(dim: int) -> None: logger.info("kappa = %.8e", kappa[i]) logger.info("median: %.8e max %.8e", np.median(kappa), np.max(kappa)) - assert np.max(kappa) < 2.0e+5 + assert np.max(kappa) < 3.0e+5 # }}} From cd7889ccf0966e20ef7f14e4de49e4d1ba8558b6 Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Fri, 11 Sep 2026 15:49:48 +0300 Subject: [PATCH 07/10] feat: add check for scalar kernel --- sumpy/kernel_rewrite.py | 117 +++++++++++++++++++++------------------- 1 file changed, 63 insertions(+), 54 deletions(-) diff --git a/sumpy/kernel_rewrite.py b/sumpy/kernel_rewrite.py index 7e4e4cef..56efab62 100644 --- a/sumpy/kernel_rewrite.py +++ b/sumpy/kernel_rewrite.py @@ -397,6 +397,61 @@ def _make_derivative_matrix( return sym.Matrix(entries) +def _get_base_kernel_matrix_lu_factorization( + base_kernel: ScalarKernel, + order: int, + *, + rng: np.random.Generator, + retries: int, + ) -> _LUDecomposition: + pde = base_kernel.get_pde_as_diff_op() + if order > pde.order: + raise NotImplementedError( + "Rewriting when the base kernel's derivatives are linearly dependent " + "is not implemented") + + dim = base_kernel.dim + + mis = list(gnitstam(order, dim)) + if order == pde.order: + pde_mis = [ident.mi for eq in pde.eqs for ident in eq] + pde_mis = [mi for mi in pde_mis if sum(mi) == order] + mis.remove(pde_mis[-1]) + + logger.debug("Removing %s to avoid linear dependent mis", pde_mis[-1]) + + # get sympy expression for the base kernel + dvec = sym.make_sym_vector("d", dim) + base_expr = base_kernel.get_expression(dvec) + + # evaluate all the needed derivatives + mi_to_derivative = _make_expr_derivatives(base_expr, dvec, mis) + + # try to LU factorize on random points + for _ in range(retries): + points = _generate_points_shells(dim, len(mis) + 1, rng=rng) + mat = _make_derivative_matrix(points, dvec, mis, mi_to_derivative) + + # NOTE: this can only happen if the points are somehow degenerate, e.g. + # points[:, 0] == points[:, 1], so it shouldn't happen under normal + # operating conditions? + try: + L, U, perm = mat.LUdecomposition() # ruff: ignore[non-lowercase-variable-in-function] + except RuntimeError: + continue + else: + # NOTE: and sympy seems to set the last row of U to 0 + if not sym.USE_SYMENGINE and all(expr == 0 for expr in U[-1, :]): + continue + + return _LUDecomposition(L, U, perm, mis, points) + + raise FactorizationFailedError( + f"failed to compute LU factorization to order {order} for {base_kernel} " + f"after {retries} retries" + ) + + def rewrite_using_base_kernel_lu( target_kernel: ScalarKernel, base_kernel: ScalarKernel, @@ -420,6 +475,14 @@ def rewrite_using_base_kernel_lu( :arg retries: maximum number of retries for each order. If the LU decomposition fails due to a poor choice of random points, it is retried several times. """ + try: + _ = base_kernel.get_pde_system_kernel() + except TypeError: + pass + else: + raise ValueError( + f"'base_kernel' cannot be part of a system: {type(base_kernel)}" + ) pde = base_kernel.get_pde_as_diff_op() if min_order is None: @@ -505,58 +568,4 @@ def rewrite_using_base_kernel_lu( return LinearOperatorRepresentation(target_kernel, base_kernel, mis, coeffs) -def _get_base_kernel_matrix_lu_factorization( - base_kernel: ScalarKernel, - order: int, - *, - rng: np.random.Generator, - retries: int, - ) -> _LUDecomposition: - pde = base_kernel.get_pde_as_diff_op() - if order > pde.order: - raise NotImplementedError( - "Rewriting when the base kernel's derivatives are linearly dependent " - "is not implemented") - - dim = base_kernel.dim - - mis = list(gnitstam(order, dim)) - if order == pde.order: - pde_mis = [ident.mi for eq in pde.eqs for ident in eq] - pde_mis = [mi for mi in pde_mis if sum(mi) == order] - mis.remove(pde_mis[-1]) - - logger.debug("Removing %s to avoid linear dependent mis", pde_mis[-1]) - - # get sympy expression for the base kernel - dvec = sym.make_sym_vector("d", dim) - base_expr = base_kernel.get_expression(dvec) - - # evaluate all the needed derivatives - mi_to_derivative = _make_expr_derivatives(base_expr, dvec, mis) - - # try to LU factorize on random points - for _ in range(retries): - points = _generate_points_shells(dim, len(mis) + 1, rng=rng) - mat = _make_derivative_matrix(points, dvec, mis, mi_to_derivative) - - # TODO: LUdecomposition in symengine is not implemented for non-square matrices - try: - L, U, perm = mat.LUdecomposition() # ruff: ignore[non-lowercase-variable-in-function] - except RuntimeError: - # NOTE: symengine seems to throw a SymEngineError -> RuntimeError when - # it fails to do the LU factorization due to rank-deficiency - continue - else: - # NOTE: and sympy seems to set the last row of U to 0 - if not sym.USE_SYMENGINE and all(expr == 0 for expr in U[-1, :]): - continue - - return _LUDecomposition(L, U, perm, mis, points) - - raise FactorizationFailedError( - f"failed to compute LU factorization to order {order} for {base_kernel} " - f"after {retries} retries" - ) - # }}} From fc5d837502e5f7500dc45c251e783b4b6eb4b807 Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Fri, 11 Sep 2026 16:26:11 +0300 Subject: [PATCH 08/10] feat: add proper check for match --- sumpy/kernel_rewrite.py | 69 ++++++++++++++++++++++++++++++++++------- 1 file changed, 58 insertions(+), 11 deletions(-) diff --git a/sumpy/kernel_rewrite.py b/sumpy/kernel_rewrite.py index 56efab62..6f7bdd66 100644 --- a/sumpy/kernel_rewrite.py +++ b/sumpy/kernel_rewrite.py @@ -397,6 +397,36 @@ def _make_derivative_matrix( return sym.Matrix(entries) +def _check_linear_combination( + target_expr: sym.Expr, + base_expr: sym.Expr, + dvec: sym.Matrix, + mis: Sequence[MultiIndex], + coeffs: Sequence[sym.Basic], + *, + rng: np.random.Generator, + rtol: float = 1.0e-8, + ) -> bool: + dim = len(dvec) + mi_to_derivative = _make_expr_derivatives(base_expr, dvec, mis) + points = _generate_points_shells(dim, len(mis) + 1, rng=rng) + + max_lhs = 0.0 + max_err = 0.0 + for i in range(points.shape[1]): + subst = dict(zip(dvec, points[:, i], strict=True)) + lhs = float(evalf(target_expr.xreplace(subst))) + + rhs = float(evalf(coeffs[0])) + for c, mi in zip(coeffs[1:], mis, strict=True): + rhs += float(evalf(c * mi_to_derivative[mi].xreplace(subst))) + + max_lhs = max(max_lhs, abs(lhs)) + max_err = max(max_err, abs(lhs - rhs)) + + return max_err <= rtol * max(max_lhs, 1.0) + + def _get_base_kernel_matrix_lu_factorization( base_kernel: ScalarKernel, order: int, @@ -502,12 +532,12 @@ def rewrite_using_base_kernel_lu( dim = base_kernel.dim dvec = sym.make_sym_vector("d", dim) target_expr = target_kernel.get_expression(dvec) + base_expr = base_kernel.get_expression(dvec) target_scaling = target_kernel.get_global_scaling_const() base_scaling = base_kernel.get_global_scaling_const() order = min_order - to_pymbolic = sym.SympyToPymbolicMapperWithSymbols() while order <= pde.order: try: lu = _get_base_kernel_matrix_lu_factorization( @@ -535,7 +565,7 @@ def rewrite_using_base_kernel_lu( postprocess=lambda x: x.expand()) # gather all non-zero coefficients from the result - const = 0 + const = sym.Integer(0) coeffs = [] mis = [] for i, coeff in enumerate(all_coeffs): @@ -544,20 +574,28 @@ def rewrite_using_base_kernel_lu( continue if i == 0: - const = to_pymbolic(simplify(coeff * target_scaling)) - logger.debug(" %s", coeff) + const = coeff + logger.debug(" %s", const) else: - mi = lu.mis[i - 1] - coeff = simplify(coeff * target_scaling / base_scaling) - - mis.append(mi) - coeffs.append(to_pymbolic(coeff)) - logger.debug(" + %s*%s.diff%s", coeff, base_kernel, mi) + mis.append(lu.mis[i - 1]) + coeffs.append(coeff) + logger.debug(" + %s*%s.diff%s", coeff, base_kernel, lu.mis[-1]) if coeffs: coeffs.insert(0, const) + + success = _check_linear_combination( + target_expr, + base_expr, + dvec, + mis, + coeffs, + rng=rng, + ) + if coeffs and success: break + coeffs = [] order += 1 if not coeffs: @@ -565,7 +603,16 @@ def rewrite_using_base_kernel_lu( f"could not express {target_kernel} in terms of {base_kernel}" ) - return LinearOperatorRepresentation(target_kernel, base_kernel, mis, coeffs) + to_pymbolic = sym.SympyToPymbolicMapperWithSymbols() + return LinearOperatorRepresentation( + target_kernel, base_kernel, mis, + tuple( + to_pymbolic(simplify( + c * (target_scaling if i == 0 else (target_scaling / base_scaling)) + )) + for i, c in enumerate(coeffs) + ), + ) # }}} From baf577f3a2d29d7aef507dd291ddf850d4da22d5 Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Fri, 11 Sep 2026 19:59:28 +0300 Subject: [PATCH 09/10] feat: add more tests --- sumpy/test/test_kernel_rewrite.py | 42 +++++++++++++++++++++++++++++++ 1 file changed, 42 insertions(+) diff --git a/sumpy/test/test_kernel_rewrite.py b/sumpy/test/test_kernel_rewrite.py index ae251e44..b44e6977 100644 --- a/sumpy/test/test_kernel_rewrite.py +++ b/sumpy/test/test_kernel_rewrite.py @@ -14,9 +14,11 @@ import sumpy.symbolic as sym from sumpy.kernel import ( BiharmonicKernel, + BrinkmanletComponentKernel, LaplaceKernel, StokesletComponentKernel, StressletComponentKernel, + YukawaKernel, ) from sumpy.kernel_rewrite import ( LinearOperatorRepresentation, @@ -146,6 +148,26 @@ def test_rewrite_using_base_kernel_lu_stokeslet_biharmonic(dim: int) -> None: # }}} +# {{{ test_rewrite_using_base_kernel_lu_brinkman_yukawa + +@pytest.mark.parametrize("dim", [2, 3]) +def test_rewrite_using_base_kernel_lu_brinkman_yukawa(dim: int) -> None: + from pytools import generate_nonnegative_integer_tuples_below as gnitb + + rng = np.random.default_rng(seed=42) + + base_kernel = YukawaKernel(dim, yukawa_lambda_name="k") + for i, j in gnitb(dim, 2): + target_kernel = BrinkmanletComponentKernel( + dim, i, j, viscosity_mu_name="mu", darcy_impermeability_name="k") + result = rewrite_using_base_kernel_lu(target_kernel, base_kernel, rng=rng) + print(result.pretty()) + check_kernel_rewrite(result) + + +# }}} + + # {{{ test_rewrite_using_base_kernel_lu_stresslet_biharmonic @pytest.mark.parametrize("dim", [2, 3]) @@ -202,6 +224,26 @@ def test_rewrite_using_base_kernel_fourier_stokeslet_biharmonic(dim: int) -> Non # }}} +# {{{ test_rewrite_using_base_kernel_fourier_brinkman_yukawa + + +@pytest.mark.parametrize("dim", [2, 3]) +def test_rewrite_using_base_kernel_fourier_brinkman_yukawa(dim: int) -> None: + from itertools import product + + base_kernel = YukawaKernel(dim, yukawa_lambda_name="k") + + for i, j in product(range(dim), repeat=2): + target_kernel = BrinkmanletComponentKernel( + dim, i, j, viscosity_mu_name="mu", darcy_impermeability_name="k") + result = rewrite_using_base_kernel_fourier(target_kernel, base_kernel) + + logger.info(result.pretty()) + check_kernel_rewrite(result) + +# }}} + + # {{{ test_rewrite_using_base_kernel_fourier_indivisible From db038c57d1dc3ff55660c07eb90840346b8fa6ec Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Fri, 18 Sep 2026 17:07:16 +0300 Subject: [PATCH 10/10] feat: enable elasticity kernel --- sumpy/kernel_rewrite.py | 20 ++++++++++++-------- sumpy/test/test_kernel_rewrite.py | 29 +++++++++++++++++++++++++---- 2 files changed, 37 insertions(+), 12 deletions(-) diff --git a/sumpy/kernel_rewrite.py b/sumpy/kernel_rewrite.py index 6f7bdd66..d7e37e51 100644 --- a/sumpy/kernel_rewrite.py +++ b/sumpy/kernel_rewrite.py @@ -415,14 +415,19 @@ def _check_linear_combination( max_err = 0.0 for i in range(points.shape[1]): subst = dict(zip(dvec, points[:, i], strict=True)) - lhs = float(evalf(target_expr.xreplace(subst))) + lhs = evalf(target_expr.xreplace(subst)) - rhs = float(evalf(coeffs[0])) + rhs = evalf(coeffs[0]) for c, mi in zip(coeffs[1:], mis, strict=True): - rhs += float(evalf(c * mi_to_derivative[mi].xreplace(subst))) + rhs += evalf(c * mi_to_derivative[mi].xreplace(subst)) - max_lhs = max(max_lhs, abs(lhs)) - max_err = max(max_err, abs(lhs - rhs)) + if lhs.is_number and rhs.is_number: + max_lhs = max(max_lhs, abs(lhs)) + max_err = max(max_err, abs(lhs - rhs)) + else: + max_err = round_expr(simplify(lhs - rhs), rtol=rtol) + if not max_err.is_zero: + return False return max_err <= rtol * max(max_lhs, 1.0) @@ -603,13 +608,12 @@ def rewrite_using_base_kernel_lu( f"could not express {target_kernel} in terms of {base_kernel}" ) - to_pymbolic = sym.SympyToPymbolicMapperWithSymbols() return LinearOperatorRepresentation( target_kernel, base_kernel, mis, tuple( - to_pymbolic(simplify( + sym.to_pymbolic(simplify( c * (target_scaling if i == 0 else (target_scaling / base_scaling)) - )) + ), symbols=True) for i, c in enumerate(coeffs) ), ) diff --git a/sumpy/test/test_kernel_rewrite.py b/sumpy/test/test_kernel_rewrite.py index b44e6977..5a9f800b 100644 --- a/sumpy/test/test_kernel_rewrite.py +++ b/sumpy/test/test_kernel_rewrite.py @@ -15,6 +15,7 @@ from sumpy.kernel import ( BiharmonicKernel, BrinkmanletComponentKernel, + ElasticityComponentKernel, LaplaceKernel, StokesletComponentKernel, StressletComponentKernel, @@ -41,11 +42,9 @@ def mi_derivative(expr: sym.Expr, x: sym.Matrix, mi: tuple[int, ...]) -> sym.Exp def check_kernel_rewrite(op: LinearOperatorRepresentation) -> None: from sumpy.kernel_rewrite import evalf, simplify - from sumpy.symbolic import PymbolicToSympyMapperWithSymbols dim = op.target_kernel.dim dvec = sym.make_sym_vector("d", dim) - to_sympy = PymbolicToSympyMapperWithSymbols() target_expr = ( op.target_kernel.get_global_scaling_const() @@ -54,8 +53,8 @@ def check_kernel_rewrite(op: LinearOperatorRepresentation) -> None: op.base_kernel.get_global_scaling_const() * op.base_kernel.get_expression(dvec)) - expr = to_sympy(op.coeffs[0]) + sum(( - to_sympy(c) * mi_derivative(base_expr, dvec, mi) + expr = sym.to_symbolic(op.coeffs[0], symbols=True) + sum(( + sym.to_symbolic(c, symbols=True) * mi_derivative(base_expr, dvec, mi) for c, mi in zip(op.coeffs[1:], op.mis, strict=True) ), sym.Integer(0)) @@ -150,6 +149,7 @@ def test_rewrite_using_base_kernel_lu_stokeslet_biharmonic(dim: int) -> None: # {{{ test_rewrite_using_base_kernel_lu_brinkman_yukawa +@pytest.mark.skip(reason="needs two bases (Laplace and Yukawa)") @pytest.mark.parametrize("dim", [2, 3]) def test_rewrite_using_base_kernel_lu_brinkman_yukawa(dim: int) -> None: from pytools import generate_nonnegative_integer_tuples_below as gnitb @@ -187,6 +187,27 @@ def test_rewrite_using_base_kernel_lu_stresslet_biharmonic(dim: int) -> None: # }}} +# {{{ test_rewrite_using_base_kernel_lu_elasticity_biharmonic + + +@pytest.mark.parametrize("dim", [2, 3]) +def test_rewrite_using_base_kernel_lu_elasticity_biharmonic(dim: int) -> None: + from pytools import generate_nonnegative_integer_tuples_below as gnitb + + rng = np.random.default_rng(seed=42) + + base_kernel = BiharmonicKernel(dim) + for i, j in gnitb(dim, 2): + target_kernel = ElasticityComponentKernel( + dim, i, j, viscosity_mu_name="mu", poisson_ratio_name="nu" + ) + result = rewrite_using_base_kernel_lu(target_kernel, base_kernel, rng=rng) + print(result.pretty()) + check_kernel_rewrite(result) + +# }}} + + # {{{ test_rewrite_using_base_kernel_fourier_laplace_biharmonic