From 4e9d1eb15d3feb1c04808fdfaedaed47d74e03c2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andreas=20Kl=C3=B6ckner?= Date: Sat, 25 Feb 2017 19:49:52 -0500 Subject: [PATCH 1/5] Revert "Merge branch 'revert-a63f699d' into 'master'" This reverts merge request !31 --- .gitlab-ci.yml | 12 ++ .test-conda-env-py3-requirements.txt | 3 + .test-conda-env-py3.yml | 16 +++ doc/misc.rst | 13 +++ requirements.txt | 1 + sumpy/assignment_collection.py | 40 +------ sumpy/codegen.py | 159 +++++++++++++++++++++++---- sumpy/cse.py | 61 ++++++---- sumpy/e2e.py | 17 +-- sumpy/e2p.py | 21 ++-- sumpy/expansion/__init__.py | 4 +- sumpy/expansion/local.py | 73 ++++++------ sumpy/expansion/multipole.py | 49 +++++---- sumpy/kernel.py | 11 +- sumpy/p2e.py | 13 +-- sumpy/p2p.py | 5 +- sumpy/qbx.py | 20 ++-- sumpy/symbolic.py | 123 +++++++++++---------- sumpy/tools.py | 56 +++++++++- sumpy/version.py | 2 +- test/test_codegen.py | 32 ++++++ test/test_cse.py | 63 +++++++---- test/test_fmm.py | 1 + 23 files changed, 524 insertions(+), 271 deletions(-) create mode 100644 .test-conda-env-py3-requirements.txt create mode 100644 .test-conda-env-py3.yml diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index 6ee07ae4..63768430 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -66,6 +66,18 @@ Python 3.6 POCL: except: - tags +Python 3.5 Conda: + script: + - export SUMPY_NO_CACHE=1 + - CONDA_ENVIRONMENT=.test-conda-env-py3.yml + - REQUIREMENTS_TXT=.test-conda-env-py3-requirements.txt + - curl -L -O -k https://gitlab.tiker.net/inducer/ci-support/raw/master/build-and-test-py-project-within-miniconda.sh + - ". ./build-and-test-py-project-within-miniconda.sh" + tags: + - linux + except: + - tags + Documentation: script: - EXTRA_INSTALL="numpy mako" diff --git a/.test-conda-env-py3-requirements.txt b/.test-conda-env-py3-requirements.txt new file mode 100644 index 00000000..01f08662 --- /dev/null +++ b/.test-conda-env-py3-requirements.txt @@ -0,0 +1,3 @@ +git+https://github.com/inducer/boxtree +git+https://github.com/inducer/pymbolic +git+https://github.com/inducer/loopy diff --git a/.test-conda-env-py3.yml b/.test-conda-env-py3.yml new file mode 100644 index 00000000..f336c9e9 --- /dev/null +++ b/.test-conda-env-py3.yml @@ -0,0 +1,16 @@ +name: test-conda-env-py3 +channels: +- inducer +- symengine/label/dev +- conda-forge +- defaults +dependencies: +- git +- conda-forge::numpy +- conda-forge::sympy +- pocl +- islpy +- pyopencl +- python=3.5 +- python-symengine=0.2.0.53.g83912b7=py35_1 +# things not in here: loopy boxtree pymbolic pyfmmlib diff --git a/doc/misc.rst b/doc/misc.rst index c6ea75c3..c77c7c22 100644 --- a/doc/misc.rst +++ b/doc/misc.rst @@ -19,6 +19,19 @@ and say:: In addition, you need to have :mod:`numpy` installed. +Symbolic backends +================= + +:mod:`sumpy` supports two symbolic backends: sympy and SymEngine. To use the +SymEngine backend, ensure that the `SymEngine library +`_ and the `SymEngine Python bindings +`_ are installed. + +By default, :mod:`sumpy` prefers using SymEngine but falls back to sympy if it +detects that SymEngine is not installed. To force the use of a particular +backend, set the environment variable `SUMPY_FORCE_SYMBOLIC_BACKEND` to +`symengine` or `sympy`. + User-visible Changes ==================== diff --git a/requirements.txt b/requirements.txt index abe39e89..6e23ae95 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1,5 +1,6 @@ numpy sympy==1.0 +git+https://github.com/inducer/pymbolic git+https://github.com/inducer/islpy git+https://github.com/pyopencl/pyopencl git+https://github.com/inducer/boxtree diff --git a/sumpy/assignment_collection.py b/sumpy/assignment_collection.py index c1e9197b..a12a17fd 100644 --- a/sumpy/assignment_collection.py +++ b/sumpy/assignment_collection.py @@ -26,7 +26,7 @@ THE SOFTWARE. """ -import sympy as sp +import sumpy.symbolic as sym import logging logger = logging.getLogger(__name__) @@ -70,7 +70,7 @@ class _SymbolGenerator(object): self.base_to_count[base] = count + 1 - return sp.Symbol(id_str) + return sym.Symbol(id_str) def __iter__(self): return self @@ -132,7 +132,7 @@ class SymbolicAssignmentCollection(object): result = set() for dep in self.assignments[var_name].atoms(): - if not isinstance(dep, sp.Symbol): + if not isinstance(dep, sym.Symbol): continue dep_name = dep.name @@ -153,7 +153,7 @@ class SymbolicAssignmentCollection(object): if root_name is None: root_name = name - self.assignments[name] = sp.sympify(expr) + self.assignments[name] = sym.sympify(expr) def assign_unique(self, name_base, expr): """Assign *expr* to a new variable whose name is based on *name_base*. @@ -177,7 +177,7 @@ class SymbolicAssignmentCollection(object): # Options here: # - checked_cse: if you mistrust the result of the cse. # Uses maxima to verify. - # - sp.cse: The sympy thing. + # - sym.cse: The sympy thing. # - sumpy.cse.cse: Based on sympy, designed to go faster. #from sumpy.symbolic import checked_cse @@ -192,41 +192,13 @@ class SymbolicAssignmentCollection(object): self.assignments[name] = new_expr for name, value in new_assignments: - assert isinstance(name, sp.Symbol) + assert isinstance(name, sym.Symbol) self.add_assignment(name.name, value) logger.info("common subexpression elimination: done after {dur:.2f} s" .format(dur=time.time() - start_time)) return new_extra_exprs - def kill_trivial_assignments(self, exprs): - logger.info("kill trivial assignments: start") - - approved_assignments = [] - rejected_assignments = [] - - from sumpy.symbolic import is_assignment_nontrivial - for name, value in six.iteritems(self.assignments): - if name in self.user_symbols or is_assignment_nontrivial(name, value): - approved_assignments.append((name, value)) - else: - rejected_assignments.append((name, value)) - - # un-substitute rejected assignments - from sumpy.symbolic import make_one_step_subst - unsubst_rej = make_one_step_subst(rejected_assignments) - - new_assignments = dict( - (name, expr.subs(unsubst_rej)) - for name, expr in approved_assignments) - - exprs = [expr.subs(unsubst_rej) for expr in exprs] - self.assignments = new_assignments - - logger.info("kill trivial assignments: done") - - return exprs - # }}} # vim: fdm=marker diff --git a/sumpy/codegen.py b/sumpy/codegen.py index 0cb14028..99e4f47d 100644 --- a/sumpy/codegen.py +++ b/sumpy/codegen.py @@ -29,6 +29,7 @@ import pyopencl as cl import pyopencl.tools # noqa import loopy as lp +import six import re from pymbolic.mapper import IdentityMapper, WalkMapper, CSECachingMapperMixin @@ -38,8 +39,7 @@ from loopy.types import NumpyType from pytools import memoize_method -from pymbolic.interop.sympy import ( - SympyToPymbolicMapper as SympyToPymbolicMapperBase) +from sumpy.symbolic import (SympyToPymbolicMapper as SympyToPymbolicMapperBase) import logging logger = logging.getLogger(__name__) @@ -58,25 +58,138 @@ Conversion of :mod:`sympy` expressions to :mod:`loopy` # {{{ sympy -> pymbolic mapper +import sumpy.symbolic as sym +_SPECIAL_FUNCTION_NAMES = frozenset(dir(sym.functions)) + + class SympyToPymbolicMapper(SympyToPymbolicMapperBase): - def __init__(self, assignments): - self.assignments = assignments - self.derivative_cse_names = set() - def map_Derivative(self, expr): # noqa - # Sympy has picked up the habit of picking arguments out of derivatives - # and pronounce them common subexpressions. Me no like. Undo it, so - # that the bessel substitutor further down can do its job. + def not_supported(self, expr): + if isinstance(expr, int): + return expr + elif getattr(expr, "is_Function", False): + func_name = SympyToPymbolicMapperBase.function_name(self, expr) + # SymEngine capitalizes the names of the special functions. + if func_name.lower() in _SPECIAL_FUNCTION_NAMES: + func_name = func_name.lower() + return prim.Variable(func_name)( + *tuple(self.rec(arg) for arg in expr.args)) + else: + return SympyToPymbolicMapperBase.not_supported(self, expr) + +# }}} + + +# {{{ trivial assignment elimination + +def make_one_step_subst(assignments): + assignments = dict(assignments) + unwanted_vars = set(six.iterkeys(assignments)) + + # Ensure no re-assignments. + assert len(unwanted_vars) == len(assignments) + + from loopy.symbolic import get_dependencies + unwanted_deps = dict( + (name, get_dependencies(value) & unwanted_vars) + for name, value in six.iteritems(assignments)) + + # {{{ compute substitution order + + toposort = [] + visited = set() + visiting = set() + + while unwanted_vars: + stack = [unwanted_vars.pop()] + + while stack: + top = stack[-1] + + if top in visiting: + visiting.remove(top) + toposort.append(top) + + if top in visited: + stack.pop() + continue + + visited.add(top) + visiting.add(top) + + for dep in unwanted_deps[top]: + # Check for no cycles. + assert dep not in visiting + stack.append(dep) - if expr.expr.is_Symbol: - # These will get removed, because loopy wont' be able to deal - # with them--they contain undefined placeholder symbols. - self.derivative_cse_names.add(expr.expr.name) + # }}} - return prim.Derivative(self.rec( - expr.expr.subs(self.assignments)), - tuple(v.name for v in expr.variables)) + # {{{ make substitution + from pymbolic import substitute + + result = {} + used_name_to_var = {} + from pymbolic import evaluate + from functools import partial + simplify = partial(evaluate, context=used_name_to_var) + + for name in toposort: + value = assignments[name] + value = substitute(value, result) + used_name_to_var.update( + (used_name, prim.Variable(used_name)) + for used_name in get_dependencies(value) + if used_name not in used_name_to_var) + + result[name] = simplify(value) + + # }}} + + return result + + +def is_assignment_nontrivial(name, value): + if prim.is_constant(value): + return False + elif isinstance(value, prim.Variable): + return False + elif (isinstance(value, prim.Product) + and len(value.children) == 2 + and sum(1 for arg in value.children if prim.is_constant(arg)) == 1 + and sum(1 for arg in value.children + if isinstance(arg, prim.Variable)) == 1): + # const*var: not good enough + return False + + return True + + +def kill_trivial_assignments(assignments, retain_names=set()): + logger.info("kill trivial assignments (plain): start") + approved_assignments = [] + rejected_assignments = [] + + for name, value in assignments: + if name in retain_names or is_assignment_nontrivial(name, value): + approved_assignments.append((name, value)) + else: + rejected_assignments.append((name, value)) + + # un-substitute rejected assignments + unsubst_rej = make_one_step_subst(rejected_assignments) + + result = [] + from pymbolic import substitute + for name, expr in approved_assignments: + r = substitute(expr, unsubst_rej) + result.append((name, r)) + + logger.info( + "kill trivial assignments (plain): done, {nrej} assignments killed" + .format(nrej=len(rejected_assignments))) + + return result # }}} @@ -318,7 +431,7 @@ class BesselDerivativeReplacer(CSECachingMapperMixin, IdentityMapper): arg, = expr.values n_derivs = len(expr.child.variables) - import sympy as sp + import sympy as sym # AS (9.1.31) # http://dlmf.nist.gov/10.6.7 @@ -329,7 +442,7 @@ class BesselDerivativeReplacer(CSECachingMapperMixin, IdentityMapper): k = n_derivs return prim.CommonSubexpression( 2**(-k)*sum( - (-1)**idx*int(sp.binomial(k, idx)) * function(i, arg) + (-1)**idx*int(sym.binomial(k, idx)) * function(i, arg) for idx, i in enumerate(range(order-k, order+k+1, 2))), "d%d_%s_%s" % (n_derivs, function.name, order_str)) else: @@ -538,17 +651,15 @@ class MathConstantRewriter(CSECachingMapperMixin, IdentityMapper): def to_loopy_insns(assignments, vector_names=set(), pymbolic_expr_maps=[], - complex_dtype=None): + complex_dtype=None, retain_names=set()): logger.info("loopy instruction generation: start") assignments = list(assignments) # convert from sympy - sympy_conv = SympyToPymbolicMapper(assignments) + sympy_conv = SympyToPymbolicMapper() assignments = [(name, sympy_conv(expr)) for name, expr in assignments] - assignments = [ - (name, expr) for name, expr in assignments - if name not in sympy_conv.derivative_cse_names - ] + + assignments = kill_trivial_assignments(assignments, retain_names) bdr = BesselDerivativeReplacer() assignments = [(name, bdr(expr)) for name, expr in assignments] diff --git a/sumpy/cse.py b/sumpy/cse.py index ad44fc41..573c78c9 100644 --- a/sumpy/cse.py +++ b/sumpy/cse.py @@ -66,8 +66,8 @@ DAMAGE. # }}} -from sympy.core import Basic, Mul, Add, Pow -from sympy.core.function import _coeff_isneg +from sumpy.symbolic import ( + Basic, Mul, Add, Pow, Symbol, _coeff_isneg, Derivative, Subs) from sympy.core.compatibility import iterable from sympy.utilities.iterables import numbered_symbols @@ -82,6 +82,10 @@ Common subexpression elimination """ +# Don't CSE child nodes of these classes. +CSE_NO_DESCEND_CLASSES = (Derivative, Subs) + + # {{{ cse pre/postprocessing def preprocess_for_cse(expr, optimizations): @@ -248,6 +252,19 @@ class FuncArgTracker(object): self.func_to_argset[func_i].update(new_args) +class Unevaluated(object): + + def __init__(self, func, args): + self.func = func + self.args = args + + def __str__(self): + return "Uneval<{}>({})".format( + self.func, ", ".join(str(a) for a in self.args)) + + __repr__ = __str__ + + def match_common_args(func_class, funcs, opt_subs): """ Recognize and extract common subexpressions of function arguments within a @@ -305,9 +322,8 @@ def match_common_args(func_class, funcs, opt_subs): diff_i = arg_tracker.func_to_argset[i].difference(com_args) if diff_i: # com_func needs to be unevaluated to allow for recursive matches. - com_func = func_class( - *arg_tracker.get_args_in_value_order(com_args), - evaluate=False) + com_func = Unevaluated( + func_class, arg_tracker.get_args_in_value_order(com_args)) com_func_number = arg_tracker.get_or_add_value_number(com_func) arg_tracker.update_func_argset(i, diff_i | set([com_func_number])) changed.add(i) @@ -334,9 +350,8 @@ def match_common_args(func_class, funcs, opt_subs): changed.add(k) if i in changed: - opt_subs[funcs[i]] = func_class( - *arg_tracker.get_args_in_value_order(arg_tracker.func_to_argset[i]), - evaluate=False) + opt_subs[funcs[i]] = Unevaluated(func_class, + arg_tracker.get_args_in_value_order(arg_tracker.func_to_argset[i])) arg_tracker.stop_arg_tracking(i) @@ -366,6 +381,9 @@ def opt_cse(exprs): if expr.is_Atom: return + if isinstance(expr, CSE_NO_DESCEND_CLASSES): + return + if iterable(expr): for item in expr: find_opts(item) @@ -382,7 +400,7 @@ def opt_cse(exprs): if _coeff_isneg(expr): neg_expr = -expr if not neg_expr.is_Atom: - opt_subs[expr] = Mul(-1, neg_expr, evaluate=False) + opt_subs[expr] = Unevaluated(Mul, (-1, neg_expr)) seen_subexp.add(neg_expr) expr = neg_expr @@ -395,7 +413,7 @@ def opt_cse(exprs): elif isinstance(expr, Pow): base, exp = expr.args if _coeff_isneg(exp): - opt_subs[expr] = Pow(Pow(base, -exp), -1, evaluate=False) + opt_subs[expr] = Unevaluated(Pow, (Pow(base, -exp), -1)) # }}} @@ -436,10 +454,10 @@ def tree_cse(exprs, symbols, opt_subs=None): excluded_symbols = set() def find_repeated(expr): - if not isinstance(expr, Basic): + if not isinstance(expr, (Basic, Unevaluated)): return - if expr.is_Atom: + if isinstance(expr, Basic) and expr.is_Atom: if expr.is_Symbol: excluded_symbols.add(expr) return @@ -457,7 +475,10 @@ def tree_cse(exprs, symbols, opt_subs=None): if expr in opt_subs: expr = opt_subs[expr] - args = expr.args + if isinstance(expr, CSE_NO_DESCEND_CLASSES): + args = () + else: + args = expr.args for arg in args: find_repeated(arg) @@ -478,7 +499,7 @@ def tree_cse(exprs, symbols, opt_subs=None): subs = dict() def rebuild(expr): - if not isinstance(expr, Basic): + if not isinstance(expr, (Basic, Unevaluated)): return expr if not expr.args: @@ -495,11 +516,11 @@ def tree_cse(exprs, symbols, opt_subs=None): if expr in opt_subs: expr = opt_subs[expr] - new_args = [rebuild(arg) for arg in expr.args] - if new_args != expr.args: - new_expr = expr.func(*new_args) - else: - new_expr = expr + new_expr = expr + if not isinstance(expr, CSE_NO_DESCEND_CLASSES): + new_args = tuple(rebuild(arg) for arg in expr.args) + if isinstance(expr, Unevaluated) or new_args != expr.args: + new_expr = expr.func(*new_args) if orig_expr in to_eliminate: try: @@ -560,7 +581,7 @@ def cse(exprs, symbols=None, optimizations=None): reduced_exprs = [preprocess_for_cse(e, optimizations) for e in exprs] if symbols is None: - symbols = numbered_symbols() + symbols = numbered_symbols(cls=Symbol) else: # In case we get passed an iterable with an __iter__ method instead of # an actual iterator. diff --git a/sumpy/e2e.py b/sumpy/e2e.py index f5b34940..66e5e975 100644 --- a/sumpy/e2e.py +++ b/sumpy/e2e.py @@ -27,7 +27,7 @@ from six.moves import range import numpy as np import loopy as lp -import sympy as sp +import sumpy.symbolic as sym from sumpy.tools import KernelCacheWrapper import logging @@ -77,10 +77,10 @@ class E2EBase(KernelCacheWrapper): self.dim = src_expansion.dim def get_translation_loopy_insns(self): - from sumpy.symbolic import make_sympy_vector - dvec = make_sympy_vector("d", self.dim) + from sumpy.symbolic import make_sym_vector + dvec = make_sym_vector("d", self.dim) - src_coeff_exprs = [sp.Symbol("src_coeff%d" % i) + src_coeff_exprs = [sym.Symbol("src_coeff%d" % i) for i in range(len(self.src_expansion))] from sumpy.assignment_collection import SymbolicAssignmentCollection @@ -93,17 +93,12 @@ class E2EBase(KernelCacheWrapper): sac.run_global_cse() - from sumpy.symbolic import kill_trivial_assignments - assignments = kill_trivial_assignments([ - (name, expr) - for name, expr in six.iteritems(sac.assignments)], - retain_names=tgt_coeff_names) - from sumpy.codegen import to_loopy_insns return to_loopy_insns( - assignments, + six.iteritems(sac.assignments), vector_names=set(["d"]), pymbolic_expr_maps=[self.tgt_expansion.get_code_transformer()], + retain_names=tgt_coeff_names, complex_dtype=np.complex128 # FIXME ) diff --git a/sumpy/e2p.py b/sumpy/e2p.py index b2e0e585..6a77a5cc 100644 --- a/sumpy/e2p.py +++ b/sumpy/e2p.py @@ -27,7 +27,7 @@ from six.moves import range import numpy as np import loopy as lp -import sympy as sp +import sumpy.symbolic as sym from sumpy.tools import KernelCacheWrapper @@ -74,15 +74,16 @@ class E2PBase(KernelCacheWrapper): assert tdr(knl) == expansion.kernel def get_loopy_insns_and_result_names(self): - from sumpy.symbolic import make_sympy_vector - bvec = make_sympy_vector("b", self.dim) + from sumpy.symbolic import make_sym_vector + bvec = make_sym_vector("b", self.dim) from sumpy.assignment_collection import SymbolicAssignmentCollection sac = SymbolicAssignmentCollection() - coeff_exprs = [sp.Symbol("coeff%d" % i) + coeff_exprs = [sym.Symbol("coeff%d" % i) for i in range(len(self.expansion.get_coefficient_identifiers()))] value = self.expansion.evaluate(coeff_exprs, bvec) + result_names = [ sac.assign_unique("result_%d_p" % i, knl.postprocess_at_target(value, bvec)) @@ -91,23 +92,19 @@ class E2PBase(KernelCacheWrapper): sac.run_global_cse() - from sumpy.symbolic import kill_trivial_assignments - assignments = kill_trivial_assignments([ - (name, expr) - for name, expr in six.iteritems(sac.assignments)], - retain_names=result_names) - from sumpy.codegen import to_loopy_insns - loopy_insns = to_loopy_insns(assignments, + loopy_insns = to_loopy_insns( + six.iteritems(sac.assignments), vector_names=set(["b"]), pymbolic_expr_maps=[self.expansion.get_code_transformer()], + retain_names=result_names, complex_dtype=np.complex128 # FIXME ) return loopy_insns, result_names def get_kernel_scaling_assignment(self): - from pymbolic.interop.sympy import SympyToPymbolicMapper + from sumpy.symbolic import SympyToPymbolicMapper sympy_conv = SympyToPymbolicMapper() return [lp.Assignment(id=None, assignee="kernel_scaling", diff --git a/sumpy/expansion/__init__.py b/sumpy/expansion/__init__.py index cc49c6ca..f50bd5fb 100644 --- a/sumpy/expansion/__init__.py +++ b/sumpy/expansion/__init__.py @@ -25,9 +25,9 @@ THE SOFTWARE. """ import numpy as np -import sympy as sp import logging from pytools import memoize_method +import sumpy.symbolic as sym from sumpy.tools import MiDerivativeTaker __doc__ = """ @@ -356,7 +356,7 @@ class HelmholtzDerivativeWrangler(LinearRecurrenceBasedDerivativeWrangler): coeffs[needed_deriv] = -1 else: - k = sp.Symbol(self.helmholtz_k_name) + k = sym.Symbol(self.helmholtz_k_name) coeffs[tuple(reduced_deriv)] = -k*k return coeffs diff --git a/sumpy/expansion/local.py b/sumpy/expansion/local.py index 3ba643cd..f9ec106a 100644 --- a/sumpy/expansion/local.py +++ b/sumpy/expansion/local.py @@ -23,7 +23,7 @@ OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. """ -import sympy as sp +import sumpy.symbolic as sym from sumpy.expansion import ( ExpansionBase, VolumeTaylorExpansion, LaplaceConformingVolumeTaylorExpansion, @@ -61,24 +61,29 @@ class LineTaylorLocalExpansion(LocalExpansionBase): raise RuntimeError("cannot use line-Taylor expansions in a setting " "where the center-target vector is not known at coefficient " "formation") - avec_line = avec + sp.Symbol("tau")*bvec + + tau = sym.Symbol("tau") + + avec_line = avec + tau*bvec line_kernel = self.kernel.get_expression(avec_line) - return [ - self.kernel.postprocess_at_target( - self.kernel.postprocess_at_source( - line_kernel.diff("tau", i), - avec), - bvec) - .subs("tau", 0) + from sumpy.tools import MiDerivativeTaker, my_syntactic_subs + deriv_taker = MiDerivativeTaker(line_kernel, (tau,)) + + return [my_syntactic_subs( + self.kernel.postprocess_at_target( + self.kernel.postprocess_at_source( + deriv_taker.diff(i), + avec), bvec), + {tau: 0}) for i in self.get_coefficient_identifiers()] def evaluate(self, coeffs, bvec): from pytools import factorial - return sum( + return sym.Add(*( coeffs[self.get_storage_index(i)] / factorial(i) - for i in self.get_coefficient_identifiers()) + for i in self.get_coefficient_identifiers())) # }}} @@ -140,7 +145,7 @@ class VolumeTaylorLocalExpansionBase(LocalExpansionBase): src_expansion.get_coefficient_identifiers()): kernel_deriv = taker.diff(add_mi(deriv, term)) local_result.append(coeff * kernel_deriv) - result.append(sp.Add(*local_result)) + result.append(sym.Add(*local_result)) else: from sumpy.tools import MiDerivativeTaker expr = src_expansion.evaluate(src_coeff_exprs, dvec) @@ -197,61 +202,61 @@ class H2DLocalExpansion(LocalExpansionBase): return list(range(-self.order, self.order+1)) def coefficients_from_source(self, avec, bvec): - from sumpy.symbolic import sympy_real_norm_2 - hankel_1 = sp.Function("hankel_1") + from sumpy.symbolic import sym_real_norm_2 + hankel_1 = sym.Function("hankel_1") - k = sp.Symbol(self.kernel.get_base_kernel().helmholtz_k_name) + k = sym.Symbol(self.kernel.get_base_kernel().helmholtz_k_name) # The coordinates are negated since avec points from source to center. - source_angle_rel_center = sp.atan2(-avec[1], -avec[0]) - avec_len = sympy_real_norm_2(avec) + source_angle_rel_center = sym.atan2(-avec[1], -avec[0]) + avec_len = sym_real_norm_2(avec) return [self.kernel.postprocess_at_source( hankel_1(l, k * avec_len) - * sp.exp(sp.I * l * source_angle_rel_center), avec) + * sym.exp(sym.I * l * source_angle_rel_center), avec) for l in self.get_coefficient_identifiers()] def evaluate(self, coeffs, bvec): - from sumpy.symbolic import sympy_real_norm_2 - bessel_j = sp.Function("bessel_j") - bvec_len = sympy_real_norm_2(bvec) - target_angle_rel_center = sp.atan2(bvec[1], bvec[0]) + from sumpy.symbolic import sym_real_norm_2 + bessel_j = sym.Function("bessel_j") + bvec_len = sym_real_norm_2(bvec) + target_angle_rel_center = sym.atan2(bvec[1], bvec[0]) - k = sp.Symbol(self.kernel.get_base_kernel().helmholtz_k_name) + k = sym.Symbol(self.kernel.get_base_kernel().helmholtz_k_name) return sum(coeffs[self.get_storage_index(l)] * self.kernel.postprocess_at_target( bessel_j(l, k * bvec_len) - * sp.exp(sp.I * l * -target_angle_rel_center), bvec) + * sym.exp(sym.I * l * -target_angle_rel_center), bvec) for l in self.get_coefficient_identifiers()) def translate_from(self, src_expansion, src_coeff_exprs, dvec): - from sumpy.symbolic import sympy_real_norm_2 + from sumpy.symbolic import sym_real_norm_2 - k = sp.Symbol(self.kernel.get_base_kernel().helmholtz_k_name) + k = sym.Symbol(self.kernel.get_base_kernel().helmholtz_k_name) if isinstance(src_expansion, H2DLocalExpansion): - dvec_len = sympy_real_norm_2(dvec) - bessel_j = sp.Function("bessel_j") - new_center_angle_rel_old_center = sp.atan2(dvec[1], dvec[0]) + dvec_len = sym_real_norm_2(dvec) + bessel_j = sym.Function("bessel_j") + new_center_angle_rel_old_center = sym.atan2(dvec[1], dvec[0]) translated_coeffs = [] for l in self.get_coefficient_identifiers(): translated_coeffs.append( sum(src_coeff_exprs[src_expansion.get_storage_index(m)] * bessel_j(m - l, k * dvec_len) - * sp.exp(sp.I * (m - l) * -new_center_angle_rel_old_center) + * sym.exp(sym.I * (m - l) * -new_center_angle_rel_old_center) for m in src_expansion.get_coefficient_identifiers())) return translated_coeffs from sumpy.expansion.multipole import H2DMultipoleExpansion if isinstance(src_expansion, H2DMultipoleExpansion): - dvec_len = sympy_real_norm_2(dvec) - hankel_1 = sp.Function("hankel_1") - new_center_angle_rel_old_center = sp.atan2(dvec[1], dvec[0]) + dvec_len = sym_real_norm_2(dvec) + hankel_1 = sym.Function("hankel_1") + new_center_angle_rel_old_center = sym.atan2(dvec[1], dvec[0]) translated_coeffs = [] for l in self.get_coefficient_identifiers(): translated_coeffs.append( sum((-1) ** l * hankel_1(m + l, k * dvec_len) - * sp.exp(sp.I * (m + l) * new_center_angle_rel_old_center) + * sym.exp(sym.I * (m + l) * new_center_angle_rel_old_center) * src_coeff_exprs[src_expansion.get_storage_index(m)] for m in src_expansion.get_coefficient_identifiers())) return translated_coeffs diff --git a/sumpy/expansion/multipole.py b/sumpy/expansion/multipole.py index fa227e89..ce3b78d8 100644 --- a/sumpy/expansion/multipole.py +++ b/sumpy/expansion/multipole.py @@ -23,7 +23,7 @@ THE SOFTWARE. """ from six.moves import range, zip -import sympy as sp # noqa +import sumpy.symbolic as sym # noqa from sumpy.expansion import ( ExpansionBase, VolumeTaylorExpansion, LaplaceConformingVolumeTaylorExpansion, @@ -63,8 +63,8 @@ class VolumeTaylorMultipoleExpansionBase(MultipoleExpansionBase): raise NotImplementedError("more than one source derivative " "not supported at present") - from sumpy.symbolic import make_sympy_vector - dir_vec = make_sympy_vector(kernel.dir_vec_name, kernel.dim) + from sumpy.symbolic import make_sym_vector + dir_vec = make_sym_vector(kernel.dir_vec_name, kernel.dim) coeff_identifiers = self.get_full_coefficient_identifiers() result = [0] * len(coeff_identifiers) @@ -91,9 +91,9 @@ class VolumeTaylorMultipoleExpansionBase(MultipoleExpansionBase): def evaluate(self, coeffs, bvec): taker = self.get_kernel_derivative_taker(bvec) - result = sum( + result = sym.Add(*tuple( coeff * taker.diff(mi) - for coeff, mi in zip(coeffs, self.get_coefficient_identifiers())) + for coeff, mi in zip(coeffs, self.get_coefficient_identifiers()))) return result def get_kernel_derivative_taker(self, bvec): @@ -143,7 +143,8 @@ class VolumeTaylorMultipoleExpansionBase(MultipoleExpansionBase): n = tgt_mi[idim] k = src_mi[idim] assert n >= k - contrib *= (sp.binomial(n, k) + from sympy import binomial + contrib *= (binomial(n, k) * dvec[idim]**(n-k)) result[i] += contrib @@ -201,31 +202,31 @@ class H2DMultipoleExpansion(MultipoleExpansionBase): return list(range(-self.order, self.order+1)) def coefficients_from_source(self, avec, bvec): - from sumpy.symbolic import sympy_real_norm_2 - bessel_j = sp.Function("bessel_j") - avec_len = sympy_real_norm_2(avec) + from sumpy.symbolic import sym_real_norm_2 + bessel_j = sym.Function("bessel_j") + avec_len = sym_real_norm_2(avec) - k = sp.Symbol(self.kernel.get_base_kernel().helmholtz_k_name) + k = sym.Symbol(self.kernel.get_base_kernel().helmholtz_k_name) # The coordinates are negated since avec points from source to center. - source_angle_rel_center = sp.atan2(-avec[1], -avec[0]) + source_angle_rel_center = sym.atan2(-avec[1], -avec[0]) return [self.kernel.postprocess_at_source( bessel_j(l, k * avec_len) * - sp.exp(sp.I * l * -source_angle_rel_center), avec) + sym.exp(sym.I * l * -source_angle_rel_center), avec) for l in self.get_coefficient_identifiers()] def evaluate(self, coeffs, bvec): - from sumpy.symbolic import sympy_real_norm_2 - hankel_1 = sp.Function("hankel_1") - bvec_len = sympy_real_norm_2(bvec) - target_angle_rel_center = sp.atan2(bvec[1], bvec[0]) + from sumpy.symbolic import sym_real_norm_2 + hankel_1 = sym.Function("hankel_1") + bvec_len = sym_real_norm_2(bvec) + target_angle_rel_center = sym.atan2(bvec[1], bvec[0]) - k = sp.Symbol(self.kernel.get_base_kernel().helmholtz_k_name) + k = sym.Symbol(self.kernel.get_base_kernel().helmholtz_k_name) return sum(coeffs[self.get_storage_index(l)] * self.kernel.postprocess_at_target( hankel_1(l, k * bvec_len) - * sp.exp(sp.I * l * target_angle_rel_center), bvec) + * sym.exp(sym.I * l * target_angle_rel_center), bvec) for l in self.get_coefficient_identifiers()) def translate_from(self, src_expansion, src_coeff_exprs, dvec): @@ -233,19 +234,19 @@ class H2DMultipoleExpansion(MultipoleExpansionBase): raise RuntimeError("do not know how to translate %s to " "multipole 2D Helmholtz Bessel expansion" % type(src_expansion).__name__) - from sumpy.symbolic import sympy_real_norm_2 - dvec_len = sympy_real_norm_2(dvec) - bessel_j = sp.Function("bessel_j") - new_center_angle_rel_old_center = sp.atan2(dvec[1], dvec[0]) + from sumpy.symbolic import sym_real_norm_2 + dvec_len = sym_real_norm_2(dvec) + bessel_j = sym.Function("bessel_j") + new_center_angle_rel_old_center = sym.atan2(dvec[1], dvec[0]) - k = sp.Symbol(self.kernel.get_base_kernel().helmholtz_k_name) + k = sym.Symbol(self.kernel.get_base_kernel().helmholtz_k_name) translated_coeffs = [] for l in self.get_coefficient_identifiers(): translated_coeffs.append( sum(src_coeff_exprs[src_expansion.get_storage_index(m)] * bessel_j(m - l, k * dvec_len) - * sp.exp(sp.I * (m - l) * new_center_angle_rel_old_center) + * sym.exp(sym.I * (m - l) * new_center_angle_rel_old_center) for m in src_expansion.get_coefficient_identifiers())) return translated_coeffs diff --git a/sumpy/kernel.py b/sumpy/kernel.py index e29ac516..ea0ea02b 100644 --- a/sumpy/kernel.py +++ b/sumpy/kernel.py @@ -264,10 +264,11 @@ class ExpressionKernel(Kernel): if self.dim != len(dist_vec): raise ValueError("dist_vec length does not match expected dimension") - expr = expr.subs([ - ("d%d" % i, dist_vec_i) + from sumpy.symbolic import Symbol + expr = expr.subs(dict( + (Symbol("d%d" % i), dist_vec_i) for i, dist_vec_i in enumerate(dist_vec) - ]) + )) return expr @@ -715,7 +716,7 @@ class DirectionalTargetDerivative(DirectionalDerivative): dim = len(bvec) assert dim == self.dim - from sumpy.symbolic import make_sympy_vector + from sumpy.symbolic import make_sym_vector as make_sympy_vector dir_vec = make_sympy_vector(self.dir_vec_name, dim) # bvec = tgt-center @@ -746,7 +747,7 @@ class DirectionalSourceDerivative(DirectionalDerivative): dimensions = len(avec) assert dimensions == self.dim - from sumpy.symbolic import make_sympy_vector + from sumpy.symbolic import make_sym_vector as make_sympy_vector dir_vec = make_sympy_vector(self.dir_vec_name, dimensions) # avec = center-src -> minus sign from chain rule diff --git a/sumpy/p2e.py b/sumpy/p2e.py index 1b037dd4..6a8fcfa9 100644 --- a/sumpy/p2e.py +++ b/sumpy/p2e.py @@ -70,8 +70,8 @@ class P2EBase(KernelCacheWrapper): self.dim = expansion.dim def get_loopy_instructions(self): - from sumpy.symbolic import make_sympy_vector - avec = make_sympy_vector("a", self.dim) + from sumpy.symbolic import make_sym_vector + avec = make_sym_vector("a", self.dim) from sumpy.assignment_collection import SymbolicAssignmentCollection sac = SymbolicAssignmentCollection() @@ -83,17 +83,12 @@ class P2EBase(KernelCacheWrapper): sac.run_global_cse() - from sumpy.symbolic import kill_trivial_assignments - assignments = kill_trivial_assignments([ - (name, expr) - for name, expr in six.iteritems(sac.assignments)], - retain_names=coeff_names) - from sumpy.codegen import to_loopy_insns return to_loopy_insns( - assignments, + six.iteritems(sac.assignments), vector_names=set(["a"]), pymbolic_expr_maps=[self.expansion.get_code_transformer()], + retain_names=coeff_names, complex_dtype=np.complex128 # FIXME ) diff --git a/sumpy/p2p.py b/sumpy/p2p.py index 3071a3b8..39a90205 100644 --- a/sumpy/p2p.py +++ b/sumpy/p2p.py @@ -69,8 +69,8 @@ class P2PBase(KernelComputation, KernelCacheWrapper): self.dim = single_valued(knl.dim for knl in self.kernels) def get_loopy_insns_and_result_names(self): - from sumpy.symbolic import make_sympy_vector - dvec = make_sympy_vector("d", self.dim) + from sumpy.symbolic import make_sym_vector + dvec = make_sym_vector("d", self.dim) from sumpy.assignment_collection import SymbolicAssignmentCollection sac = SymbolicAssignmentCollection() @@ -90,6 +90,7 @@ class P2PBase(KernelComputation, KernelCacheWrapper): vector_names=set(["d"]), pymbolic_expr_maps=[ knl.get_code_transformer() for knl in self.kernels], + retain_names=result_names, complex_dtype=np.complex128 # FIXME ) diff --git a/sumpy/qbx.py b/sumpy/qbx.py index 0b3dc035..8892739e 100644 --- a/sumpy/qbx.py +++ b/sumpy/qbx.py @@ -27,7 +27,7 @@ import six from six.moves import range, zip import numpy as np import loopy as lp -import sympy as sp +import sumpy.symbolic as sym from pytools import memoize_method from pymbolic import parse, var @@ -64,7 +64,7 @@ def expand(expansion_nr, sac, expansion, avec, bvec): coefficients = expansion.coefficients_from_source(avec, bvec) assigned_coeffs = [ - sp.Symbol( + sym.Symbol( sac.assign_unique("expn%dcoeff%s" % ( expansion_nr, stringify_expn_index(i)), coefficients[expansion.get_storage_index(i)])) @@ -117,10 +117,10 @@ class LayerPotentialBase(KernelComputation, KernelCacheWrapper): @memoize_method def get_kernel(self): - from sumpy.symbolic import make_sympy_vector + from sumpy.symbolic import make_sym_vector - avec = make_sympy_vector("a", self.dim) - bvec = make_sympy_vector("b", self.dim) + avec = make_sym_vector("a", self.dim) + bvec = make_sym_vector("b", self.dim) from sumpy.assignment_collection import SymbolicAssignmentCollection sac = SymbolicAssignmentCollection() @@ -134,17 +134,13 @@ class LayerPotentialBase(KernelComputation, KernelCacheWrapper): sac.run_global_cse() - from sumpy.symbolic import kill_trivial_assignments - assignments = kill_trivial_assignments([ - (name, expr.subs("tau", 0)) - for name, expr in six.iteritems(sac.assignments)], - retain_names=result_names) - from sumpy.codegen import to_loopy_insns - loopy_insns = to_loopy_insns(assignments, + loopy_insns = to_loopy_insns( + six.iteritems(sac.assignments), vector_names=set(["a", "b"]), pymbolic_expr_maps=[ expn.kernel.get_code_transformer() for expn in self.expansions], + retain_names=result_names, complex_dtype=np.complex128 # FIXME ) diff --git a/sumpy/symbolic.py b/sumpy/symbolic.py index 86c2666d..85dbb278 100644 --- a/sumpy/symbolic.py +++ b/sumpy/symbolic.py @@ -27,68 +27,77 @@ import six from six.moves import range from six.moves import zip -import sympy as sp import numpy as np from pymbolic.mapper import IdentityMapper as IdentityMapperBase -from pymbolic.interop.sympy import PymbolicToSympyMapper import pymbolic.primitives as prim import logging logger = logging.getLogger(__name__) -# {{{ trivial assignment elimination +# {{{ symbolic backend -def make_one_step_subst(assignments): - unwanted_vars = set(sp.Symbol(name) for name, value in assignments) +def _find_symbolic_backend(): + global USE_SYMENGINE - result = {} - for name, value in assignments: - while value.atoms() & unwanted_vars: - value = value.subs(assignments) + try: + import symengine # noqa + symengine_found = True + except ImportError: + symengine_found = False - result[sp.Symbol(name)] = value - - return result + ALLOWED_BACKENDS = ("sympy", "symengine") # noqa + BACKEND_ENV_VAR = "SUMPY_FORCE_SYMBOLIC_BACKEND" # noqa + import os + backend = os.environ.get(BACKEND_ENV_VAR) + if backend is not None: + if backend not in ALLOWED_BACKENDS: + raise RuntimeError( + "%s value is unrecognized: '%s' " + "(allowed values are %s)" % ( + BACKEND_ENV_VAR, + backend, + ", ".join("'%s'" % val for val in ALLOWED_BACKENDS))) -def is_assignment_nontrivial(name, value): - if value.is_Number: - return False - elif isinstance(value, sp.Symbol): - return False - elif (isinstance(value, sp.Mul) - and len(value.args) == 2 - and sum(1 for arg in value.args if arg.is_Number) == 1 - and sum(1 for arg in value.args if isinstance(arg, sp.Symbol)) == 1): - # const*var: not good enough - return False + if backend == "symengine" and not symengine_found: + from warnings import warn + warn("%s=symengine was specified, but could not find symengine. " + "Using sympy." % BACKEND_ENV_VAR, RuntimeWarning) - return True + USE_SYMENGINE = backend == "symengine" and symengine_found + else: + USE_SYMENGINE = symengine_found -def kill_trivial_assignments(assignments, retain_names=set()): - logger.info("kill trivial assignments (plain): start") - approved_assignments = [] - rejected_assignments = [] +_find_symbolic_backend() - for name, value in assignments: - if name in retain_names or is_assignment_nontrivial(name, value): - approved_assignments.append((name, value)) - else: - rejected_assignments.append((name, value)) +# }}} - # un-substitute rejected assignments - unsubst_rej = make_one_step_subst(rejected_assignments) +# Symbolic API common to SymEngine and sympy. +# Before adding a function here, make sure it's present in both modules. +SYMBOLIC_API = """ +Add Basic Mul Pow exp sqrt symbols sympify cos sin atan2 Function Symbol +Derivative Integer Matrix Subs I pi functions""".split() - result = [(name, expr.xreplace(unsubst_rej)) - for name, expr in approved_assignments] +if USE_SYMENGINE: + from symengine import sympy_compat as sym + from pymbolic.interop.symengine import ( + PymbolicToSymEngineMapper as PymbolicToSympyMapper, + SymEngineToPymbolicMapper as SympyToPymbolicMapper) +else: + import sympy as sym + from pymbolic.interop.sympy import ( + PymbolicToSympyMapper, SympyToPymbolicMapper) - logger.info("kill trivial assignments (plain): done") +for _apifunc in SYMBOLIC_API: + globals()[_apifunc] = getattr(sym, _apifunc) - return result -# }}} +def _coeff_isneg(a): + if a.is_Mul: + a = a.args[0] + return a.is_Number and a.is_negative # {{{ debugging of sympy CSE via Maxima @@ -111,7 +120,6 @@ def _get_assignments_in_maxima(assignments, prefix=""): from pymbolic.maxima import MaximaStringifyMapper mstr = MaximaStringifyMapper() - from pymbolic.interop.sympy import SympyToPymbolicMapper s2p = SympyToPymbolicMapper() dkill = _DerivativeKiller() @@ -119,12 +127,12 @@ def _get_assignments_in_maxima(assignments, prefix=""): def write_assignment(name): symbols = [atm for atm in assignments[name].atoms() - if isinstance(atm, sp.Symbol) + if isinstance(atm, sym.Symbol) and atm.name in my_variable_names] - for sym in symbols: - if sym.name not in written_assignments: - write_assignment(sym.name) + for symb in symbols: + if symb.name not in written_assignments: + write_assignment(symb.name) result.append("%s%s : %s;" % ( prefix, name, mstr(dkill(s2p( @@ -143,7 +151,7 @@ def checked_cse(exprs, symbols=None): if symbols is not None: kwargs["symbols"] = symbols - new_assignments, new_exprs = sp.cse(exprs, **kwargs) + new_assignments, new_exprs = sym.cse(exprs, **kwargs) max_old = _get_assignments_in_maxima(dict( ("old_expr%d" % i, expr) @@ -171,8 +179,8 @@ def checked_cse(exprs, symbols=None): # }}} -def sympy_real_norm_2(x): - return sp.sqrt((x.T*x)[0, 0]) +def sym_real_norm_2(x): + return sym.sqrt((x.T*x)[0, 0]) def pymbolic_real_norm_2(x): @@ -180,9 +188,9 @@ def pymbolic_real_norm_2(x): return var("sqrt")(np.dot(x, x)) -def make_sympy_vector(name, components): - return sp.Matrix( - [sp.Symbol("%s%d" % (name, i)) for i in range(components)]) +def make_sym_vector(name, components): + return sym.Matrix( + [sym.Symbol("%s%d" % (name, i)) for i in range(components)]) def vector_subs(expr, old, new): @@ -194,8 +202,8 @@ def vector_subs(expr, old, new): def find_power_of(base, prod): - remdr = sp.Wild("remdr") - power = sp.Wild("power") + remdr = sym.Wild("remdr") + power = sym.Wild("power") result = prod.match(remdr*base**power) if result is None: return 0 @@ -205,17 +213,16 @@ def find_power_of(base, prod): class PymbolicToSympyMapperWithSymbols(PymbolicToSympyMapper): def map_variable(self, expr): if expr.name == "I": - return sp.I + return sym.I elif expr.name == "pi": - return sp.pi + return sym.pi else: return PymbolicToSympyMapper.map_variable(self, expr) def map_subscript(self, expr): if isinstance(expr.aggregate, prim.Variable) and isinstance(expr.index, int): - return sp.Symbol("%s%d" % (expr.aggregate.name, expr.index)) + return sym.Symbol("%s%d" % (expr.aggregate.name, expr.index)) else: - raise RuntimeError("do not know how to translate '%s' to sympy" - % expr) + self.raise_conversion_error(expr) # vim: fdm=marker diff --git a/sumpy/tools.py b/sumpy/tools.py index 60687510..f5014295 100644 --- a/sumpy/tools.py +++ b/sumpy/tools.py @@ -28,6 +28,7 @@ THE SOFTWARE. from pytools import memoize_method import numpy as np +import sumpy.symbolic as sym import logging logger = logging.getLogger(__name__) @@ -57,6 +58,7 @@ def mi_power(vector, mi): class MiDerivativeTaker(object): def __init__(self, expr, var_list): + assert isinstance(expr, sym.Basic) self.var_list = var_list empty_mi = (0,) * len(var_list) self.cache_by_mi = {empty_mi: expr} @@ -107,7 +109,7 @@ class LinearRecurrenceBasedMiDerivativeTaker(MiDerivativeTaker): expr = self.cache_by_mi[mi] except KeyError: from six import iteritems - from sympy import Add + from sumpy.symbolic import Add closest_mi = self.get_closest_cached_mi(mi) expr = self.cache_by_mi[closest_mi] @@ -264,7 +266,7 @@ class KernelComputation(object): self.name = name or self.default_name def get_kernel_scaling_assignments(self): - from pymbolic.interop.sympy import SympyToPymbolicMapper + from sumpy.symbolic import SympyToPymbolicMapper sympy_conv = SympyToPymbolicMapper() import loopy as lp @@ -383,4 +385,54 @@ class KernelCacheWrapper(object): return knl + +def my_syntactic_subs(expr, subst_dict): + # Workaround for differing substitution semantics between sympy and symengine. + # FIXME: This is a hack. + from sumpy.symbolic import Basic, Subs, Derivative, USE_SYMENGINE + + if not isinstance(expr, Basic): + return expr + + elif expr.is_Symbol: + return subst_dict.get(expr, expr) + + elif isinstance(expr, Subs): + new_point = tuple(my_syntactic_subs(p, subst_dict) for p in expr.point) + + import six + new_subst_dict = dict( + (var, subs) for var, subs in six.iteritems(subst_dict) + if var not in expr.variables) + + new_expr = my_syntactic_subs(expr.expr, new_subst_dict) + + if new_point != expr.point or new_expr != expr.expr: + return Subs(new_expr, expr.variables, new_point) + + return expr + + elif isinstance(expr, Derivative): + new_expr = my_syntactic_subs(expr.expr, subst_dict) + new_variables = my_syntactic_subs(expr.variables, subst_dict) + + if new_expr != expr.expr or any(new_var != var for new_var, var in + zip(new_variables, expr.variables)): + # FIXME in SymEngine + if USE_SYMENGINE: + return Derivative(new_expr, new_variables) + else: + return Derivative(new_expr, *new_variables) + + return expr + + else: + new_args = tuple(my_syntactic_subs(arg, subst_dict) for arg in expr.args) + + if any(new_arg != arg for arg, new_arg in zip(expr.args, new_args)): + return expr.func(*new_args) + + return expr + + # vim: fdm=marker diff --git a/sumpy/version.py b/sumpy/version.py index be874f65..f0497a83 100644 --- a/sumpy/version.py +++ b/sumpy/version.py @@ -25,4 +25,4 @@ VERSION = (2016, 1) VERSION_STATUS = "beta1" VERSION_TEXT = ".".join(str(x) for x in VERSION) + VERSION_STATUS -KERNEL_VERSION = 13 +KERNEL_VERSION = 14 diff --git a/test/test_codegen.py b/test/test_codegen.py index 59fe0340..c7b5e0bc 100644 --- a/test/test_codegen.py +++ b/test/test_codegen.py @@ -29,6 +29,38 @@ import logging logger = logging.getLogger(__name__) +def test_kill_trivial_assignments(): + from pymbolic import var + x, y, t0, t1, t2 = [var(s) for s in "x y t0 t1 t2".split()] + + assignments = ( + ("t0", 6), + ("t1", -t0), + ("t2", 6*x), + ("nt", x**y), + # users of trivial assignments + ("u0", t0 + 1), + ("u1", t1 + 1), + ("u2", t2 + 1), + ) + + from sumpy.codegen import kill_trivial_assignments + result = kill_trivial_assignments( + assignments, + retain_names=("u0", "u1", "u2")) + + from pymbolic.primitives import Sum + + def _s(*vals): + return Sum(vals) + + assert result == [ + ('nt', x**y), + ('u0', _s(6, 1)), + ('u1', _s(-6, 1)), + ('u2', _s(6*x, 1))] + + def test_symbolic_assignment_name_uniqueness(): # https://gitlab.tiker.net/inducer/sumpy/issues/13 from sumpy.assignment_collection import SymbolicAssignmentCollection diff --git a/test/test_cse.py b/test/test_cse.py index d2904176..4f957aaa 100644 --- a/test/test_cse.py +++ b/test/test_cse.py @@ -69,12 +69,13 @@ DAMAGE. import pytest import sys -from sympy import (Add, Pow, exp, sqrt, symbols, sympify, S, cos, sin, Eq, - Function, Tuple, CRootOf, IndexedBase, Idx, Piecewise) -from sympy.simplify.cse_opts import sub_pre, sub_post -from sympy.functions.special.hyper import meijerg -from sympy.simplify import cse_opts +from sumpy.symbolic import ( + Add, Pow, exp, sqrt, symbols, sympify, cos, sin, Function, USE_SYMENGINE) +if not USE_SYMENGINE: + from sympy.simplify.cse_opts import sub_pre, sub_post + from sympy.functions.special.hyper import meijerg + from sympy.simplify import cse_opts from sumpy.cse import ( cse, preprocess_for_cse, postprocess_for_cse) @@ -83,6 +84,9 @@ from sumpy.cse import ( w, x, y, z = symbols('w,x,y,z') x0, x1, x2, x3, x4, x5, x6, x7, x8, x9, x10, x11, x12 = symbols('x:13') +sympyonly = ( + pytest.mark.skipif(USE_SYMENGINE, reason="uses a sympy-only feature")) + # Dummy "optimization" functions for testing. @@ -122,6 +126,7 @@ def test_cse_single(): assert reduced == [sqrt(x0) + x0**2] +@sympyonly def test_cse_not_possible(): # No substitution possible. e = Add(x, y) @@ -142,6 +147,7 @@ def test_nested_substitution(): assert reduced == [sqrt(x0) + x0**2] +@sympyonly def test_subtraction_opt(): # Make sure subtraction is optimized. e = (x - y)*(z - y) + exp((x - y)*(z - y)) @@ -177,12 +183,13 @@ def test_multiple_expressions(): rsubsts, _ = cse(reversed(l)) assert substs == rsubsts assert reduced == [x1, x1 + z, x0] - l = [(x - z)*(y - z), x - z, y - z] + f = Function("f") + l = [f(x - z, y - z), x - z, y - z] substs, reduced = cse(l) rsubsts, _ = cse(reversed(l)) assert substs == [(x0, -z), (x1, x + x0), (x2, x0 + y)] assert rsubsts == [(x0, -z), (x1, x0 + y), (x2, x + x0)] - assert reduced == [x1*x2, x1, x2] + assert reduced == [f(x1, x2), x1, x2] l = [w*y + w + x + y + z, w*x*y] assert cse(l) == ([(x0, w*y)], [w + x + x0 + y + z, x*x0]) assert cse([x + y, x + y + z]) == ([(x0, x + y)], [x0, z + x0]) @@ -195,26 +202,31 @@ def test_issue_4203(): assert cse(sin(x**x)/x**x) == ([(x0, x**x)], [sin(x0)/x0]) -def test_dont_cse_tuples(): - from sympy import Subs +def test_dont_cse_subs(): + from sumpy.symbolic import Subs f = Function("f") g = Function("g") name_val, (expr,) = cse( - Subs(f(x, y), (x, y), (0, 1)) - + Subs(g(x, y), (x, y), (0, 1))) + Subs(f(x, y), (x, y), (0, x + y)) + + Subs(g(x, y), (x, y), (0, x + y))) assert name_val == [] - assert expr == (Subs(f(x, y), (x, y), (0, 1)) - + Subs(g(x, y), (x, y), (0, 1))) + assert expr == Subs(f(x, y), (x, y), (0, x + y)) + \ + Subs(g(x, y), (x, y), (0, x + y)) - name_val, (expr,) = cse( - Subs(f(x, y), (x, y), (0, x + y)) - + Subs(g(x, y), (x, y), (0, x + y))) - assert name_val == [(x0, x + y)] - assert expr == Subs(f(x, y), (x, y), (0, x0)) + \ - Subs(g(x, y), (x, y), (0, x0)) +def test_dont_cse_derivative(): + from sumpy.symbolic import Derivative + f = Function("f") + + # FIXME + deriv = Derivative(f(x+y), (x,)) if USE_SYMENGINE else Derivative(f(x+y), x) + + name_val, (expr,) = cse(x + y + deriv) + + assert name_val == [] + assert expr == x + y + deriv def test_pow_invpow(): @@ -238,9 +250,11 @@ def test_pow_invpow(): ([(x0, x**(2*y))], [x0 + 1/x0]) +@sympyonly def test_issue_4499(): # previously, this gave 16 constants from sympy.abc import a, b + from sympy import Tuple, S B = Function('B') # noqa G = Function('G') # noqa t = Tuple(* @@ -256,7 +270,9 @@ def test_issue_4499(): assert len(c[0]) == 11 +@sympyonly def test_issue_6169(): + from sympy import CRootOf r = CRootOf(x**6 - 4*x**5 - 2, 1) assert cse(r) == ([], [r]) # and a check that the right thing is done with the new @@ -264,7 +280,9 @@ def test_issue_6169(): assert sub_post(sub_pre((-x - y)*z - x - y)) == -z*(x + y) - x - y +@sympyonly def test_cse_Indexed(): # noqa + from sympy import IndexedBase, Idx len_y = 5 y = IndexedBase('y', shape=(len_y,)) x = IndexedBase('x', shape=(len_y,)) @@ -277,7 +295,9 @@ def test_cse_Indexed(): # noqa assert len(replacements) > 0 +@sympyonly def test_Piecewise(): # noqa + from sympy import Piecewise, Eq f = Piecewise((-z + x*y, Eq(y, 0)), (-z - x*y, True)) ans = cse(f) actual_ans = ([(x0, -z), (x1, x*y)], @@ -290,7 +310,7 @@ def test_name_conflict(): z2 = x2 + x3 l = [cos(z1) + z1, cos(z2) + z2, x0 + x2] substs, reduced = cse(l) - assert [e.subs(reversed(substs)) for e in reduced] == l + assert [e.subs(dict(substs)) for e in reduced] == l def test_name_conflict_cust_symbols(): @@ -298,7 +318,7 @@ def test_name_conflict_cust_symbols(): z2 = x2 + x3 l = [cos(z1) + z1, cos(z2) + z2, x0 + x2] substs, reduced = cse(l, symbols("x:10")) - assert [e.subs(reversed(substs)) for e in reduced] == l + assert [e.subs(dict(substs)) for e in reduced] == l def test_symbols_exhausted_error(): @@ -308,6 +328,7 @@ def test_symbols_exhausted_error(): print(cse(l, symbols=sym)) +@sympyonly def test_issue_7840(): # daveknippers' example C393 = sympify( # noqa diff --git a/test/test_fmm.py b/test/test_fmm.py index 2e2bfea4..33d5ddad 100644 --- a/test/test_fmm.py +++ b/test/test_fmm.py @@ -63,6 +63,7 @@ else: (l2d_level_to_order_lookup, ()), ]) def test_level_to_order_lookup(ctx_getter, lookup_func, extra_args): + pytest.importorskip("pyfmmlib") ctx = ctx_getter() queue = cl.CommandQueue(ctx) -- GitLab From 11f1697650e2458cbd10362ffea2939cb5472a1e Mon Sep 17 00:00:00 2001 From: Matt Wala Date: Fri, 19 May 2017 15:32:58 -0500 Subject: [PATCH 2/5] LineTaylorLocalExpansion: Reinstate the old derivative taking behavior when using sympy, so code generation doesn't blow up. Also add a test the uses LineTaylorLocalExpansion. --- sumpy/expansion/local.py | 33 ++++++++++--- test/test_qbx.py | 104 +++++++++++++++++++++++++++++++++++++++ 2 files changed, 129 insertions(+), 8 deletions(-) create mode 100644 test/test_qbx.py diff --git a/sumpy/expansion/local.py b/sumpy/expansion/local.py index f9ec106a..02a47fed 100644 --- a/sumpy/expansion/local.py +++ b/sumpy/expansion/local.py @@ -68,16 +68,33 @@ class LineTaylorLocalExpansion(LocalExpansionBase): line_kernel = self.kernel.get_expression(avec_line) - from sumpy.tools import MiDerivativeTaker, my_syntactic_subs - deriv_taker = MiDerivativeTaker(line_kernel, (tau,)) + from sumpy.symbolic import USE_SYMENGINE + + if USE_SYMENGINE: + from sumpy.tools import MiDerivativeTaker, my_syntactic_subs + deriv_taker = MiDerivativeTaker(line_kernel, (tau,)) + + return [my_syntactic_subs( + self.kernel.postprocess_at_target( + self.kernel.postprocess_at_source( + deriv_taker.diff(i), + avec), bvec), + {tau: 0}) + for i in self.get_coefficient_identifiers()] + else: + # Workaround for sympy. The automatic distribution after + # single-variable diff makes the expressions very large + # (https://github.com/sympy/sympy/issues/4596), so avoid doing + # single variable diff. + # + # See also https://gitlab.tiker.net/inducer/pytential/merge_requests/12 - return [my_syntactic_subs( - self.kernel.postprocess_at_target( + return [self.kernel.postprocess_at_target( self.kernel.postprocess_at_source( - deriv_taker.diff(i), - avec), bvec), - {tau: 0}) - for i in self.get_coefficient_identifiers()] + line_kernel.diff("tau", i), avec), + bvec) + .subs("tau", 0) + for i in self.get_coefficient_identifiers()] def evaluate(self, coeffs, bvec): from pytools import factorial diff --git a/test/test_qbx.py b/test/test_qbx.py new file mode 100644 index 00000000..a6c2319b --- /dev/null +++ b/test/test_qbx.py @@ -0,0 +1,104 @@ +from __future__ import division, absolute_import, print_function + +__copyright__ = "Copyright (C) 2017 Matt Wala" + +__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 numpy as np +import sys + +import pyopencl as cl + +import logging +logger = logging.getLogger(__name__) + + +try: + import faulthandler +except ImportError: + pass +else: + faulthandler.enable() + + +def test_direct(ctx_getter): + # This evaluates a single layer potential on a circle. + logging.basicConfig(level=logging.INFO) + + ctx = ctx_getter() + queue = cl.CommandQueue(ctx) + + from sumpy.kernel import LaplaceKernel + lknl = LaplaceKernel(2) + + order = 12 + + from sumpy.qbx import LayerPotential + from sumpy.expansion.local import LineTaylorLocalExpansion + + lpot = LayerPotential(ctx, [LineTaylorLocalExpansion(lknl, order)]) + + mode_nr = 25 + + from pytools.convergence import EOCRecorder + + eocrec = EOCRecorder() + + for n in [200, 300, 400]: + t = np.linspace(0, 2 * np.pi, n, endpoint=False) + unit_circle = np.exp(1j * t) + unit_circle = np.array([unit_circle.real, unit_circle.imag]) + + sigma = np.cos(mode_nr * t) + eigval = 1/(2*mode_nr) + + result_ref = eigval * sigma + + h = 2 * np.pi / n + + targets = unit_circle + sources = unit_circle + + radius = 7 * h + centers = unit_circle * (1 - radius) + + strengths = (sigma * h,) + evt, (result_qbx,) = lpot(queue, targets, sources, centers, strengths) + + eocrec.add_data_point(h, np.max(np.abs(result_ref - result_qbx))) + + print(eocrec) + + slack = 1.5 + assert eocrec.order_estimate() > order - slack + + +# You can test individual routines by typing +# $ python test_kernels.py 'test_p2p(cl.create_some_context)' + +if __name__ == "__main__": + if len(sys.argv) > 1: + exec(sys.argv[1]) + else: + from py.test.cmdline import main + main([__file__]) + +# vim: fdm=marker -- GitLab From 6de49f501434144524c031e06c4afb3251f2ed54 Mon Sep 17 00:00:00 2001 From: Matt Wala Date: Sat, 20 May 2017 15:19:08 -0500 Subject: [PATCH 3/5] Fix test_qbx (also, see #2). --- test/test_qbx.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/test/test_qbx.py b/test/test_qbx.py index a6c2319b..aa86edf8 100644 --- a/test/test_qbx.py +++ b/test/test_qbx.py @@ -26,6 +26,8 @@ import numpy as np import sys import pyopencl as cl +from pyopencl.tools import ( # noqa + pytest_generate_tests_for_pyopencl as pytest_generate_tests) import logging logger = logging.getLogger(__name__) -- GitLab From d3f99ed01676bd9565dcb240df96201772f4b1af Mon Sep 17 00:00:00 2001 From: Matt Wala Date: Sat, 20 May 2017 15:21:45 -0500 Subject: [PATCH 4/5] Add a regression test for LineTaylorLocalExpansion. --- test/test_codegen.py | 27 +++++++++++++++++++++++++++ 1 file changed, 27 insertions(+) diff --git a/test/test_codegen.py b/test/test_codegen.py index c7b5e0bc..fb165598 100644 --- a/test/test_codegen.py +++ b/test/test_codegen.py @@ -78,6 +78,33 @@ def test_symbolic_assignment_name_uniqueness(): assert len(sac.assignments) == 3 +def test_line_taylor_coeff_growth(): + # Regression test for LineTaylorLocalExpansion. + # See https://gitlab.tiker.net/inducer/pytential/merge_requests/12 + from sumpy.kernel import LaplaceKernel + from sumpy.expansion.local import LineTaylorLocalExpansion + from sumpy.symbolic import make_sym_vector, SympyToPymbolicMapper + + import numpy as np + + order = 10 + expn = LineTaylorLocalExpansion(LaplaceKernel(2), order) + avec = make_sym_vector("a", 2) + bvec = make_sym_vector("b", 2) + coeffs = expn.coefficients_from_source(avec, bvec) + + sym2pymbolic = SympyToPymbolicMapper() + coeffs_pymbolic = [sym2pymbolic(c) for c in coeffs] + + from pymbolic.mapper.flop_counter import FlopCounter + flop_counter = FlopCounter() + counts = [flop_counter(c) for c in coeffs_pymbolic] + + indices = np.arange(1, order + 2) + max_order = 2 + assert np.polyfit(np.log(indices), np.log(counts), deg=1)[0] < max_order + + # You can test individual routines by typing # $ python test_fmm.py 'test_sumpy_fmm(cl.create_some_context)' -- GitLab From ff160e08986d7e79661e30777e77cd442983451c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andreas=20Kl=C3=B6ckner?= Date: Thu, 25 May 2017 21:36:45 -0400 Subject: [PATCH 5/5] Pin conda-forge pocl to 0.13 (0.14 is broken without LLVM 4.0) --- .test-conda-env-py3.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.test-conda-env-py3.yml b/.test-conda-env-py3.yml index f336c9e9..bba7fcb0 100644 --- a/.test-conda-env-py3.yml +++ b/.test-conda-env-py3.yml @@ -8,7 +8,7 @@ dependencies: - git - conda-forge::numpy - conda-forge::sympy -- pocl +- pocl=0.13 - islpy - pyopencl - python=3.5 -- GitLab