Commit f81d83b7 authored by Michael Kuron's avatar Michael Kuron
Browse files

Merge branch 'master' of i10git.cs.fau.de:pycodegen/pystencils into oldsympy

parents 16d731b5 46b7b88d
import os
import runpy
import sys
import warnings
import tempfile
import warnings
import nbformat
import pytest
......@@ -34,7 +34,7 @@ def add_path_to_ignore(path):
collect_ignore = [os.path.join(SCRIPT_FOLDER, "doc", "conf.py"),
os.path.join(SCRIPT_FOLDER, "pystencils", "opencl", "opencl.autoinit")]
os.path.join(SCRIPT_FOLDER, "pystencils", "opencl", "opencl.autoinit")]
add_path_to_ignore('pystencils_tests/benchmark')
add_path_to_ignore('_local_tmp')
......@@ -98,8 +98,6 @@ for root, sub_dirs, files in os.walk('.'):
collect_ignore.append(f)
class IPythonMockup:
def run_line_magic(self, *args, **kwargs):
pass
......
......@@ -5,7 +5,7 @@ from . import stencil as stencil
from .assignment import Assignment, assignment_from_stencil
from .data_types import TypedSymbol
from .datahandling import create_data_handling
from .display_utils import show_code, to_dot
from .display_utils import show_code, get_code_obj, get_code_str, to_dot
from .field import Field, FieldType, fields
from .kernel_decorator import kernel
from .kernelcreation import create_indexed_kernel, create_kernel, create_staggered_kernel
......@@ -25,7 +25,7 @@ __all__ = ['Field', 'FieldType', 'fields',
'TypedSymbol',
'make_slice',
'create_kernel', 'create_indexed_kernel', 'create_staggered_kernel',
'show_code', 'to_dot',
'show_code', 'to_dot', 'get_code_obj', 'get_code_str',
'AssignmentCollection',
'Assignment',
'assignment_from_stencil',
......
# -*- coding: utf-8 -*-
import numpy as np
import sympy as sp
from sympy.printing.latex import LatexPrinter
......@@ -24,9 +23,20 @@ def assignment_str(assignment):
if Assignment:
_old_new = sp.codegen.ast.Assignment.__new__
def _Assignment__new__(cls, lhs, rhs, *args, **kwargs):
if isinstance(lhs, (list, tuple, sp.Matrix)) and isinstance(rhs, (list, tuple, sp.Matrix)):
assert len(lhs) == len(rhs), f'{lhs} and {rhs} must have same length when performing vector assignment!'
return tuple(_old_new(cls, a, b, *args, **kwargs) for a, b in zip(lhs, rhs))
return _old_new(cls, lhs, rhs, *args, **kwargs)
Assignment.__str__ = assignment_str
Assignment.__new__ = _Assignment__new__
LatexPrinter._print_Assignment = print_assignment_latex
sp.MutableDenseMatrix.__hash__ = lambda self: hash(tuple(self))
else:
# back port for older sympy versions that don't have Assignment yet
......
from typing import Any, Dict, Optional
from typing import Any, Dict, Optional, Union
import sympy as sp
......@@ -35,10 +35,10 @@ def highlight_cpp(code: str):
return HTML(highlight(code, CppLexer(), HtmlFormatter()))
def show_code(ast: KernelFunction, custom_backend=None):
def get_code_obj(ast: Union[KernelFunction, KernelWrapper], custom_backend=None):
"""Returns an object to display generated code (C/C++ or CUDA)
Can either be displayed as HTML in Jupyter notebooks or printed as normal string.
Can either be displayed as HTML in Jupyter notebooks or printed as normal string.
"""
from pystencils.backends.cbackend import generate_c
......@@ -65,3 +65,17 @@ def show_code(ast: KernelFunction, custom_backend=None):
def __repr__(self):
return generate_c(self.ast, dialect=dialect, custom_backend=custom_backend)
return CodeDisplay(ast)
def get_code_str(ast, custom_backend=None):
return str(get_code_obj(ast, custom_backend))
def show_code(ast: Union[KernelFunction, KernelWrapper], custom_backend=None):
code = get_code_obj(ast, custom_backend)
try:
from IPython.display import display
display(code)
except Exception:
print(code)
"""
Light-weight wrapper around a compiled kernel
"""
import pystencils
class KernelWrapper:
def __init__(self, kernel, parameters, ast_node):
"""
Light-weight wrapper around a compiled kernel.
Can be called while still providing access to underlying AST.
"""
def __init__(self, kernel, parameters, ast_node: pystencils.astnodes.KernelFunction):
self.kernel = kernel
self.parameters = parameters
self.ast = ast_node
......@@ -16,4 +19,4 @@ class KernelWrapper:
@property
def code(self):
return str(pystencils.show_code(self.ast))
return pystencils.get_code_str(self.ast)
import itertools
from copy import copy
from typing import Any, Dict, Iterable, Iterator, List, Optional, Sequence, Set, Union
......@@ -43,6 +44,11 @@ class AssignmentCollection:
subexpressions = [Assignment(k, v)
for k, v in subexpressions.items()]
main_assignments = list(itertools.chain.from_iterable(
[(a if isinstance(a, Iterable) else [a]) for a in main_assignments]))
subexpressions = list(itertools.chain.from_iterable(
[(a if isinstance(a, Iterable) else [a]) for a in subexpressions]))
self.main_assignments = main_assignments
self.subexpressions = subexpressions
......
......@@ -18,8 +18,7 @@ def test_type_interference():
ast = pystencils.create_kernel(assignments)
code = str(pystencils.show_code(ast))
print(code)
code = str(pystencils.get_code_str(ast))
assert 'double a' in code
assert 'uint16_t b' in code
assert 'uint16_t f' in code
......
......@@ -17,16 +17,14 @@ def test_address_of():
}, {})
ast = pystencils.create_kernel(assignments)
code = pystencils.show_code(ast)
print(code)
pystencils.show_code(ast)
assignments = pystencils.AssignmentCollection({
y[0, 0]: cast_func(address_of(x[0, 0]), create_type('int64'))
}, {})
ast = pystencils.create_kernel(assignments)
code = pystencils.show_code(ast)
print(code)
pystencils.show_code(ast)
def test_address_of_with_cse():
......@@ -39,9 +37,8 @@ def test_address_of_with_cse():
}, {})
ast = pystencils.create_kernel(assignments)
code = pystencils.show_code(ast)
pystencils.show_code(ast)
assignments_cse = sympy_cse(assignments)
ast = pystencils.create_kernel(assignments_cse)
code = pystencils.show_code(ast)
print(code)
pystencils.show_code(ast)
import pytest
import sympy as sp
from pystencils import Assignment, AssignmentCollection
......@@ -40,3 +41,39 @@ def test_free_and_defined_symbols():
print(ac)
print(ac.__repr__)
def test_vector_assignments():
"""From #17 (https://i10git.cs.fau.de/pycodegen/pystencils/issues/17)"""
import pystencils as ps
import sympy as sp
a, b, c = sp.symbols("a b c")
assignments = ps.Assignment(sp.Matrix([a,b,c]), sp.Matrix([1,2,3]))
print(assignments)
def test_wrong_vector_assignments():
"""From #17 (https://i10git.cs.fau.de/pycodegen/pystencils/issues/17)"""
import pystencils as ps
import sympy as sp
a, b = sp.symbols("a b")
with pytest.raises(AssertionError,
match=r'Matrix(.*) and Matrix(.*) must have same length when performing vector assignment!'):
ps.Assignment(sp.Matrix([a,b]), sp.Matrix([1,2,3]))
def test_vector_assignment_collection():
"""From #17 (https://i10git.cs.fau.de/pycodegen/pystencils/issues/17)"""
import pystencils as ps
import sympy as sp
a, b, c = sp.symbols("a b c")
y, x = sp.Matrix([a,b,c]), sp.Matrix([1,2,3])
assignments = ps.AssignmentCollection({y: x})
print(assignments)
assignments = ps.AssignmentCollection([ps.Assignment(y,x)])
print(assignments)
......@@ -52,7 +52,7 @@ def test_complex_numbers(assignment, scalar_dtypes, target):
ast = pystencils.create_kernel(assignment,
target=target,
data_type=scalar_dtypes)
code = str(pystencils.show_code(ast))
code = pystencils.get_code_str(ast)
print(code)
assert "Not supported" not in code
......@@ -98,7 +98,7 @@ def test_complex_numbers_64(assignment, target):
ast = pystencils.create_kernel(assignment,
target=target,
data_type='double')
code = str(pystencils.show_code(ast))
code = pystencils.get_code_str(ast)
print(code)
assert "Not supported" not in code
......
......@@ -63,7 +63,7 @@ def test_boundary_check(with_cse):
print(assignments)
kernel_checked = ps.create_kernel(assignments, ghost_layers=0).compile()
print(ps.show_code(kernel_checked))
ps.show_code(kernel_checked)
# No SEGFAULT, please!!
kernel_checked(f=f_arr, g=g_arr)
......@@ -20,8 +20,8 @@ def test_cuda_known_functions():
})
ast = pystencils.create_kernel(assignments, 'gpu')
print(pystencils.show_code(ast))
pytest.importorskip('pycuda')
pystencils.show_code(ast)
kernel = ast.compile()
assert(kernel is not None)
......@@ -35,7 +35,7 @@ def test_cuda_but_not_c():
})
ast = pystencils.create_kernel(assignments, 'cpu')
print(pystencils.show_code(ast))
pystencils.show_code(ast)
def test_cuda_unknown():
......@@ -46,5 +46,4 @@ def test_cuda_unknown():
})
ast = pystencils.create_kernel(assignments, 'gpu')
code = str(pystencils.show_code(ast))
print(code)
pystencils.show_code(ast)
......@@ -3,6 +3,7 @@ from subprocess import CalledProcessError
import pytest
import sympy
import pycuda.driver
import pystencils
import pystencils.cpu.cpujit
import pystencils.gpucuda.cudajit
......@@ -31,7 +32,7 @@ def test_custom_backends_cpu():
z[0, 0], x[0, 0] * sympy.log(x[0, 0] * y[0, 0]))], [])
ast = pystencils.create_kernel(normal_assignments, target='cpu')
print(pystencils.show_code(ast, ScreamingBackend()))
pystencils.show_code(ast, ScreamingBackend())
with pytest.raises(CalledProcessError):
pystencils.cpu.cpujit.make_python_function(ast, custom_backend=ScreamingBackend())
......@@ -46,6 +47,6 @@ def test_custom_backends_gpu():
z[0, 0], x[0, 0] * sympy.log(x[0, 0] * y[0, 0]))], [])
ast = pystencils.create_kernel(normal_assignments, target='gpu')
print(pystencils.show_code(ast, ScreamingGpuBackend()))
pystencils.show_code(ast, ScreamingGpuBackend())
with pytest.raises(pycuda.driver.CompileError):
pystencils.gpucuda.cudajit.make_python_function(ast, custom_backend=ScreamingGpuBackend())
......@@ -12,7 +12,7 @@ def test_fast_sqrt():
assert len(insert_fast_sqrts(expr).atoms(fast_sqrt)) == 1
assert len(insert_fast_sqrts([expr])[0].atoms(fast_sqrt)) == 1
ast = ps.create_kernel(ps.Assignment(g[0, 0], insert_fast_sqrts(expr)), target='gpu')
code_str = str(ps.show_code(ast))
code_str = ps.get_code_str(ast)
assert '__fsqrt_rn' in code_str
expr = ps.Assignment(sp.Symbol("tmp"), 3 / sp.sqrt(f[0, 0] + f[1, 0]))
......@@ -21,7 +21,7 @@ def test_fast_sqrt():
ac = ps.AssignmentCollection([expr], [])
assert len(insert_fast_sqrts(ac).main_assignments[0].atoms(fast_inv_sqrt)) == 1
ast = ps.create_kernel(insert_fast_sqrts(ac), target='gpu')
code_str = str(ps.show_code(ast))
code_str = ps.get_code_str(ast)
assert '__frsqrt_rn' in code_str
......@@ -34,5 +34,5 @@ def test_fast_divisions():
assert len(insert_fast_divisions(expr).atoms(fast_division)) == 1
ast = ps.create_kernel(ps.Assignment(g[0, 0], insert_fast_divisions(expr)), target='gpu')
code_str = str(ps.show_code(ast))
code_str = ps.get_code_str(ast)
assert '__fdividef' in code_str
......@@ -49,7 +49,7 @@ def test_interpolation():
print(assignments)
ast = pystencils.create_kernel(assignments)
print(ast)
print(pystencils.show_code(ast))
pystencils.show_code(ast)
kernel = ast.compile()
pyconrad.imshow(lenna)
......@@ -69,7 +69,7 @@ def test_scale_interpolation():
print(assignments)
ast = pystencils.create_kernel(assignments)
print(ast)
print(pystencils.show_code(ast))
pystencils.show_code(ast)
kernel = ast.compile()
out = np.zeros_like(lenna)
......@@ -81,9 +81,9 @@ def test_scale_interpolation():
['border',
'clamp',
pytest.param('warp', marks=pytest.mark.xfail(
reason="requires interpolation-refactoring branch")),
reason="Fails on newer SymPy version due to complex conjugate()")),
pytest.param('mirror', marks=pytest.mark.xfail(
reason="requires interpolation-refactoring branch")),
reason="Fails on newer SymPy version due to complex conjugate()")),
])
def test_rotate_interpolation(address_mode):
"""
......@@ -100,7 +100,7 @@ def test_rotate_interpolation(address_mode):
print(assignments)
ast = pystencils.create_kernel(assignments)
print(ast)
print(pystencils.show_code(ast))
pystencils.show_code(ast)
kernel = ast.compile()
out = np.zeros_like(lenna)
......@@ -108,90 +108,93 @@ def test_rotate_interpolation(address_mode):
pyconrad.imshow(out, "out " + address_mode)
@pytest.mark.parametrize('dtype', (np.int32, np.float32, np.float64))
@pytest.mark.parametrize('address_mode', ('border', 'wrap', 'clamp', 'mirror'))
@pytest.mark.parametrize('use_textures', ('use_textures', False))
def test_rotate_interpolation_gpu(dtype, address_mode, use_textures):
@pytest.mark.parametrize('address_mode', ['border', 'wrap', 'clamp', 'mirror'])
def test_rotate_interpolation_gpu(address_mode):
pytest.importorskip('pycuda')
import pycuda.autoinit # NOQA
import pycuda.gpuarray as gpuarray
import pycuda.autoinit # noqa
rotation_angle = sympy.pi / 5
scale = 1
if dtype == np.int32:
lenna_gpu = gpuarray.to_gpu(
np.ascontiguousarray(lenna * 255, dtype))
else:
lenna_gpu = gpuarray.to_gpu(
np.ascontiguousarray(lenna, dtype))
x_f, y_f = pystencils.fields('x,y: %s [2d]' % type_map[dtype], ghost_layers=0)
transformed = scale * \
sympy.rot_axis3(rotation_angle)[:2, :2] * sympy.Matrix((x_, y_)) - sympy.Matrix([2, 2])
assignments = pystencils.AssignmentCollection({
y_f.center(): LinearInterpolator(x_f, address_mode=address_mode).at(transformed)
})
print(assignments)
ast = pystencils.create_kernel(assignments, target='gpu', use_textures_for_interpolation=use_textures)
print(ast)
print(pystencils.show_code(ast))
kernel = ast.compile()
out = gpuarray.zeros_like(lenna_gpu)
kernel(x=lenna_gpu, y=out)
pyconrad.imshow(out,
f"out {address_mode} texture:{use_textures} {type_map[dtype]}")
@pytest.mark.parametrize('address_mode', ['border', 'wrap',
pytest.param('warp', marks=pytest.mark.xfail(
reason="% printed as fmod on old sympy")),
pytest.param('mirror', marks=pytest.mark.xfail(
reason="% printed as fmod on old sympy")),
])
@pytest.mark.parametrize('dtype', [np.float64, np.float32, np.int32])
@pytest.mark.parametrize('use_textures', ('use_textures', False,))
def test_shift_interpolation_gpu(address_mode, dtype, use_textures):
sver = sympy.__version__.split(".")
if (int(sver[0]) == 1 and int(sver[1]) < 2) and address_mode in ['mirror', 'warp']:
pytest.skip()
previous_result = None
for dtype in [np.int32, np.float32, np.float64]:
if dtype == np.int32:
lenna_gpu = gpuarray.to_gpu(
np.ascontiguousarray(lenna * 255, dtype))
else:
lenna_gpu = gpuarray.to_gpu(
np.ascontiguousarray(lenna, dtype))
for use_textures in [True, False]:
x_f, y_f = pystencils.fields('x,y: %s [2d]' % type_map[dtype], ghost_layers=0)
transformed = scale * sympy.rot_axis3(rotation_angle)[:2, :2] * \
sympy.Matrix((x_, y_)) - sympy.Matrix([2, 2])
assignments = pystencils.AssignmentCollection({
y_f.center(): LinearInterpolator(x_f, address_mode=address_mode).at(transformed)
})
print(assignments)
ast = pystencils.create_kernel(assignments, target='gpu', use_textures_for_interpolation=use_textures)
print(ast)
pystencils.show_code(ast)
kernel = ast.compile()
out = gpuarray.zeros_like(lenna_gpu)
kernel(x=lenna_gpu, y=out)
pyconrad.imshow(out,
f"out {address_mode} texture:{use_textures} {type_map[dtype]}")
skimage.io.imsave(f"/tmp/out {address_mode} texture:{use_textures} {type_map[dtype]}.tif",
np.ascontiguousarray(out.get(), np.float32))
if previous_result is not None:
try:
assert np.allclose(previous_result[4:-4, 4:-4], out.get()[4:-4, 4:-4], rtol=100, atol=1e-3)
except AssertionError: # NOQA
print("Max error: %f" % np.max(previous_result - out.get()))
# pyconrad.imshow(previous_result - out.get(), "Difference image")
# raise e
previous_result = out.get()
@pytest.mark.parametrize('address_mode', ['border', 'wrap', 'clamp', 'mirror'])
def test_shift_interpolation_gpu(address_mode):
pytest.importorskip('pycuda')
import pycuda.autoinit # NOQA
import pycuda.gpuarray as gpuarray
import pycuda.autoinit # noqa
rotation_angle = 0 # sympy.pi / 5
scale = 1
# shift = - sympy.Matrix([1.5, 1.5])
shift = sympy.Matrix((0.0, 0.0))
if dtype == np.int32:
lenna_gpu = gpuarray.to_gpu(
np.ascontiguousarray(lenna * 255, dtype))
else:
lenna_gpu = gpuarray.to_gpu(
np.ascontiguousarray(lenna, dtype))
x_f, y_f = pystencils.fields('x,y: %s [2d]' % type_map[dtype], ghost_layers=0)
if use_textures:
transformed = scale * sympy.rot_axis3(rotation_angle)[:2, :2] * sympy.Matrix((x_, y_)) + shift
else:
transformed = scale * sympy.rot_axis3(rotation_angle)[:2, :2] * sympy.Matrix((x_, y_)) + shift
assignments = pystencils.AssignmentCollection({
y_f.center(): LinearInterpolator(x_f, address_mode=address_mode).at(transformed)
})
# print(assignments)
ast = pystencils.create_kernel(assignments, target='gpu', use_textures_for_interpolation=use_textures)
# print(ast)
print(pystencils.show_code(ast))
kernel = ast.compile()
out = gpuarray.zeros_like(lenna_gpu)
kernel(x=lenna_gpu, y=out)
pyconrad.imshow(out,
f"out {address_mode} texture:{use_textures} {type_map[dtype]}")
for dtype in [np.float64, np.float32, np.int32]:
if dtype == np.int32:
lenna_gpu = gpuarray.to_gpu(
np.ascontiguousarray(lenna * 255, dtype))
else:
lenna_gpu = gpuarray.to_gpu(
np.ascontiguousarray(lenna, dtype))
for use_textures in [True, False]:
x_f, y_f = pystencils.fields('x,y: %s [2d]' % type_map[dtype], ghost_layers=0)
if use_textures:
transformed = scale * sympy.rot_axis3(rotation_angle)[:2, :2] * sympy.Matrix((x_, y_)) + shift
else:
transformed = scale * sympy.rot_axis3(rotation_angle)[:2, :2] * sympy.Matrix((x_, y_)) + shift
assignments = pystencils.AssignmentCollection({
y_f.center(): LinearInterpolator(x_f, address_mode=address_mode).at(transformed)
})
# print(assignments)
ast = pystencils.create_kernel(assignments, target='gpu', use_textures_for_interpolation=use_textures)
# print(ast)
pystencils.show_code(ast)
kernel = ast.compile()
out = gpuarray.zeros_like(lenna_gpu)
kernel(x=lenna_gpu, y=out)
pyconrad.imshow(out,
f"out {address_mode} texture:{use_textures} {type_map[dtype]}")
skimage.io.imsave(f"/tmp/out {address_mode} texture:{use_textures} {type_map[dtype]}.tif",
np.ascontiguousarray(out.get(), np.float32))
@pytest.mark.parametrize('address_mode', ['border', 'clamp'])
......@@ -210,7 +213,7 @@ def test_rotate_interpolation_size_change(address_mode):
print(assignments)
ast = pystencils.create_kernel(assignments)
print(ast)
print(pystencils.show_code(ast))
pystencils.show_code(ast)
kernel = ast.compile()
out = np.zeros((100, 150), np.float64)
......@@ -219,7 +222,7 @@ def test_rotate_interpolation_size_change(address_mode):
@pytest.mark.parametrize('address_mode, target',
itertools.product(['border', 'wrap', 'clamp', 'mirror'], ['cpu']))
itertools.product(['border', 'wrap', 'clamp', 'mirror'], ['cpu', 'gpu']))
def test_field_interpolated(address_mode, target):
x_f, y_f = pystencils.fields('x,y: float64 [2d]')
......@@ -227,19 +230,13 @@ def test_field_interpolated(address_mode, target):
y_f.center(): x_f.interpolated_access([0.5 * x_ + 2.7, 0.25 * y_ + 7.2], address_mode=address_mode)
})
print(assignments)
ast = pystencils.create_kernel(assignments, target=target)
ast = pystencils.create_kernel(assignments)
print(ast)
print(pystencils.show_code(ast))
pystencils.show_code(ast)
kernel = ast.compile()
out = np.zeros_like(lenna)
kernel(x=lenna, y=out)
pyconrad.imshow(out, "out " + address_mode)
def test_spatial_derivative():
x, y = pystencils.fields('x, y: float32[2d]')
tx, ty = pystencils.fields('t_x, t_y: float32[2d]')
diff = sympy.diff(x.interpolated_access((tx.center, ty.center)), tx.center)
print("diff: " + str(diff))
kernel(x=lenna, y=out)
pyconrad.imshow(out, "out " + address_mode)
import numpy as np
from pystencils import show_code
from pystencils import get_code_obj
from pystencils.astnodes import Block, KernelFunction, SympyAssignment
from pystencils.cpu import make_python_function