Support for kernels with no known base kernel

sumpy/derivative_taker.py

@@ -36,12 +36,16 @@
 """
 
 from pytools.tag import tag_dataclass
+from pytools import (single_valued,
+    generate_nonnegative_integer_tuples_summing_to_at_most as gnitstam)
+from math import floor
 
 import numpy as np
 import sumpy.symbolic as sym
-from sumpy.tools import add_to_sac, add_mi
+from sumpy.tools import add_to_sac, add_mi, nullspace
+from sumpy.kernel import Kernel
 
-from typing import Dict, Tuple, Any
+from typing import Dict, Tuple, Any, Mapping, Text
 
 import logging
 
@@ -404,6 +408,172 @@ def diff_derivative_coeff_dict(derivative_coeff_dict: DerivativeCoeffDict,
     return {derivative: coeff for derivative, coeff in
             new_derivative_coeff_dict.items() if coeff != 0}
 
+
+def _get_sympy_kernel_expression(kernel: Kernel,
+        kernel_arguments: Mapping[Text, Any]) -> sym.Basic:
+    """Convert a :mod:`pymbolic` expression to :mod:`sympy` expression
+    after substituting kernel arguments.
+    For eg: `exp(I*k*r)/r` with `{k: 1}` is converted to the sympy expression
+    `exp(I*r)/r`
+    """
+    from pymbolic.mapper.substitutor import substitute
+    from sumpy.symbolic import PymbolicToSympyMapperWithSymbols
+
+    expr = substitute(kernel.get_base_kernel().expression, kernel_arguments)
+    expr = PymbolicToSympyMapperWithSymbols()(expr)
+
+    dvec = sym.make_sym_vector("d", kernel.dim)
+    res = kernel.postprocess_at_target(kernel.postprocess_at_source(
+        expr, dvec), dvec)
+    return res
+
+
+def evalf(expr: sym.Basic, dps: float):
+    """evaluate an expression numerically using ``dps``
+    number of digits.
+    """
+    from sumpy.symbolic import USE_SYMENGINE
+    if USE_SYMENGINE:
+        import symengine
+        prec = int(symengine.log(10**dps, 2))
+        return expr.n(prec=prec)
+    else:
+        return expr.n(n=dps)
+
+
+def chop(expr: sym.Basic, tol: float) -> sym.Basic:
+    """Given a symbolic expression, remove all occurences of numbers
+    with absolute value less than a given tolerance and replace floating
+    point numbers that are close to an integer up to a given relative
+    tolerance by the integer.
+    """
+    nums = expr.atoms(sym.Number)
+    replace_dict = {}
+    for num in nums:
+        if float(abs(num)) < tol:
+            replace_dict[num] = 0
+        else:
+            new_num = float(num)
+            if abs((int(new_num) - new_num)/new_num) < tol:
+                new_num = int(new_num)
+            replace_dict[num] = new_num
+    return expr.xreplace(replace_dict)
+
+
+def get_deriv_sample(kernel, order, samples, kernel_arguments, atol):
+    dim = kernel.dim
+    sym_vec = sym.make_sym_vector("d", dim)
+    base_expr = _get_sympy_kernel_expression(kernel,
+        dict(kernel_arguments))
+
+    mis = sorted(gnitstam(order, dim), key=sum)
+    assert samples.shape[0] == dim
+
+    exprs = []
+    for mi in mis:
+        expr = base_expr
+        for var_idx, nderivs in enumerate(mi):
+            if nderivs == 0:
+                continue
+            expr = expr.diff(sym_vec[var_idx], nderivs)
+        exprs.append(expr)
+
+    dps = -sym.log(atol, 10)
+    mat = []
+    for isample in range(samples.shape[1]):
+        row = []
+        for ideriv in range(len(mis)):
+            expr = exprs[ideriv]
+            replace_dict = dict(zip(sym_vec, samples[:, isample]))
+            eval_expr = evalf(expr.xreplace(replace_dict), dps)
+            row.append(eval_expr)
+        mat.append(row)
+    mat = sym.Matrix(mat)
+
+    return mat, mis
+
+
+def get_dependent_columns(matrix, atol):
+    import sympy
+    m = matrix.T
+    l, u, p = sympy.Matrix(m).LUdecomposition(
+            iszerofunc=lambda x: abs(x) < atol)
+    nrows = m.shape[0]
+    idxs = list(range(nrows))
+    for i, j in p:
+        idxs[i], idxs[j] = idxs[j], idxs[i]
+
+    nonzero_rows = 0
+    for i in range(nrows - 1, -1, -1):
+        if not all(abs(elem) < atol for elem in u[i, :]):
+            nonzero_rows = i + 1
+            break
+
+    return idxs[nonzero_rows:]
+
+
+def get_pde_operators(kernels, order, kernel_arguments, atol=1e-30):
+    import sympy
+    from sumpy.expansion.diff_op import diff, make_identity_diff_op
+    dim = single_valued(kernel.dim for kernel in kernels)
+
+    mis = sorted(gnitstam(order, dim), key=sum)
+
+    # (-1, -1, -1) represents a constant
+    # ((0,0,0) would be "function with no derivatives")
+    # mis.append((-1,)*dim)
+
+    n = len(kernels)
+    nsamples = int(floor(len(mis) * n / (n-1))) + n + 1
+    rand = np.random.randint(1, 10**15, (dim, nsamples))
+    rand = rand.astype(object)
+    for i in range(rand.shape[0]):
+        for j in range(rand.shape[1]):
+            rand[i, j] = sym.sympify(rand[i, j])/10**15
+
+    derivs_evaluated = [
+        get_deriv_sample(kernel, order, rand, kernel_arguments, atol)[0]
+        for kernel in kernels]
+
+    for mat in derivs_evaluated:
+        dep_cols = get_dependent_columns(mat, atol * 1e10)
+        zeros = [0]*mat.shape[0]
+        for col in dep_cols:
+            mat[:, col] = zeros
+
+    full_mat = sym.zeros((n - 1) * nsamples, len(mis) * n)
+    assert full_mat.shape[0] > full_mat.shape[1]
+    for i in range(1, n):
+        full_mat[(i - 1)*nsamples:i*nsamples, :len(mis)] = derivs_evaluated[i]
+        full_mat[(i - 1)*nsamples:i*nsamples, i*len(mis):(i + 1)*len(mis)] = \
+                -derivs_evaluated[0]
+
+    ns = nullspace(full_mat.tolist(), atol * 1e10)
+    for col in range(ns.shape[1]):
+        for i in range(n):
+            if all(abs(elem) < atol * 1e10 for elem in
+                    ns[i*len(mis):(i + 1)*len(mis), col]):
+                break
+        else:
+            ops = sym.Matrix(ns[:, col].tolist()).reshape(n, len(mis))
+            ops = ops.applyfunc(lambda x: sympy.nsimplify(x, tolerance=atol*1e10))
+            id_op = make_identity_diff_op(dim, 1)
+            diff_ops = []
+            for i in range(n):
+                diff_op = None
+                for mi_idx, coeff in enumerate(ops[i, :]):
+                    if coeff == 0:
+                        continue
+                    mi = mis[mi_idx]
+                    if not diff_op:
+                        diff_op = coeff * diff(id_op, mi)
+                    else:
+                        diff_op += coeff * diff(id_op, mi)
+                diff_ops.append(diff_op)
+            return diff_ops
+
+    raise RuntimeError("Could not find PDE operators")
+
 # }}}
 
 # vim: fdm=marker

sumpy/expansion/diff_op.py

@@ -60,7 +60,8 @@ def __init__(self, dim, *eqs):
                   to a coefficient.
         """
         self.dim = dim
-        self.eqs = tuple(eqs)
+        self.eqs = tuple([pmap(dict(filter(lambda p: p[1] != 0, eq.items())))
+                for eq in tuple(eqs)])
 
     def __eq__(self, other):
         return self.dim == other.dim and self.eqs == other.eqs
@@ -105,6 +106,9 @@ def __add__(self, other_diff_op):
     def __sub__(self, other_diff_op):
         return self + (-1)*other_diff_op
 
+    def __neg__(self):
+        return self*(-1)
+
     def __repr__(self):
         return f"LinearPDESystemOperator({self.dim}, {repr(self.eqs)})"
 

sumpy/symbolic.py

@@ -111,7 +111,9 @@ def _find_symbolic_backend():
 Symbol = sym.Symbol
 Derivative = sym.Derivative
 Integer = sym.Integer
+Rational = sym.Rational
 Matrix = sym.Matrix
+zeros = sym.zeros
 Subs = sym.Subs
 I = sym.I  # noqa: E741
 pi = sym.pi

sumpy/tools.py

@@ -528,6 +528,11 @@ def nullspace(m, atol=0):
                     vec[piv_col] -= int(mat[piv_row, pos])
                 else:
                     vec[piv_col] -= mat[piv_row, pos]
+        for i in range(len(vec)):
+            if isinstance(vec[i], sym.Basic) and vec[i].is_number and \
+                    abs(vec[i]) < atol:
+                vec[i] = 0
+
         n.append(vec)
     return np.array(n, dtype=object).T
 

sumpy/typing.py

@@ -0,0 +1,31 @@
+__copyright__ = "Copyright (C) Isuru Fernando"
+
+__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.
+"""
+
+
+from typing import Union
+import numpy as np
+from pymbolic.primitives import Expression
+
+IntegralT = Union[int, np.integer]
+FloatT = Union[float, complex, np.floating, np.complexfloating]
+
+ExpressionT = Union[IntegralT, FloatT, Expression]

test/test_misc.py

@@ -44,10 +44,12 @@
     StressletKernel,
     ElasticityKernel,
     LineOfCompressionKernel,
+    AxisTargetDerivative,
     ExpressionKernel)
 from sumpy.expansion.diff_op import (
     make_identity_diff_op, concat, as_scalar_pde, diff,
     gradient, divergence, laplacian, curl)
+from sumpy.derivative_taker import get_pde_operators
 
 import logging
 logger = logging.getLogger(__name__)
@@ -516,6 +518,42 @@ def get_pde_as_diff_op(self):
 # }}}
 
 
+# {{{ test_get_pde_operators
+
+def test_get_pde_operators_laplace_biharmonic():
+    dim = 3
+    laplace = LaplaceKernel(dim)
+    biharmonic = BiharmonicKernel(dim)
+    id_op = make_identity_diff_op(dim, 1)
+
+    op1, op2 = get_pde_operators([laplace, biharmonic], 2, {})
+    assert op1 == laplacian(id_op) * sym.Rational(1, 2)
+    assert op2 == id_op
+
+    d_biharmonic = AxisTargetDerivative(1, biharmonic)
+    op1, op2, op3 = get_pde_operators(
+            [laplace, biharmonic, d_biharmonic], 2, {})
+    assert op1 == laplacian(id_op) * sym.Rational(1, 2)
+    assert op2 == id_op
+    assert op3 == diff(id_op, [0, 1, 0][:dim])
+
+
+def test_get_pde_operators_stokes():
+    for dim in (2, 3):
+        stokes00 = StokesletKernel(dim, 0, 0)
+        stokes01 = StokesletKernel(dim, 0, 1)
+        stokes11 = StokesletKernel(dim, 1, 1)
+
+        id_op = make_identity_diff_op(dim, 1)
+
+        op1, op2, op3 = get_pde_operators([stokes00, stokes01, stokes11], 2, {})
+
+        assert op1 == laplacian(id_op) - diff(id_op, [2, 0, 0][:dim])
+        assert op2 == -diff(id_op, [1, 1, 0][:dim])
+        assert op3 == laplacian(id_op) - diff(id_op, [0, 2, 0][:dim])
+# }}}
+
+
 # You can test individual routines by typing
 # $ python test_misc.py 'test_pde_check_kernels(_acf,
 #       KernelInfo(HelmholtzKernel(2), k=5), order=5)'