diff --git a/requirements.txt b/requirements.txt index a26c33e3b462457551693c89262966a504c1eb28..ff9ca871b7caebcaf92e6fda6c5a4f18b00b1ab2 100644 --- a/requirements.txt +++ b/requirements.txt @@ -8,4 +8,5 @@ git+https://github.com/inducer/islpy.git#egg=islpy git+https://github.com/inducer/pyopencl.git#egg=pyopencl git+https://github.com/inducer/boxtree.git#egg=boxtree git+https://github.com/inducer/loopy.git#egg=loopy +git+https://github.com/inducer/arraycontext.git#egg=arraycontext git+https://github.com/inducer/pyfmmlib.git#egg=pyfmmlib diff --git a/setup.py b/setup.py index f9f0be9bef525795257f584617ffd49893fdee46..19ca094fc5dfff531d95b574bb6e5898c2c3757e 100644 --- a/setup.py +++ b/setup.py @@ -101,6 +101,7 @@ setup( "pytools>=2021.1.1", "loopy>=2021.1", "boxtree>=2018.1", + "arraycontext", "pytest>=2.3", "pyrsistent>=0.16.0", "dataclasses>=0.7;python_version<='3.6'", diff --git a/sumpy/array_context.py b/sumpy/array_context.py new file mode 100644 index 0000000000000000000000000000000000000000..760f7d5d8362b2f171e26530c72be66316ab7fa3 --- /dev/null +++ b/sumpy/array_context.py @@ -0,0 +1,77 @@ +__copyright__ = "Copyright (C) 2022 Alexandru Fikl" + +__license__ = """ +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in +all copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN +THE SOFTWARE. +""" + +from boxtree.array_context import PyOpenCLArrayContext as PyOpenCLArrayContextBase +from arraycontext.pytest import ( + _PytestPyOpenCLArrayContextFactoryWithClass, + register_pytest_array_context_factory) + +__doc__ = """ +.. autoclass:: PyOpenCLArrayContext +""" + + +# {{{ PyOpenCLArrayContext + +class PyOpenCLArrayContext(PyOpenCLArrayContextBase): + def transform_loopy_program(self, t_unit): + default_ep = t_unit.default_entrypoint + options = default_ep.options + + if not (options.return_dict and options.no_numpy): + raise ValueError("Loopy kernel passed to call_loopy must " + "have return_dict and no_numpy options set. " + "Did you use arraycontext.make_loopy_program " + "to create this kernel?") + + return super().transform_loopy_program(t_unit) + +# }}} + + +# {{{ pytest + +def _acf(): + import pyopencl as cl + ctx = cl.create_some_context() + queue = cl.CommandQueue(ctx) + + return PyOpenCLArrayContext(queue, force_device_scalars=True) + + +class PytestPyOpenCLArrayContextFactory( + _PytestPyOpenCLArrayContextFactoryWithClass): + actx_class = PyOpenCLArrayContext + + def __call__(self): + # NOTE: prevent any cache explosions during testing! + from sympy.core.cache import clear_cache + clear_cache() + + return super().__call__() + + +register_pytest_array_context_factory( + "sumpy.pyopencl", + PytestPyOpenCLArrayContextFactory) + +# }}} diff --git a/sumpy/expansion/multipole.py b/sumpy/expansion/multipole.py index fd5a3b1d1ac3ca8fb1f05bb48b9d845a0bf4ff97..1924e6c7f0784f00a210d0818cb94b3682118c3b 100644 --- a/sumpy/expansion/multipole.py +++ b/sumpy/expansion/multipole.py @@ -203,7 +203,7 @@ class VolumeTaylorMultipoleExpansionBase(MultipoleExpansionBase): # └───⬠↑ # └─────┘ # - # For the second hyperplane, data is propogated rightwards first + # For the second hyperplane, data is propagated rightwards first # and then upwards second which is opposite to that of the first # hyperplane. # diff --git a/sumpy/qbx.py b/sumpy/qbx.py index b479e1db31b59d691a2c602080f5f7d3ef80055f..de2598881fc067dbe24e1b4e8f61b6510f2a3849 100644 --- a/sumpy/qbx.py +++ b/sumpy/qbx.py @@ -104,7 +104,7 @@ class LayerPotentialBase(KernelComputation, KernelCacheWrapper): # In LineTaylorLocalExpansion.evaluate, we can't run # postprocess_at_target because the coefficients are assigned # symbols and postprocess with a derivative will make them zero. - # Instead run postprocess here before the coeffients are assigned. + # Instead run postprocess here before the coefficients are assigned. coefficients = [tgt_knl.postprocess_at_target(coeff, bvec) for coeff in coefficients] diff --git a/test/test_codegen.py b/test/test_codegen.py index 815b92dd9dccf8111cd5f7f1204f99e9485911de..9c092e58efc85ed52e128226518609ab908ebf71 100644 --- a/test/test_codegen.py +++ b/test/test_codegen.py @@ -20,13 +20,15 @@ OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. """ - +import pytest import sys import logging logger = logging.getLogger(__name__) +# {{{ test_symbolic_assignment_name_uniqueness + def test_symbolic_assignment_name_uniqueness(): # https://gitlab.tiker.net/inducer/sumpy/issues/13 from sumpy.assignment_collection import SymbolicAssignmentCollection @@ -43,6 +45,10 @@ def test_symbolic_assignment_name_uniqueness(): assert len(sac.assignments) == 3 +# }}} + + +# {{{ test_line_taylor_coeff_growth def test_line_taylor_coeff_growth(): # Regression test for LineTaylorLocalExpansion. @@ -70,15 +76,16 @@ def test_line_taylor_coeff_growth(): max_order = 2 assert np.polyfit(np.log(indices), np.log(counts), deg=1)[0] < max_order +# }}} + # You can test individual routines by typing -# $ python test_fmm.py "test_sumpy_fmm(cl.create_some_context)" +# $ python test_codegen.py 'test_line_taylor_coeff_growth()' if __name__ == "__main__": if len(sys.argv) > 1: exec(sys.argv[1]) else: - from pytest import main - main([__file__]) + pytest.main([__file__]) # vim: fdm=marker diff --git a/test/test_cse.py b/test/test_cse.py index ba9986bed133a9155e78ee644a0a459d16e5862f..68352c9059e8587cd8c2da63a32075ef100841cc 100644 --- a/test/test_cse.py +++ b/test/test_cse.py @@ -67,26 +67,22 @@ DAMAGE. import pytest import sys -from sumpy.symbolic import ( - Add, Pow, exp, sqrt, symbols, sympify, cos, sin, Function, USE_SYMENGINE) +import sumpy.symbolic as sym +from sumpy.cse import cse, preprocess_for_cse, postprocess_for_cse -if not USE_SYMENGINE: +if not sym.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) +import logging +logger = logging.getLogger(__name__) - -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") +w, x, y, z = sym.symbols("w,x,y,z") +x0, x1, x2, x3, x4, x5, x6, x7, x8, x9, x10, x11, x12 = sym.symbols("x:13") sympyonly = ( - pytest.mark.skipif(USE_SYMENGINE, reason="uses a sympy-only feature")) - - -# Dummy "optimization" functions for testing. + pytest.mark.skipif(sym.USE_SYMENGINE, reason="uses a sympy-only feature")) def opt1(expr): @@ -97,6 +93,8 @@ def opt2(expr): return expr*z +# {{{ test_preprocess_for_cse + def test_preprocess_for_cse(): assert preprocess_for_cse(x, [(opt1, None)]) == x + y assert preprocess_for_cse(x, [(None, opt1)]) == x @@ -105,6 +103,10 @@ def test_preprocess_for_cse(): assert preprocess_for_cse( x, [(opt1, None), (opt2, None)]) == (x + y)*z +# }}} + + +# {{{ test_postprocess_for_cse def test_postprocess_for_cse(): assert postprocess_for_cse(x, [(opt1, None)]) == x @@ -115,19 +117,27 @@ def test_postprocess_for_cse(): assert postprocess_for_cse( x, [(None, opt1), (None, opt2)]) == x*z + y +# }}} + + +# {{{ test_cse_single def test_cse_single(): # Simple substitution. - e = Add(Pow(x + y, 2), sqrt(x + y)) + e = sym.Add(sym.Pow(x + y, 2), sym.sqrt(x + y)) substs, reduced = cse([e]) assert substs == [(x0, x + y)] - assert reduced == [sqrt(x0) + x0**2] + assert reduced == [sym.sqrt(x0) + x0**2] + +# }}} + +# {{{ @sympyonly def test_cse_not_possible(): # No substitution possible. - e = Add(x, y) + e = sym.Add(x, y) substs, reduced = cse([e]) assert substs == [] assert reduced == [x + y] @@ -136,34 +146,46 @@ def test_cse_not_possible(): + meijerg((1, 3), (y, 4), (5,), [], x)) assert cse(eq) == ([], [eq]) +# }}} + + +# {{{ test_nested_substitution def test_nested_substitution(): # Substitution within a substitution. - e = Add(Pow(w*x + y, 2), sqrt(w*x + y)) + e = sym.Add(sym.Pow(w*x + y, 2), sym.sqrt(w*x + y)) substs, reduced = cse([e]) assert substs == [(x0, w*x + y)] - assert reduced == [sqrt(x0) + x0**2] + assert reduced == [sym.sqrt(x0) + x0**2] +# }}} + + +# {{{ test_subtraction_opt @sympyonly def test_subtraction_opt(): # Make sure subtraction is optimized. - e = (x - y)*(z - y) + exp((x - y)*(z - y)) + e = (x - y)*(z - y) + sym.exp((x - y)*(z - y)) substs, reduced = cse( [e], optimizations=[(cse_opts.sub_pre, cse_opts.sub_post)]) assert substs == [(x0, (x - y)*(y - z))] - assert reduced == [-x0 + exp(-x0)] - e = -(x - y)*(z - y) + exp(-(x - y)*(z - y)) + assert reduced == [-x0 + sym.exp(-x0)] + e = -(x - y)*(z - y) + sym.exp(-(x - y)*(z - y)) substs, reduced = cse( [e], optimizations=[(cse_opts.sub_pre, cse_opts.sub_post)]) assert substs == [(x0, (x - y)*(y - z))] - assert reduced == [x0 + exp(x0)] + assert reduced == [x0 + sym.exp(x0)] # issue 4077 n = -1 + 1/x e = n/x/(-n)**2 - 1/n/x assert cse(e, optimizations=[(cse_opts.sub_pre, cse_opts.sub_post)]) == \ ([], [0]) +# }}} + + +# {{{ test_multiple_expressions def test_multiple_expressions(): e1 = (x + y)*z @@ -181,7 +203,7 @@ def test_multiple_expressions(): rsubsts, _ = cse(reversed(l_)) assert substs == rsubsts assert reduced == [x1, x1 + z, x0] - f = Function("f") + f = sym.Function("f") l_ = [f(x - z, y - z), x - z, y - z] substs, reduced = cse(l_) rsubsts, _ = cse(reversed(l_)) @@ -195,31 +217,46 @@ def test_multiple_expressions(): assert cse([x*y, z + x*y, x*y*z + 3]) == \ ([(x0, x*y)], [x0, z + x0, 3 + x0*z]) +# }}} + + +# {{{ test_issue_4203 def test_issue_4203(): - assert cse(sin(x**x)/x**x) == ([(x0, x**x)], [sin(x0)/x0]) + assert cse(sym.sin(x**x)/x**x) == ([(x0, x**x)], [sym.sin(x0)/x0]) + +# }}} +# {{{ test_dont_cse_subs + def test_dont_cse_subs(): - f = Function("f") - g = Function("g") + f = sym.Function("f") + g = sym.Function("g") name_val, (expr,) = cse(f(x+y).diff(x) + g(x+y).diff(x)) assert name_val == [] +# }}} + + +# {{{ test_dont_cse_derivative def test_dont_cse_derivative(): - from sumpy.symbolic import Derivative - f = Function("f") + f = sym.Function("f") - deriv = Derivative(f(x+y), x) + deriv = sym.Derivative(f(x+y), x) name_val, (expr,) = cse(x + y + deriv) assert name_val == [] assert expr == x + y + deriv +# }}} + + +# {{{ test_pow_invpow def test_pow_invpow(): assert cse(1/x**2 + x**2) == \ @@ -228,39 +265,64 @@ def test_pow_invpow(): ([(x0, x**2), (x1, 1/x0)], [x0 + x1*(x1 + 1)]) assert cse(1/x**2 + (1 + 1/x**2)*x**2) == \ ([(x0, x**2), (x1, 1/x0)], [x0*(x1 + 1) + x1]) - assert cse(cos(1/x**2) + sin(1/x**2)) == \ - ([(x0, x**(-2))], [sin(x0) + cos(x0)]) - assert cse(cos(x**2) + sin(x**2)) == \ - ([(x0, x**2)], [sin(x0) + cos(x0)]) + assert cse(sym.cos(1/x**2) + sym.sin(1/x**2)) == \ + ([(x0, x**(-2))], [sym.sin(x0) + sym.cos(x0)]) + assert cse(sym.cos(x**2) + sym.sin(x**2)) == \ + ([(x0, x**2)], [sym.sin(x0) + sym.cos(x0)]) assert cse(y/(2 + x**2) + z/x**2/y) == \ ([(x0, x**2)], [y/(x0 + 2) + z/(x0*y)]) - assert cse(exp(x**2) + x**2*cos(1/x**2)) == \ - ([(x0, x**2)], [x0*cos(1/x0) + exp(x0)]) + assert cse(sym.exp(x**2) + x**2*sym.cos(1/x**2)) == \ + ([(x0, x**2)], [x0*sym.cos(1/x0) + sym.exp(x0)]) assert cse((1 + 1/x**2)/x**2) == \ ([(x0, x**(-2))], [x0*(x0 + 1)]) assert cse(x**(2*y) + x**(-2*y)) == \ ([(x0, x**(2*y))], [x0 + 1/x0]) +# }}} + + +# {{{ test_issue_4499 @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( - *(a, a + S(1)/2, 2*a, b, 2*a - b + 1, (sqrt(z)/2)**(-2*a + 1)*B(2*a - - b, sqrt(z))*B(b - 1, sqrt(z))*G(b)*G(2*a - b + 1), # noqa - sqrt(z)*(sqrt(z)/2)**(-2*a + 1)*B(b, sqrt(z))*B(2*a - b, # noqa - sqrt(z))*G(b)*G(2*a - b + 1), sqrt(z)*(sqrt(z)/2)**(-2*a + 1)*B(b - 1, - sqrt(z))*B(2*a - b + 1, sqrt(z))*G(b)*G(2*a - b + 1), - (sqrt(z)/2)**(-2*a + 1)*B(b, sqrt(z))*B(2*a - b + 1, # noqa - sqrt(z))*G(b)*G(2*a - b + 1), 1, 0, S(1)/2, z/2, -b + 1, -2*a + b, # noqa - -2*a)) # noqa + + B = sym.Function("B") # noqa: N806 + G = sym.Function("G") # noqa: N806 + t = Tuple(*( + a, + a + S(1)/2, + 2*a, + b, + 2*a - b + 1, + (sym.sqrt(z)/2)**(-2*a + 1) + * B(2*a-b, sym.sqrt(z)) + * B(b - 1, sym.sqrt(z))*G(b)*G(2*a - b + 1), + sym.sqrt(z)*(sym.sqrt(z)/2)**(-2*a + 1) + * B(b, sym.sqrt(z)) + * B(2*a - b, sym.sqrt(z))*G(b)*G(2*a - b + 1), + sym.sqrt(z)*(sym.sqrt(z)/2)**(-2*a + 1) + * B(b - 1, sym.sqrt(z)) + * B(2*a - b + 1, sym.sqrt(z))*G(b)*G(2*a - b + 1), + (sym.sqrt(z)/2)**(-2*a + 1) + * B(b, sym.sqrt(z)) + * B(2*a - b + 1, sym.sqrt(z))*G(b)*G(2*a - b + 1), + 1, + 0, + S(1)/2, + z/2, + -b + 1, + -2*a + b, + -2*a)) c = cse(t) assert len(c[0]) == 11 +# }}} + + +# {{{ test_issue_6169 @sympyonly def test_issue_6169(): @@ -271,14 +333,17 @@ def test_issue_6169(): # mechanism assert sub_post(sub_pre((-x - y)*z - x - y)) == -z*(x + y) - x - y +# }}} + + +# {{{ test_cse_Indexed @sympyonly -def test_cse_Indexed(): # noqa +def test_cse_indexed(): from sympy import IndexedBase, Idx len_y = 5 y = IndexedBase("y", shape=(len_y,)) x = IndexedBase("x", shape=(len_y,)) - Dy = IndexedBase("Dy", shape=(len_y-1,)) # noqa i = Idx("i", len_y-1) expr1 = (y[i+1]-y[i])/(x[i+1]-x[i]) @@ -286,9 +351,13 @@ def test_cse_Indexed(): # noqa replacements, reduced_exprs = cse([expr1, expr2]) assert len(replacements) > 0 +# }}} + + +# {{{ test_Piecewise @sympyonly -def test_Piecewise(): # noqa +def test_piecewise(): from sympy import Piecewise, Eq f = Piecewise((-z + x*y, Eq(y, 0)), (-z - x*y, True)) ans = cse(f) @@ -296,41 +365,57 @@ def test_Piecewise(): # noqa [Piecewise((x0+x1, Eq(y, 0)), (x0 - x1, True))]) assert ans == actual_ans +# }}} + + +# {{{ test_name_conflict def test_name_conflict(): z1 = x0 + y z2 = x2 + x3 - l_ = [cos(z1) + z1, cos(z2) + z2, x0 + x2] + l_ = [sym.cos(z1) + z1, sym.cos(z2) + z2, x0 + x2] substs, reduced = cse(l_) assert [e.subs(dict(substs)) for e in reduced] == l_ +# }}} + + +# {{{ test_name_conflict_cust_symbols def test_name_conflict_cust_symbols(): z1 = x0 + y z2 = x2 + x3 - l_ = [cos(z1) + z1, cos(z2) + z2, x0 + x2] - substs, reduced = cse(l_, symbols("x:10")) + l_ = [sym.cos(z1) + z1, sym.cos(z2) + z2, x0 + x2] + substs, reduced = cse(l_, sym.symbols("x:10")) assert [e.subs(dict(substs)) for e in reduced] == l_ +# }}} + + +# {{{ test_symbols_exhausted_error def test_symbols_exhausted_error(): - l_ = cos(x+y)+x+y+cos(w+y)+sin(w+y) - sym = [x, y, z] + l_ = sym.cos(x+y)+x+y+sym.cos(w+y)+sym.sin(w+y) + s = [x, y, z] with pytest.raises(ValueError): - print(cse(l_, symbols=sym)) + logger.info("cse:\n%s", cse(l_, symbols=s)) + +# }}} + +# {{{ test_issue_7840 @sympyonly def test_issue_7840(): # daveknippers' example - C393 = sympify( # noqa + C393 = sym.sympify( # noqa: N806 "Piecewise((C391 - 1.65, C390 < 0.5), (Piecewise((C391 - 1.65, \ C391 > 2.35), (C392, True)), True))" ) - C391 = sympify( # noqa + C391 = sym.sympify( # noqa: N806 "Piecewise((2.05*C390**(-1.03), C390 < 0.5), (2.5*C390**(-0.625), True))" ) - C393 = C393.subs("C391",C391) # noqa + C393 = C393.subs("C391", C391) # noqa: N806 # simple substitution sub = {} sub["C390"] = 0.703451854 @@ -345,7 +430,7 @@ def test_issue_7840(): assert ss_answer == cse_answer # GitRay's example - expr = sympify( + expr = sym.sympify( "Piecewise((Symbol('ON'), Equality(Symbol('mode'), Symbol('ON'))), \ (Piecewise((Piecewise((Symbol('OFF'), StrictLessThan(Symbol('x'), \ Symbol('threshold'))), (Symbol('ON'), S.true)), Equality(Symbol('mode'), \ @@ -357,6 +442,10 @@ def test_issue_7840(): # there should not be any replacements assert len(substitutions) < 1 +# }}} + + +# {{{ test_recursive_matching def test_recursive_matching(): assert cse([x+y, 2+x+y, x+y+z, 3+x+y+z]) == \ @@ -372,12 +461,16 @@ def test_recursive_matching(): assert cse([2*x*x, x*x*y, x*x*y*w, x*x*y*w*x0, x*x*y*w*x2]) == \ ([(x1, x**2), (x3, x1*y), (x4, w*x3)], [2*x1, x3, x4, x0*x4, x2*x4]) +# }}} + + +# You can test individual routines by typing +# $ python test_cse.py 'test_recursive_matching()' if __name__ == "__main__": if len(sys.argv) > 1: exec(sys.argv[1]) else: - from pytest import main - main([__file__]) + pytest.main([__file__]) # vim: fdm=marker