From 552391553883572469fc53ff86c69060aaa84b19 Mon Sep 17 00:00:00 2001 From: Matt Wala Date: Fri, 27 Jan 2017 16:43:35 -0600 Subject: [PATCH 01/33] [WIP] Initial SymEngine integration. --- sumpy/assignment_collection.py | 40 ++--------- sumpy/codegen.py | 28 +++++++- sumpy/cse.py | 51 ++++++++------ sumpy/e2e.py | 10 +-- sumpy/e2p.py | 10 +-- sumpy/expansion/__init__.py | 4 +- sumpy/expansion/local.py | 53 +++++++------- sumpy/expansion/multipole.py | 49 ++++++------- sumpy/kernel.py | 65 ++++++++--------- sumpy/p2e.py | 5 +- sumpy/p2p.py | 4 +- sumpy/qbx.py | 10 +-- sumpy/symbolic.py | 123 +++++++++++++++++++++++++-------- sumpy/tools.py | 4 +- test/test_codegen.py | 66 ++++++++++++++++++ test/test_cse.py | 32 ++++++--- 16 files changed, 350 insertions(+), 204 deletions(-) create mode 100644 test/test_codegen.py diff --git a/sumpy/assignment_collection.py b/sumpy/assignment_collection.py index 840b04da..8aa07d81 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__) @@ -63,7 +63,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 @@ -125,7 +125,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 @@ -146,7 +146,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*. @@ -170,7 +170,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 @@ -185,41 +185,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..59cc9071 100644 --- a/sumpy/codegen.py +++ b/sumpy/codegen.py @@ -38,7 +38,7 @@ from loopy.types import NumpyType from pytools import memoize_method -from pymbolic.interop.sympy import ( +from sumpy.symbolic import (USE_SYMENGINE, SympyToPymbolicMapper as SympyToPymbolicMapperBase) import logging @@ -58,6 +58,10 @@ 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 @@ -77,6 +81,24 @@ class SympyToPymbolicMapper(SympyToPymbolicMapperBase): expr.expr.subs(self.assignments)), tuple(v.name for v in expr.variables)) + def not_supported(self, expr): + if isinstance(expr, int): + return expr + elif getattr(expr, "is_Function", False): + if USE_SYMENGINE: + try: + func_name = expr.get_name() + except AttributeError: + func_name = type(expr).__name__ + else: + func_name = type(expr).__name__ + # 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) # }}} @@ -318,7 +340,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 +351,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: diff --git a/sumpy/cse.py b/sumpy/cse.py index 350a2c8d..def1c7c6 100644 --- a/sumpy/cse.py +++ b/sumpy/cse.py @@ -66,8 +66,7 @@ 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 from sympy.core.compatibility import iterable from sympy.utilities.iterables import numbered_symbols @@ -192,19 +191,16 @@ class FuncArgTracker(object): (k, v) for k, v in count_map.items() if v >= threshold) - def get_subset_candidates(self, argset, restrict_to_funcset=None): + def get_subset_candidates(self, argset, restrict_to_funcset=frozenset()): """ Return a set of functions each of which whose argument list contains - `argset`, optionally filtered only to contain functions in - `restrict_to_funcset`. + `argset.` """ iarg = iter(argset) indices = set( - fi for fi in self.arg_to_funcset[next(iarg)]) - - if restrict_to_funcset is not None: - indices &= restrict_to_funcset + fi for fi in self.arg_to_funcset[next(iarg)] + if fi in restrict_to_funcset) for arg in iarg: indices &= self.arg_to_funcset[arg] @@ -227,6 +223,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 @@ -284,9 +293,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) @@ -313,9 +321,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) @@ -361,7 +368,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 @@ -374,7 +381,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)) # }}} @@ -415,10 +422,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 +464,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: @@ -475,7 +482,7 @@ def tree_cse(exprs, symbols, opt_subs=None): expr = opt_subs[expr] new_args = [rebuild(arg) for arg in expr.args] - if new_args != expr.args: + if isinstance(expr, Unevaluated) or new_args != expr.args: new_expr = expr.func(*new_args) else: new_expr = expr @@ -539,7 +546,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..f93a8682 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 @@ -202,6 +202,8 @@ class E2EFromCSR(E2EBase): silenced_warnings="write_race(write_expn*)", default_offset=lp.auto) + print(loopy_knl) + loopy_knl = lp.fix_parameters(loopy_knl, dim=self.dim) for expn in [self.src_expansion, self.tgt_expansion]: diff --git a/sumpy/e2p.py b/sumpy/e2p.py index b2e0e585..118f1743 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,13 +74,13 @@ 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 = [ @@ -107,7 +107,7 @@ class E2PBase(KernelCacheWrapper): 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..e958c2cb 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,7 +61,7 @@ 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 + avec_line = avec + sym.Symbol("tau")*bvec line_kernel = self.kernel.get_expression(avec_line) @@ -117,8 +117,9 @@ class VolumeTaylorLocalExpansionBase(LocalExpansionBase): type(self).__name__, self.order)) + import time from sumpy.expansion.multipole import VolumeTaylorMultipoleExpansionBase - if isinstance(src_expansion, VolumeTaylorMultipoleExpansionBase): + if 1 and isinstance(src_expansion, VolumeTaylorMultipoleExpansionBase): # We know the general form of the multipole expansion is: # # coeff0 * diff(kernel, mi0) + coeff1 * diff(kernel, mi1) + ... @@ -140,7 +141,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 +198,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 4a3987b5..0121cf6b 100644 --- a/sumpy/kernel.py +++ b/sumpy/kernel.py @@ -27,6 +27,7 @@ from six.moves import range, zip import loopy as lp import numpy as np from pymbolic.mapper import IdentityMapper, CSECachingMapperMixin +import sumpy.symbolic as sym from sumpy.symbolic import pymbolic_real_norm_2 from pymbolic.primitives import make_sym_vector from pymbolic import var @@ -263,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) + import sumpy.symbolic as sym + expr = expr.subs(dict( + (sym.Symbol("d%d" % i), dist_vec_i) for i, dist_vec_i in enumerate(dist_vec) - ]) + )) return expr @@ -506,22 +508,26 @@ class StressletKernel(ExpressionKernel): if dim == 2: d = make_sym_vector("d", dim) + n = make_sym_vector(stresslet_vector_name, dim) r = pymbolic_real_norm_2(d) expr = ( - -var("log")(r)*(1 if icomp == jcomp else 0) - + - d[icomp]*d[jcomp]/r**2 + sum(n[axis]*d[axis] for axis in range(dim)) + * + d[icomp]*d[jcomp]/r**4 ) - scaling = 1/(4*var("pi")) + scaling = 1/(var("pi")) + elif dim == 3: d = make_sym_vector("d", dim) + n = make_sym_vector(stresslet_vector_name, dim) r = pymbolic_real_norm_2(d) expr = ( - (1/r)*(1 if icomp == jcomp else 0) - + - d[icomp]*d[jcomp]/r**3 + sum(n[axis]*d[axis] for axis in range(dim)) + * + d[icomp]*d[jcomp]/r**5 ) - scaling = -1/(8*var("pi")) + scaling = -3/(4*var("pi")) + elif dim is None: expr = None scaling = None @@ -558,28 +564,15 @@ class StressletKernel(ExpressionKernel): return [ KernelArgument( loopy_arg=lp.ValueArg(self.viscosity_mu_name, np.float64), - )] - - def get_source_args(self): - return [ + ), KernelArgument( loopy_arg=lp.GlobalArg(self.stresslet_vector_name, None, shape=(self.dim, "nsources"), - dim_tags="sep,C"))] - - def postprocess_at_source(self, expr, avec): - dimensions = len(avec) - assert dimensions == self.dim - - from sumpy.symbolic import make_sympy_vector - n = make_sympy_vector(self.stresslet_vector_name, dimensions) + dim_tags="sep,C")) + ] - # avec = center-src -> minus sign from chain rule - return sum(-n[axis]*expr.diff(avec[axis]) - for axis in range(dimensions)) - - def get_code_transformer(self): + def get_code_tranformer(self): from sumpy.codegen import VectorComponentRewriter vcr = VectorComponentRewriter([self.stresslet_vector_name]) from pymbolic.primitives import Variable @@ -746,12 +739,12 @@ class DirectionalTargetDerivative(DirectionalDerivative): dim = len(bvec) assert dim == self.dim - from sumpy.symbolic import make_sympy_vector - dir_vec = make_sympy_vector(self.dir_vec_name, dim) + from sumpy.symbolic import make_sym_vector + dir_vec = make_sym_vector(self.dir_vec_name, dim) # bvec = tgt-center - return sum(dir_vec[axis]*expr.diff(bvec[axis]) - for axis in range(dim)) + return sym.Add(*tuple(dir_vec[axis]*expr.diff(bvec[axis]) + for axis in range(dim))) mapper_method = "map_directional_target_derivative" @@ -777,12 +770,12 @@ class DirectionalSourceDerivative(DirectionalDerivative): dimensions = len(avec) assert dimensions == self.dim - from sumpy.symbolic import make_sympy_vector - dir_vec = make_sympy_vector(self.dir_vec_name, dimensions) + from sumpy.symbolic import make_sym_vector + dir_vec = make_sym_vector(self.dir_vec_name, dimensions) # avec = center-src -> minus sign from chain rule - return sum(-dir_vec[axis]*expr.diff(avec[axis]) - for axis in range(dimensions)) + return sym.Add(*tuple(-dir_vec[axis]*expr.diff(avec[axis]) + for axis in range(dimensions))) mapper_method = "map_directional_source_derivative" diff --git a/sumpy/p2e.py b/sumpy/p2e.py index 1b037dd4..4f687ac5 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() @@ -161,6 +161,7 @@ class P2EFromSingleBox(P2EBase): loopy_knl = self.expansion.prepare_loopy_kernel(loopy_knl) loopy_knl = lp.tag_inames(loopy_knl, "idim*:unr") + print(loopy_knl) return loopy_knl def get_optimized_kernel(self): diff --git a/sumpy/p2p.py b/sumpy/p2p.py index 3071a3b8..a8b92b91 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() diff --git a/sumpy/qbx.py b/sumpy/qbx.py index 0b3dc035..5012e1a6 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() diff --git a/sumpy/symbolic.py b/sumpy/symbolic.py index 86c2666d..565fa8c5 100644 --- a/sumpy/symbolic.py +++ b/sumpy/symbolic.py @@ -27,27 +27,92 @@ 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__) +# {{{ symbolic backend + +def _find_symbolic_backend(): + global USE_SYMENGINE + + try: + import symengine # noqa + symengine_found = True + except ImportError: + symengine_found = False + + ALLOWED_BACKENDS = ("sympy", "symengine") # noqa + + import os + backend = os.environ.get("SUMPY_FORCE_SYMBOLIC_BACKEND") + if backend is not None: + if backend not in ALLOWED_BACKENDS: + raise RuntimeError( + "SUMPY_FORCE_SYMBOLIC_BACKEND value is unrecognized: '%s' " + "(allowed values are %s)" % ( + backend, + ", ".join("'%s'" % val for val in ALLOWED_BACKENDS))) + + if backend == "symengine" and not symengine_found: + from warnings import warn + warn("SUMPY_FORCE_SYMBOLIC_BACKEND=symengine was specified, but " + "could not find symengine. Using sympy.", RuntimeWarning) + + USE_SYMENGINE = backend == "symengine" and symengine_found + else: + USE_SYMENGINE = symengine_found + + +_find_symbolic_backend() + +# }}} + +# 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 +Integer Matrix Subs I pi functions""".split() + +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) + +for _apifunc in SYMBOLIC_API: + globals()[_apifunc] = getattr(sym, _apifunc) + + +def _coeff_isneg(a): + if a.is_Mul: + a = a.args[0] + try: + return a.is_Number and a < 0 + except: + return False + + # {{{ trivial assignment elimination def make_one_step_subst(assignments): - unwanted_vars = set(sp.Symbol(name) for name, value in assignments) + unwanted_vars = set(sym.Symbol(name) for name, value in assignments) result = {} - for name, value in assignments: - while value.atoms() & unwanted_vars: + assignments = dict((sym.Symbol(name), value) for name, value in assignments) + for name, value in assignments.items(): + while value.atoms(sym.Symbol) & unwanted_vars: value = value.subs(assignments) - result[sp.Symbol(name)] = value + result[name] = value return result @@ -55,12 +120,12 @@ def make_one_step_subst(assignments): def is_assignment_nontrivial(name, value): if value.is_Number: return False - elif isinstance(value, sp.Symbol): + elif isinstance(value, sym.Symbol): return False - elif (isinstance(value, sp.Mul) + elif (isinstance(value, sym.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): + and sum(1 for arg in value.args if isinstance(arg, sym.Symbol)) == 1): # const*var: not good enough return False @@ -81,8 +146,10 @@ def kill_trivial_assignments(assignments, retain_names=set()): # un-substitute rejected assignments unsubst_rej = make_one_step_subst(rejected_assignments) - result = [(name, expr.xreplace(unsubst_rej)) - for name, expr in approved_assignments] + result = [] + for name, expr in approved_assignments: + r = expr.xreplace(unsubst_rej) + result.append((name, r)) logger.info("kill trivial assignments (plain): done") @@ -111,7 +178,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 +185,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 +209,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 +237,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 +246,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 +260,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 +271,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..7d0738a7 100644 --- a/sumpy/tools.py +++ b/sumpy/tools.py @@ -107,7 +107,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 +264,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 diff --git a/test/test_codegen.py b/test/test_codegen.py new file mode 100644 index 00000000..fcb1dd58 --- /dev/null +++ b/test/test_codegen.py @@ -0,0 +1,66 @@ +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 sys + +import logging +logger = logging.getLogger(__name__) + + +def test_kill_trivial_assignments(): + from sumpy.symbolic import kill_trivial_assignments, symbols, sympify + x, y, nt = symbols("x, y, nt") + t0, t1, t2 = symbols("t0:3") + u0, u1, u2 = symbols("u0:3") + + assignments = ( + ("t0", sympify(6)), + ("t1", -t0), + ("t2", 6*x), + ("nt", x**y), + # users of trivial assignments + ("u0", t0 + 1), + ("u1", t1 + 1), + ("u2", t2 + 1) + ) + + result = kill_trivial_assignments( + assignments, + retain_names=("u0", "u1", "u2")) + + assert result == [('nt', x**y), ('u0', 7), ('u1', -5), ('u2', 1 + 6*x)] + + +# You can test individual routines by typing +# $ python test_fmm.py 'test_sumpy_fmm(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 diff --git a/test/test_cse.py b/test/test_cse.py index d2904176..a13ed6ab 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)) @@ -195,8 +201,9 @@ def test_issue_4203(): assert cse(sin(x**x)/x**x) == ([(x0, x**x)], [sin(x0)/x0]) +@sympyonly def test_dont_cse_tuples(): - from sympy import Subs + from sumpy.symbolic import Subs f = Function("f") g = Function("g") @@ -238,9 +245,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 +265,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 +275,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 +290,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 +305,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 +313,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 +323,7 @@ def test_symbols_exhausted_error(): print(cse(l, symbols=sym)) +@sympyonly def test_issue_7840(): # daveknippers' example C393 = sympify( # noqa -- GitLab From 6ae214053785166ef98b03e0dc5b21a750bd0020 Mon Sep 17 00:00:00 2001 From: Matt Wala Date: Fri, 27 Jan 2017 16:46:01 -0600 Subject: [PATCH 02/33] Undo debugging changes to sumpy.expansion.local --- sumpy/expansion/local.py | 53 ++++++++++++++++++++-------------------- 1 file changed, 26 insertions(+), 27 deletions(-) diff --git a/sumpy/expansion/local.py b/sumpy/expansion/local.py index e958c2cb..3ba643cd 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 sumpy.symbolic as sym +import sympy as sp from sumpy.expansion import ( ExpansionBase, VolumeTaylorExpansion, LaplaceConformingVolumeTaylorExpansion, @@ -61,7 +61,7 @@ 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 + sym.Symbol("tau")*bvec + avec_line = avec + sp.Symbol("tau")*bvec line_kernel = self.kernel.get_expression(avec_line) @@ -117,9 +117,8 @@ class VolumeTaylorLocalExpansionBase(LocalExpansionBase): type(self).__name__, self.order)) - import time from sumpy.expansion.multipole import VolumeTaylorMultipoleExpansionBase - if 1 and isinstance(src_expansion, VolumeTaylorMultipoleExpansionBase): + if isinstance(src_expansion, VolumeTaylorMultipoleExpansionBase): # We know the general form of the multipole expansion is: # # coeff0 * diff(kernel, mi0) + coeff1 * diff(kernel, mi1) + ... @@ -141,7 +140,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(sym.Add(*local_result)) + result.append(sp.Add(*local_result)) else: from sumpy.tools import MiDerivativeTaker expr = src_expansion.evaluate(src_coeff_exprs, dvec) @@ -198,61 +197,61 @@ class H2DLocalExpansion(LocalExpansionBase): return list(range(-self.order, self.order+1)) def coefficients_from_source(self, avec, bvec): - from sumpy.symbolic import sym_real_norm_2 - hankel_1 = sym.Function("hankel_1") + from sumpy.symbolic import sympy_real_norm_2 + hankel_1 = sp.Function("hankel_1") - k = sym.Symbol(self.kernel.get_base_kernel().helmholtz_k_name) + k = sp.Symbol(self.kernel.get_base_kernel().helmholtz_k_name) # The coordinates are negated since avec points from source to center. - source_angle_rel_center = sym.atan2(-avec[1], -avec[0]) - avec_len = sym_real_norm_2(avec) + source_angle_rel_center = sp.atan2(-avec[1], -avec[0]) + avec_len = sympy_real_norm_2(avec) return [self.kernel.postprocess_at_source( hankel_1(l, k * avec_len) - * sym.exp(sym.I * l * source_angle_rel_center), avec) + * sp.exp(sp.I * l * source_angle_rel_center), avec) for l in self.get_coefficient_identifiers()] def evaluate(self, coeffs, bvec): - 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]) + 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]) - k = sym.Symbol(self.kernel.get_base_kernel().helmholtz_k_name) + k = sp.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) - * sym.exp(sym.I * l * -target_angle_rel_center), bvec) + * sp.exp(sp.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 sym_real_norm_2 + from sumpy.symbolic import sympy_real_norm_2 - k = sym.Symbol(self.kernel.get_base_kernel().helmholtz_k_name) + k = sp.Symbol(self.kernel.get_base_kernel().helmholtz_k_name) if isinstance(src_expansion, H2DLocalExpansion): - 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]) + 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]) 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) - * sym.exp(sym.I * (m - l) * -new_center_angle_rel_old_center) + * sp.exp(sp.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 = sym_real_norm_2(dvec) - hankel_1 = sym.Function("hankel_1") - new_center_angle_rel_old_center = sym.atan2(dvec[1], dvec[0]) + 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]) translated_coeffs = [] for l in self.get_coefficient_identifiers(): translated_coeffs.append( sum((-1) ** l * hankel_1(m + l, k * dvec_len) - * sym.exp(sym.I * (m + l) * new_center_angle_rel_old_center) + * sp.exp(sp.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 -- GitLab From 5010870c053f81c321a33991db5b4dc9106380fc Mon Sep 17 00:00:00 2001 From: Matt Wala Date: Fri, 27 Jan 2017 18:25:33 -0600 Subject: [PATCH 03/33] Revert "Undo debugging changes to sumpy.expansion.local" This reverts commit 6ae214053785166ef98b03e0dc5b21a750bd0020. --- sumpy/expansion/local.py | 53 ++++++++++++++++++++-------------------- 1 file changed, 27 insertions(+), 26 deletions(-) diff --git a/sumpy/expansion/local.py b/sumpy/expansion/local.py index 3ba643cd..e958c2cb 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,7 +61,7 @@ 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 + avec_line = avec + sym.Symbol("tau")*bvec line_kernel = self.kernel.get_expression(avec_line) @@ -117,8 +117,9 @@ class VolumeTaylorLocalExpansionBase(LocalExpansionBase): type(self).__name__, self.order)) + import time from sumpy.expansion.multipole import VolumeTaylorMultipoleExpansionBase - if isinstance(src_expansion, VolumeTaylorMultipoleExpansionBase): + if 1 and isinstance(src_expansion, VolumeTaylorMultipoleExpansionBase): # We know the general form of the multipole expansion is: # # coeff0 * diff(kernel, mi0) + coeff1 * diff(kernel, mi1) + ... @@ -140,7 +141,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 +198,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 -- GitLab From b29a252b5e72f40e4d683b74f321e57e79b37774 Mon Sep 17 00:00:00 2001 From: Matt Wala Date: Fri, 27 Jan 2017 18:32:21 -0600 Subject: [PATCH 04/33] Actually fix local.py --- sumpy/expansion/local.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/sumpy/expansion/local.py b/sumpy/expansion/local.py index e958c2cb..bdcd9043 100644 --- a/sumpy/expansion/local.py +++ b/sumpy/expansion/local.py @@ -117,9 +117,8 @@ class VolumeTaylorLocalExpansionBase(LocalExpansionBase): type(self).__name__, self.order)) - import time from sumpy.expansion.multipole import VolumeTaylorMultipoleExpansionBase - if 1 and isinstance(src_expansion, VolumeTaylorMultipoleExpansionBase): + if isinstance(src_expansion, VolumeTaylorMultipoleExpansionBase): # We know the general form of the multipole expansion is: # # coeff0 * diff(kernel, mi0) + coeff1 * diff(kernel, mi1) + ... -- GitLab From 8621651e46788250a458beb13301cc936d948330 Mon Sep 17 00:00:00 2001 From: Matt Wala Date: Sat, 28 Jan 2017 17:06:10 -0600 Subject: [PATCH 05/33] Move kill_trivial_assignments() over to pymbolic. When doing substitution, SymEngine seems to collapse Subs nodes. This was leading to problems for sumpy, which relies on having Subs nodes in the expressions after killing trivial assignments. This change moves the job to pymbolic, which appears to preserve Subs nodes. --- sumpy/codegen.py | 143 +++++++++++++++++++++++++++++++++++++++---- sumpy/e2e.py | 9 +-- sumpy/e2p.py | 11 ++-- sumpy/p2e.py | 9 +-- sumpy/p2p.py | 1 + sumpy/qbx.py | 10 +-- sumpy/symbolic.py | 59 +----------------- sumpy/tools.py | 1 + test/test_codegen.py | 22 ++++--- 9 files changed, 161 insertions(+), 104 deletions(-) diff --git a/sumpy/codegen.py b/sumpy/codegen.py index 59cc9071..32492edc 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 sumpy.symbolic import (USE_SYMENGINE, - SympyToPymbolicMapper as SympyToPymbolicMapperBase) +from sumpy.symbolic import (SympyToPymbolicMapper as SympyToPymbolicMapperBase) import logging logger = logging.getLogger(__name__) @@ -64,7 +64,8 @@ _SPECIAL_FUNCTION_NAMES = frozenset(dir(sym.functions)) class SympyToPymbolicMapper(SympyToPymbolicMapperBase): def __init__(self, assignments): - self.assignments = assignments + self.assignments = dict( + (sym.Symbol(name), value) for name, value in assignments) self.derivative_cse_names = set() def map_Derivative(self, expr): # noqa @@ -85,13 +86,7 @@ class SympyToPymbolicMapper(SympyToPymbolicMapperBase): if isinstance(expr, int): return expr elif getattr(expr, "is_Function", False): - if USE_SYMENGINE: - try: - func_name = expr.get_name() - except AttributeError: - func_name = type(expr).__name__ - else: - func_name = type(expr).__name__ + 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() @@ -103,6 +98,130 @@ class SympyToPymbolicMapper(SympyToPymbolicMapperBase): # }}} +# {{{ 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) + + # }}} + + # {{{ make substitution + + from pymbolic import substitute + + result = {} + + for name in toposort: + value = assignments[name] + value = substitute(value, result) + + result[name] = value + + # }}} + + # {{{ simplify substitution + + used_names = set(dep + for value in six.itervalues(result) + for dep in get_dependencies(value)) + + used_name_to_var = dict( + (used_name, prim.Variable(used_name)) for used_name in used_names) + + from pymbolic import evaluate + from functools import partial + simplify = partial(evaluate, context=used_name_to_var) + + for name, value in six.iteritems(result): + 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 + +# }}} + + # {{{ bessel handling BESSEL_PREAMBLE = """//CL// @@ -560,7 +679,7 @@ 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) @@ -572,6 +691,8 @@ def to_loopy_insns(assignments, vector_names=set(), pymbolic_expr_maps=[], 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/e2e.py b/sumpy/e2e.py index f93a8682..b8091140 100644 --- a/sumpy/e2e.py +++ b/sumpy/e2e.py @@ -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 118f1743..6a77a5cc 100644 --- a/sumpy/e2p.py +++ b/sumpy/e2p.py @@ -83,6 +83,7 @@ class E2PBase(KernelCacheWrapper): 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,16 +92,12 @@ 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 ) diff --git a/sumpy/p2e.py b/sumpy/p2e.py index 4f687ac5..e85213d5 100644 --- a/sumpy/p2e.py +++ b/sumpy/p2e.py @@ -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 a8b92b91..39a90205 100644 --- a/sumpy/p2p.py +++ b/sumpy/p2p.py @@ -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 5012e1a6..8892739e 100644 --- a/sumpy/qbx.py +++ b/sumpy/qbx.py @@ -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 565fa8c5..0fcda1fb 100644 --- a/sumpy/symbolic.py +++ b/sumpy/symbolic.py @@ -76,7 +76,7 @@ _find_symbolic_backend() # 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 -Integer Matrix Subs I pi functions""".split() +Derivative Integer Matrix Subs I pi functions""".split() if USE_SYMENGINE: from symengine import sympy_compat as sym @@ -101,63 +101,6 @@ def _coeff_isneg(a): return False -# {{{ trivial assignment elimination - -def make_one_step_subst(assignments): - unwanted_vars = set(sym.Symbol(name) for name, value in assignments) - - result = {} - assignments = dict((sym.Symbol(name), value) for name, value in assignments) - for name, value in assignments.items(): - while value.atoms(sym.Symbol) & unwanted_vars: - value = value.subs(assignments) - - result[name] = value - - return result - - -def is_assignment_nontrivial(name, value): - if value.is_Number: - return False - elif isinstance(value, sym.Symbol): - return False - elif (isinstance(value, sym.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, sym.Symbol)) == 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 = [] - for name, expr in approved_assignments: - r = expr.xreplace(unsubst_rej) - result.append((name, r)) - - logger.info("kill trivial assignments (plain): done") - - return result - -# }}} - - # {{{ debugging of sympy CSE via Maxima class _DerivativeKiller(IdentityMapperBase): diff --git a/sumpy/tools.py b/sumpy/tools.py index 7d0738a7..4c952e06 100644 --- a/sumpy/tools.py +++ b/sumpy/tools.py @@ -74,6 +74,7 @@ class MiDerivativeTaker(object): for next_deriv, next_mi in self.get_derivative_taking_sequence( current_mi, mi): expr = expr.diff(next_deriv) + print("taking diff", expr) self.cache_by_mi[next_mi] = expr return expr diff --git a/test/test_codegen.py b/test/test_codegen.py index fcb1dd58..376a8fcf 100644 --- a/test/test_codegen.py +++ b/test/test_codegen.py @@ -30,27 +30,35 @@ logger = logging.getLogger(__name__) def test_kill_trivial_assignments(): - from sumpy.symbolic import kill_trivial_assignments, symbols, sympify - x, y, nt = symbols("x, y, nt") - t0, t1, t2 = symbols("t0:3") - u0, u1, u2 = symbols("u0:3") + from pymbolic import var + x, y, t0, t1, t2 = [var(s) for s in "x y t0 t1 t2".split()] assignments = ( - ("t0", sympify(6)), + ("t0", 6), ("t1", -t0), ("t2", 6*x), ("nt", x**y), # users of trivial assignments ("u0", t0 + 1), ("u1", t1 + 1), - ("u2", t2 + 1) + ("u2", t2 + 1), ) + from sumpy.codegen import kill_trivial_assignments result = kill_trivial_assignments( assignments, retain_names=("u0", "u1", "u2")) - assert result == [('nt', x**y), ('u0', 7), ('u1', -5), ('u2', 1 + 6*x)] + 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))] # You can test individual routines by typing -- GitLab From 4ae80d1c01c0a55e50901aa564233040fc971f1d Mon Sep 17 00:00:00 2001 From: Matt Wala Date: Sun, 29 Jan 2017 17:16:06 -0600 Subject: [PATCH 06/33] Fold simplification into substitution. --- sumpy/codegen.py | 26 ++++++++------------------ 1 file changed, 8 insertions(+), 18 deletions(-) diff --git a/sumpy/codegen.py b/sumpy/codegen.py index 32492edc..01385642 100644 --- a/sumpy/codegen.py +++ b/sumpy/codegen.py @@ -147,29 +147,19 @@ def make_one_step_subst(assignments): 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] = value - - # }}} - - # {{{ simplify substitution - - used_names = set(dep - for value in six.itervalues(result) - for dep in get_dependencies(value)) - - used_name_to_var = dict( - (used_name, prim.Variable(used_name)) for used_name in used_names) - - from pymbolic import evaluate - from functools import partial - simplify = partial(evaluate, context=used_name_to_var) - - for name, value in six.iteritems(result): result[name] = simplify(value) # }}} -- GitLab From abb38ae1e21255620177a3c41c649c3f5b176f4f Mon Sep 17 00:00:00 2001 From: Matt Wala Date: Sun, 29 Jan 2017 17:52:08 -0600 Subject: [PATCH 07/33] Bump kernel version. --- sumpy/version.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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 -- GitLab From 8df873e46f6ce43734e75e9c5e3b9029e8eab324 Mon Sep 17 00:00:00 2001 From: Matt Wala Date: Sun, 29 Jan 2017 17:59:41 -0600 Subject: [PATCH 08/33] Remove debugging print statements that crept in. --- sumpy/e2e.py | 2 -- sumpy/p2e.py | 1 - sumpy/tools.py | 1 - 3 files changed, 4 deletions(-) diff --git a/sumpy/e2e.py b/sumpy/e2e.py index b8091140..66e5e975 100644 --- a/sumpy/e2e.py +++ b/sumpy/e2e.py @@ -197,8 +197,6 @@ class E2EFromCSR(E2EBase): silenced_warnings="write_race(write_expn*)", default_offset=lp.auto) - print(loopy_knl) - loopy_knl = lp.fix_parameters(loopy_knl, dim=self.dim) for expn in [self.src_expansion, self.tgt_expansion]: diff --git a/sumpy/p2e.py b/sumpy/p2e.py index e85213d5..6a8fcfa9 100644 --- a/sumpy/p2e.py +++ b/sumpy/p2e.py @@ -156,7 +156,6 @@ class P2EFromSingleBox(P2EBase): loopy_knl = self.expansion.prepare_loopy_kernel(loopy_knl) loopy_knl = lp.tag_inames(loopy_knl, "idim*:unr") - print(loopy_knl) return loopy_knl def get_optimized_kernel(self): diff --git a/sumpy/tools.py b/sumpy/tools.py index 4c952e06..7d0738a7 100644 --- a/sumpy/tools.py +++ b/sumpy/tools.py @@ -74,7 +74,6 @@ class MiDerivativeTaker(object): for next_deriv, next_mi in self.get_derivative_taking_sequence( current_mi, mi): expr = expr.diff(next_deriv) - print("taking diff", expr) self.cache_by_mi[next_mi] = expr return expr -- GitLab From 185605229d24e5d0d3b080c9b91ec9acce0a0226 Mon Sep 17 00:00:00 2001 From: Matt Wala Date: Tue, 31 Jan 2017 20:13:44 -0600 Subject: [PATCH 09/33] CSE: Fix unnecessary rebuilds due to tuple != list. --- sumpy/cse.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/sumpy/cse.py b/sumpy/cse.py index def1c7c6..a073983a 100644 --- a/sumpy/cse.py +++ b/sumpy/cse.py @@ -481,7 +481,7 @@ 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] + 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) else: -- GitLab From bc97b5352548bdfb45764336d702925974b3b067 Mon Sep 17 00:00:00 2001 From: Matt Wala Date: Wed, 1 Feb 2017 23:09:31 -0600 Subject: [PATCH 10/33] test_cse: Fix failure due to expression ordering difference in SymEngine. --- test/test_cse.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/test/test_cse.py b/test/test_cse.py index a13ed6ab..c9f0264e 100644 --- a/test/test_cse.py +++ b/test/test_cse.py @@ -183,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]) -- GitLab From 661f58514353a360d18a4548d21331a7f3bab568 Mon Sep 17 00:00:00 2001 From: Matt Wala Date: Wed, 1 Feb 2017 23:10:43 -0600 Subject: [PATCH 11/33] CSE: Don't pull things out of derivatives, Subs (closes #14). --- sumpy/codegen.py | 22 ---------------------- sumpy/cse.py | 25 ++++++++++++++++++------- test/test_cse.py | 30 +++++++++++++++++------------- 3 files changed, 35 insertions(+), 42 deletions(-) diff --git a/sumpy/codegen.py b/sumpy/codegen.py index 01385642..8eef9892 100644 --- a/sumpy/codegen.py +++ b/sumpy/codegen.py @@ -63,24 +63,6 @@ _SPECIAL_FUNCTION_NAMES = frozenset(dir(sym.functions)) class SympyToPymbolicMapper(SympyToPymbolicMapperBase): - def __init__(self, assignments): - self.assignments = dict( - (sym.Symbol(name), value) for name, value in 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. - - 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)) def not_supported(self, expr): if isinstance(expr, int): @@ -676,10 +658,6 @@ def to_loopy_insns(assignments, vector_names=set(), pymbolic_expr_maps=[], # convert from sympy sympy_conv = SympyToPymbolicMapper(assignments) 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) diff --git a/sumpy/cse.py b/sumpy/cse.py index a073983a..6a0270e7 100644 --- a/sumpy/cse.py +++ b/sumpy/cse.py @@ -66,7 +66,8 @@ DAMAGE. # }}} -from sumpy.symbolic import Basic, Mul, Add, Pow, Symbol, _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 @@ -81,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): @@ -352,6 +357,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) @@ -443,7 +451,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) @@ -481,11 +492,11 @@ def tree_cse(exprs, symbols, opt_subs=None): if expr in opt_subs: expr = opt_subs[expr] - 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) - 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: diff --git a/test/test_cse.py b/test/test_cse.py index c9f0264e..4f957aaa 100644 --- a/test/test_cse.py +++ b/test/test_cse.py @@ -189,7 +189,7 @@ def test_multiple_expressions(): 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 == [f(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]) @@ -202,27 +202,31 @@ def test_issue_4203(): assert cse(sin(x**x)/x**x) == ([(x0, x**x)], [sin(x0)/x0]) -@sympyonly -def test_dont_cse_tuples(): +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(): -- GitLab From 415ae2a6638f6f72c959f406c40f802f942a12da Mon Sep 17 00:00:00 2001 From: Matt Wala Date: Wed, 1 Feb 2017 23:13:21 -0600 Subject: [PATCH 12/33] Fix SympyToPymbolicMapper api change. --- sumpy/codegen.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/sumpy/codegen.py b/sumpy/codegen.py index 8eef9892..99e4f47d 100644 --- a/sumpy/codegen.py +++ b/sumpy/codegen.py @@ -656,7 +656,7 @@ def to_loopy_insns(assignments, vector_names=set(), pymbolic_expr_maps=[], assignments = list(assignments) # convert from sympy - sympy_conv = SympyToPymbolicMapper(assignments) + sympy_conv = SympyToPymbolicMapper() assignments = [(name, sympy_conv(expr)) for name, expr in assignments] assignments = kill_trivial_assignments(assignments, retain_names) -- GitLab From ffc1575f36a86374eb83902bc83fe793863acaf6 Mon Sep 17 00:00:00 2001 From: Matt Wala Date: Thu, 2 Feb 2017 11:54:18 -0600 Subject: [PATCH 13/33] Update _coeff_isneg to use is_negative. --- sumpy/symbolic.py | 5 +---- 1 file changed, 1 insertion(+), 4 deletions(-) diff --git a/sumpy/symbolic.py b/sumpy/symbolic.py index 0fcda1fb..5cbd3df2 100644 --- a/sumpy/symbolic.py +++ b/sumpy/symbolic.py @@ -95,10 +95,7 @@ for _apifunc in SYMBOLIC_API: def _coeff_isneg(a): if a.is_Mul: a = a.args[0] - try: - return a.is_Number and a < 0 - except: - return False + return a.is_Number and a.is_negative # {{{ debugging of sympy CSE via Maxima -- GitLab From 9348fa80a49b376f9a779dc993ada4df13e731f4 Mon Sep 17 00:00:00 2001 From: Matt Wala Date: Thu, 2 Feb 2017 13:05:56 -0600 Subject: [PATCH 14/33] Add an assert to MiDerivativeTaker to ensure we haven't fallen back to sympy. --- sumpy/tools.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/sumpy/tools.py b/sumpy/tools.py index 7d0738a7..51cc52be 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} -- GitLab From 1dfb00960346680aa72907f60bce62fc52eb00d0 Mon Sep 17 00:00:00 2001 From: Matt Wala Date: Wed, 25 Jan 2017 18:01:30 -0600 Subject: [PATCH 15/33] CSE: get_subset_candidates() shouldn't default to an empty whitelist. --- sumpy/cse.py | 11 +++++++---- 1 file changed, 7 insertions(+), 4 deletions(-) diff --git a/sumpy/cse.py b/sumpy/cse.py index 6a0270e7..3318fba0 100644 --- a/sumpy/cse.py +++ b/sumpy/cse.py @@ -196,16 +196,19 @@ class FuncArgTracker(object): (k, v) for k, v in count_map.items() if v >= threshold) - def get_subset_candidates(self, argset, restrict_to_funcset=frozenset()): + def get_subset_candidates(self, argset, restrict_to_funcset=None): """ Return a set of functions each of which whose argument list contains - `argset.` + `argset`, optionally filtered only to contain functions in + `restrict_to_funcset`. """ iarg = iter(argset) indices = set( - fi for fi in self.arg_to_funcset[next(iarg)] - if fi in restrict_to_funcset) + fi for fi in self.arg_to_funcset[next(iarg)]) + + if restrict_to_funcset is not None: + indices &= restrict_to_funcset for arg in iarg: indices &= self.arg_to_funcset[arg] -- GitLab From 623e1fd16fb8310d15dd7668d537eab901402f0a Mon Sep 17 00:00:00 2001 From: Matt Wala Date: Thu, 2 Feb 2017 14:23:23 -0600 Subject: [PATCH 16/33] Document symbolic backends. --- doc/misc.rst | 13 +++++++++++++ 1 file changed, 13 insertions(+) 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 ==================== -- GitLab From ef66a1491f4f23a2dbf7d7648451e136ab4b50d1 Mon Sep 17 00:00:00 2001 From: Matt Wala Date: Fri, 17 Feb 2017 01:33:40 -0600 Subject: [PATCH 17/33] Fix subs in symengine. --- sumpy/kernel.py | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/sumpy/kernel.py b/sumpy/kernel.py index e29ac516..60f53011 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 -- GitLab From 45540e2ed41ba674fa5047c78163096deec99b05 Mon Sep 17 00:00:00 2001 From: Matt Wala Date: Fri, 17 Feb 2017 01:33:57 -0600 Subject: [PATCH 18/33] Skip pyfmmlib requiring tests. --- test/test_fmm.py | 1 + 1 file changed, 1 insertion(+) 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 d456ae49343067ba57fd0209eed4e8a6383a1279 Mon Sep 17 00:00:00 2001 From: Matt Wala Date: Fri, 17 Feb 2017 01:38:30 -0600 Subject: [PATCH 19/33] Initial attempt at CI infrastructure for SymEngine. --- .gitlab-ci.yml | 12 ++++++++++++ .test-py3-requirements.txt | 3 +++ .test-py3.yml | 16 ++++++++++++++++ 3 files changed, 31 insertions(+) create mode 100644 .test-py3-requirements.txt create mode 100644 .test-py3.yml diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index 6ee07ae4..675732e5 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-py3.yml + - REQUIREMENTS_TXT=.test-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-py3-requirements.txt b/.test-py3-requirements.txt new file mode 100644 index 00000000..01f08662 --- /dev/null +++ b/.test-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-py3.yml b/.test-py3.yml new file mode 100644 index 00000000..b1de4606 --- /dev/null +++ b/.test-py3.yml @@ -0,0 +1,16 @@ +name: 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 -- GitLab From b0e7f8765b5804afa621e34a5fd94212a77cd99c Mon Sep 17 00:00:00 2001 From: Matt Wala Date: Fri, 17 Feb 2017 01:49:33 -0600 Subject: [PATCH 20/33] [ci skip] Add py2.7 symengine tests as well. --- .gitlab-ci.yml | 12 ++++++++++++ .test-py2-requirements.txt | 3 +++ .test-py2.yml | 16 ++++++++++++++++ 3 files changed, 31 insertions(+) create mode 100644 .test-py2-requirements.txt create mode 100644 .test-py2.yml diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index 675732e5..913eb8d0 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -66,6 +66,18 @@ Python 3.6 POCL: except: - tags +Python 2.7 Conda: + script: + - export SUMPY_NO_CACHE=1 + - CONDA_ENVIRONMENT=.test-py2.yml + - REQUIREMENTS_TXT=.test-py2-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 + Python 3.5 Conda: script: - export SUMPY_NO_CACHE=1 diff --git a/.test-py2-requirements.txt b/.test-py2-requirements.txt new file mode 100644 index 00000000..01f08662 --- /dev/null +++ b/.test-py2-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-py2.yml b/.test-py2.yml new file mode 100644 index 00000000..1725a66a --- /dev/null +++ b/.test-py2.yml @@ -0,0 +1,16 @@ +name: py2 +channels: +- inducer +- symengine/label/dev +- conda-forge +- defaults +dependencies: +- git +- conda-forge::numpy +- conda-forge::sympy +- pocl +- islpy +- pyopencl +- python=2.7 +- python-symengine=0.2.0.53.g83912b7=py35_1 +# things not in here: loopy boxtree pymbolic pyfmmlib -- GitLab From d907320485f9be5ac23f189db6bda187a68d19ea Mon Sep 17 00:00:00 2001 From: Matt Wala Date: Fri, 17 Feb 2017 01:58:46 -0600 Subject: [PATCH 21/33] [ci skip] make_sympy_vector fixes. --- sumpy/kernel.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/sumpy/kernel.py b/sumpy/kernel.py index 60f53011..ea0ea02b 100644 --- a/sumpy/kernel.py +++ b/sumpy/kernel.py @@ -716,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 @@ -747,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 -- GitLab From d46ff9301480a102d2fb7f26b4d64fb6a03fde78 Mon Sep 17 00:00:00 2001 From: Matt Wala Date: Fri, 17 Feb 2017 12:51:49 -0600 Subject: [PATCH 22/33] Pull magic string out. --- sumpy/symbolic.py | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/sumpy/symbolic.py b/sumpy/symbolic.py index 5cbd3df2..85dbb278 100644 --- a/sumpy/symbolic.py +++ b/sumpy/symbolic.py @@ -47,21 +47,23 @@ def _find_symbolic_backend(): symengine_found = False ALLOWED_BACKENDS = ("sympy", "symengine") # noqa + BACKEND_ENV_VAR = "SUMPY_FORCE_SYMBOLIC_BACKEND" # noqa import os - backend = os.environ.get("SUMPY_FORCE_SYMBOLIC_BACKEND") + backend = os.environ.get(BACKEND_ENV_VAR) if backend is not None: if backend not in ALLOWED_BACKENDS: raise RuntimeError( - "SUMPY_FORCE_SYMBOLIC_BACKEND value is unrecognized: '%s' " + "%s value is unrecognized: '%s' " "(allowed values are %s)" % ( + BACKEND_ENV_VAR, backend, ", ".join("'%s'" % val for val in ALLOWED_BACKENDS))) if backend == "symengine" and not symengine_found: from warnings import warn - warn("SUMPY_FORCE_SYMBOLIC_BACKEND=symengine was specified, but " - "could not find symengine. Using sympy.", RuntimeWarning) + warn("%s=symengine was specified, but could not find symengine. " + "Using sympy." % BACKEND_ENV_VAR, RuntimeWarning) USE_SYMENGINE = backend == "symengine" and symengine_found else: -- GitLab From ade468f580b0902690955bb3f3744bc5aa287133 Mon Sep 17 00:00:00 2001 From: Matt Wala Date: Sun, 19 Feb 2017 18:14:52 -0600 Subject: [PATCH 23/33] Test requirements: specify pymbolic version. --- .test-py2-requirements.txt | 3 ++- .test-py3-requirements.txt | 3 ++- 2 files changed, 4 insertions(+), 2 deletions(-) diff --git a/.test-py2-requirements.txt b/.test-py2-requirements.txt index 01f08662..cbd5ec7e 100644 --- a/.test-py2-requirements.txt +++ b/.test-py2-requirements.txt @@ -1,3 +1,4 @@ git+https://github.com/inducer/boxtree -git+https://github.com/inducer/pymbolic git+https://github.com/inducer/loopy +# The pymbolic version needs SymEngine support +git+https://github.com/inducer/pymbolic@081714baf4476b42ca75a1324377914dd6356965 diff --git a/.test-py3-requirements.txt b/.test-py3-requirements.txt index 01f08662..cbd5ec7e 100644 --- a/.test-py3-requirements.txt +++ b/.test-py3-requirements.txt @@ -1,3 +1,4 @@ git+https://github.com/inducer/boxtree -git+https://github.com/inducer/pymbolic git+https://github.com/inducer/loopy +# The pymbolic version needs SymEngine support +git+https://github.com/inducer/pymbolic@081714baf4476b42ca75a1324377914dd6356965 -- GitLab From 41f37f9ccd164835d527a366133cf1b10c0af786 Mon Sep 17 00:00:00 2001 From: Matt Wala Date: Sun, 19 Feb 2017 18:25:32 -0600 Subject: [PATCH 24/33] Fix requirements.txt --- requirements.txt | 2 ++ 1 file changed, 2 insertions(+) diff --git a/requirements.txt b/requirements.txt index abe39e89..e28502fb 100644 --- a/requirements.txt +++ b/requirements.txt @@ -5,3 +5,5 @@ git+https://github.com/pyopencl/pyopencl git+https://github.com/inducer/boxtree git+https://github.com/inducer/loopy git+https://github.com/inducer/pyfmmlib +# The pymbolic version needs SymEngine support +git+https://github.com/inducer/pymbolic@081714baf4476b42ca75a1324377914dd6356965 -- GitLab From 2fa9851b9171b72def97600e5cc6d1881c013b76 Mon Sep 17 00:00:00 2001 From: Matt Wala Date: Sun, 19 Feb 2017 18:26:25 -0600 Subject: [PATCH 25/33] Revert "Test requirements: specify pymbolic version." This reverts commit ade468f580b0902690955bb3f3744bc5aa287133. --- .test-py2-requirements.txt | 3 +-- .test-py3-requirements.txt | 3 +-- 2 files changed, 2 insertions(+), 4 deletions(-) diff --git a/.test-py2-requirements.txt b/.test-py2-requirements.txt index cbd5ec7e..01f08662 100644 --- a/.test-py2-requirements.txt +++ b/.test-py2-requirements.txt @@ -1,4 +1,3 @@ git+https://github.com/inducer/boxtree +git+https://github.com/inducer/pymbolic git+https://github.com/inducer/loopy -# The pymbolic version needs SymEngine support -git+https://github.com/inducer/pymbolic@081714baf4476b42ca75a1324377914dd6356965 diff --git a/.test-py3-requirements.txt b/.test-py3-requirements.txt index cbd5ec7e..01f08662 100644 --- a/.test-py3-requirements.txt +++ b/.test-py3-requirements.txt @@ -1,4 +1,3 @@ git+https://github.com/inducer/boxtree +git+https://github.com/inducer/pymbolic git+https://github.com/inducer/loopy -# The pymbolic version needs SymEngine support -git+https://github.com/inducer/pymbolic@081714baf4476b42ca75a1324377914dd6356965 -- GitLab From 804989439c0666e20d6a771daa0962803f25096f Mon Sep 17 00:00:00 2001 From: Matt Wala Date: Sun, 19 Feb 2017 18:47:10 -0600 Subject: [PATCH 26/33] Revert "Fix requirements.txt" This reverts commit 41f37f9ccd164835d527a366133cf1b10c0af786. --- requirements.txt | 2 -- 1 file changed, 2 deletions(-) diff --git a/requirements.txt b/requirements.txt index e28502fb..abe39e89 100644 --- a/requirements.txt +++ b/requirements.txt @@ -5,5 +5,3 @@ git+https://github.com/pyopencl/pyopencl git+https://github.com/inducer/boxtree git+https://github.com/inducer/loopy git+https://github.com/inducer/pyfmmlib -# The pymbolic version needs SymEngine support -git+https://github.com/inducer/pymbolic@081714baf4476b42ca75a1324377914dd6356965 -- GitLab From 7b7f40c256bd66a7dd080b71231a528d089d2704 Mon Sep 17 00:00:00 2001 From: Matt Wala Date: Sun, 19 Feb 2017 18:51:34 -0600 Subject: [PATCH 27/33] Requirements.txt: specify pymbolic. --- requirements.txt | 1 + 1 file changed, 1 insertion(+) 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 -- GitLab From 4bb726a3ba0792ddde3ca915079f033527e99e18 Mon Sep 17 00:00:00 2001 From: Matt Wala Date: Sun, 19 Feb 2017 19:03:36 -0600 Subject: [PATCH 28/33] Remove py2.7 conda build, as it doesn't seem to work now since pocl in coda is py3.5+. --- .gitlab-ci.yml | 12 ------------ .test-py2-requirements.txt | 3 --- .test-py2.yml | 16 ---------------- 3 files changed, 31 deletions(-) delete mode 100644 .test-py2-requirements.txt delete mode 100644 .test-py2.yml diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index 913eb8d0..675732e5 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -66,18 +66,6 @@ Python 3.6 POCL: except: - tags -Python 2.7 Conda: - script: - - export SUMPY_NO_CACHE=1 - - CONDA_ENVIRONMENT=.test-py2.yml - - REQUIREMENTS_TXT=.test-py2-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 - Python 3.5 Conda: script: - export SUMPY_NO_CACHE=1 diff --git a/.test-py2-requirements.txt b/.test-py2-requirements.txt deleted file mode 100644 index 01f08662..00000000 --- a/.test-py2-requirements.txt +++ /dev/null @@ -1,3 +0,0 @@ -git+https://github.com/inducer/boxtree -git+https://github.com/inducer/pymbolic -git+https://github.com/inducer/loopy diff --git a/.test-py2.yml b/.test-py2.yml deleted file mode 100644 index 1725a66a..00000000 --- a/.test-py2.yml +++ /dev/null @@ -1,16 +0,0 @@ -name: py2 -channels: -- inducer -- symengine/label/dev -- conda-forge -- defaults -dependencies: -- git -- conda-forge::numpy -- conda-forge::sympy -- pocl -- islpy -- pyopencl -- python=2.7 -- python-symengine=0.2.0.53.g83912b7=py35_1 -# things not in here: loopy boxtree pymbolic pyfmmlib -- GitLab From 95b2d3554af942e2a31a552844932cc7f3aa74cc Mon Sep 17 00:00:00 2001 From: Matt Wala Date: Thu, 23 Feb 2017 18:34:07 -0600 Subject: [PATCH 29/33] Fix LineTaylorLocalExpansion to work with SymEngine (still needed: tests in sumpy, but that's a long-standing bug: see #2). --- sumpy/expansion/local.py | 25 ++++++++++++-------- sumpy/tools.py | 49 ++++++++++++++++++++++++++++++++++++++++ 2 files changed, 64 insertions(+), 10 deletions(-) diff --git a/sumpy/expansion/local.py b/sumpy/expansion/local.py index bdcd9043..f9ec106a 100644 --- a/sumpy/expansion/local.py +++ b/sumpy/expansion/local.py @@ -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 + sym.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())) # }}} diff --git a/sumpy/tools.py b/sumpy/tools.py index 51cc52be..af125566 100644 --- a/sumpy/tools.py +++ b/sumpy/tools.py @@ -385,4 +385,53 @@ 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) + else: + return expr + + # vim: fdm=marker -- GitLab From 8c2c3252f0b39f0ec498ca3e2488f0cc22e1e3ee Mon Sep 17 00:00:00 2001 From: Matt Wala Date: Thu, 23 Feb 2017 19:49:03 -0600 Subject: [PATCH 30/33] [ci skip] Code layout consistency. --- sumpy/tools.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/sumpy/tools.py b/sumpy/tools.py index af125566..f5014295 100644 --- a/sumpy/tools.py +++ b/sumpy/tools.py @@ -428,10 +428,11 @@ def my_syntactic_subs(expr, subst_dict): 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) - else: - return expr + + return expr # vim: fdm=marker -- GitLab From e5c92edb9d1f9e5573d76af519e6828ba36dd982 Mon Sep 17 00:00:00 2001 From: Matt Wala Date: Thu, 23 Feb 2017 23:17:53 -0600 Subject: [PATCH 31/33] Change name of conda testing environment. --- .gitlab-ci.yml | 4 ++-- ...3-requirements.txt => .test-conda-env-py3-requirements.txt | 0 .test-py3.yml => .test-conda-env-py3.yml | 2 +- 3 files changed, 3 insertions(+), 3 deletions(-) rename .test-py3-requirements.txt => .test-conda-env-py3-requirements.txt (100%) rename .test-py3.yml => .test-conda-env-py3.yml (91%) diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index 675732e5..63768430 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -69,8 +69,8 @@ Python 3.6 POCL: Python 3.5 Conda: script: - export SUMPY_NO_CACHE=1 - - CONDA_ENVIRONMENT=.test-py3.yml - - REQUIREMENTS_TXT=.test-py3-requirements.txt + - 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: diff --git a/.test-py3-requirements.txt b/.test-conda-env-py3-requirements.txt similarity index 100% rename from .test-py3-requirements.txt rename to .test-conda-env-py3-requirements.txt diff --git a/.test-py3.yml b/.test-conda-env-py3.yml similarity index 91% rename from .test-py3.yml rename to .test-conda-env-py3.yml index b1de4606..f336c9e9 100644 --- a/.test-py3.yml +++ b/.test-conda-env-py3.yml @@ -1,4 +1,4 @@ -name: py3 +name: test-conda-env-py3 channels: - inducer - symengine/label/dev -- GitLab From 3dd84f013ae2f819a4f7b0addc8292cfd264231b Mon Sep 17 00:00:00 2001 From: Matt Wala Date: Fri, 24 Feb 2017 02:23:12 -0600 Subject: [PATCH 32/33] Add an Apple build. --- .gitlab-ci.yml | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index 63768430..041d43dc 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -78,6 +78,18 @@ Python 3.5 Conda: except: - tags +Python 3.5 Apple: + 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: + - apple + except: + - tags + Documentation: script: - EXTRA_INSTALL="numpy mako" -- GitLab From c37e148b038fe8a54f56a590498c0d4b131b386e Mon Sep 17 00:00:00 2001 From: Matt Wala Date: Fri, 24 Feb 2017 16:12:58 -0600 Subject: [PATCH 33/33] Revert "Add an Apple build." This reverts commit 3dd84f013ae2f819a4f7b0addc8292cfd264231b. --- .gitlab-ci.yml | 12 ------------ 1 file changed, 12 deletions(-) diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index 041d43dc..63768430 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -78,18 +78,6 @@ Python 3.5 Conda: except: - tags -Python 3.5 Apple: - 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: - - apple - except: - - tags - Documentation: script: - EXTRA_INSTALL="numpy mako" -- GitLab