From 995cc9ef8b24668772a6b3d94cff26c11d87ba8e Mon Sep 17 00:00:00 2001 From: Stephan Seitz <stephan.seitz@fau.de> Date: Thu, 16 Jan 2020 13:46:37 +0100 Subject: [PATCH] Big fat interpolation mode refactoring --- pystencils/interpolation_astnodes.py | 27 ++++++----- pystencils/transformations.py | 37 +++++++++------ pystencils_tests/test_interpolation.py | 64 +++++++++++--------------- 3 files changed, 66 insertions(+), 62 deletions(-) diff --git a/pystencils/interpolation_astnodes.py b/pystencils/interpolation_astnodes.py index 3ecc2a70a..46946194a 100644 --- a/pystencils/interpolation_astnodes.py +++ b/pystencils/interpolation_astnodes.py @@ -77,13 +77,14 @@ class Interpolator(object): allow_textures=True): super().__init__() - self.field = parent_field.new_field_with_different_name(parent_field.name) + self.field = parent_field self.field.field_type = pystencils.field.FieldType.CUSTOM self.address_mode = address_mode self.use_normalized_coordinates = use_normalized_coordinates - hash_str = hashlib.md5(f'{self.field}_{address_mode}_{interpolation_mode}'.encode()).hexdigest() + hash_str = hashlib.md5( + f'{self.field}_{address_mode}_{self.field.dtype}_{interpolation_mode}'.encode()).hexdigest() self.symbol = TypedSymbol('dummy_symbol_carrying_field' + self.field.name + hash_str, - 'dummy_symbol_carrying_field' + self.field.name + hash_str) + 'dummy_symbol_carrying_field' + self.field.name + self.field.name + hash_str) self.symbol.field = self.field self.symbol.interpolator = self self.allow_textures = allow_textures @@ -97,7 +98,8 @@ class Interpolator(object): def _hashable_contents(self): return (str(self.address_mode), str(type(self)), - self.symbol, + str(self.symbol.name), + hash(self.symbol.field), self.address_mode, self.use_normalized_coordinates) @@ -114,7 +116,7 @@ class Interpolator(object): return self.__str__() def __hash__(self): - return hash(self._hashable_contents) + return hash(self._hashable_contents()) @property def reproducible_hash(self): @@ -161,8 +163,8 @@ class InterpolatorAccess(TypedSymbol): obj.interpolator = symbol.interpolator return obj - def __hash__(self): - return hash((self.symbol, self.field, tuple(self.offsets), self.interpolator)) + def _hashable_contents(self): + return super()._hashable_content() + ((self.symbol, self.field, tuple(self.offsets), self.symbol.interpolator)) def __str__(self): return '%s_interpolator(%s)' % (self.field.name, ', '.join(str(o) for o in self.offsets)) @@ -431,12 +433,13 @@ class TextureCachedField(Interpolator): @property def _hashable_contents(self): - return (type(self), + return (str(self.address_mode), + str(type(self)), + str(self.symbol.name), + hash(self.symbol.field), self.address_mode, - self.filter_mode, - self.read_as_integer, - self.interpolation_mode, - self.use_normalized_coordinates) + self.use_normalized_coordinates, + 'T') def __hash__(self): return hash(self._hashable_contents) diff --git a/pystencils/transformations.py b/pystencils/transformations.py index a81d0fb0e..3b9b59ad0 100644 --- a/pystencils/transformations.py +++ b/pystencils/transformations.py @@ -1327,8 +1327,9 @@ def implement_interpolations(ast_node: ast.Node, def can_use_hw_interpolation(i): return (use_hardware_interpolation_for_f32 + and implement_by_texture_accesses and i.dtype == FLOAT32_T - and isinstance(i.interpolator, TextureCachedField)) + and isinstance(i.symbol.interpolator, TextureCachedField)) if implement_by_texture_accesses: @@ -1337,27 +1338,35 @@ def implement_interpolations(ast_node: ast.Node, try: import pycuda.driver as cuda texture = TextureCachedField.from_interpolator(i.interpolator) - i.interpolator = texture i.symbol.interpolator = texture if can_use_hw_interpolation(i): - texture.filter_mode = cuda.filter_mode.LINEAR + i.symbol.interpolator.filter_mode = cuda.filter_mode.LINEAR else: - texture.filter_mode = cuda.filter_mode.POINT - texture.read_as_integer = True + i.symbol.interpolator.filter_mode = cuda.filter_mode.POINT + i.symbol.interpolator.read_as_integer = True except Exception: pass ast_node.subs({old_i: i}) - from pystencils.math_optimizations import ReplaceOptim, optimize_ast + # from pystencils.math_optimizations import ReplaceOptim, optimize_ast - ImplementInterpolationByStencils = ReplaceOptim(lambda e: isinstance(e, InterpolatorAccess) - and not can_use_hw_interpolation(i), - lambda e: e.implementation_with_stencils() - ) + # ImplementInterpolationByStencils = ReplaceOptim(lambda e: isinstance(e, InterpolatorAccess) + # and not can_use_hw_interpolation(i), + # lambda e: e.implementation_with_stencils() + # ) - RemoveConjugate = ReplaceOptim(lambda e: isinstance(e, sp.conjugate), - lambda e: e.args[0] - ) - optimize_ast(ast_node, [RemoveConjugate, ImplementInterpolationByStencils]) + # RemoveConjugate = ReplaceOptim(lambda e: isinstance(e, sp.conjugate), + # lambda e: e.args[0] + # ) + if vectorize: + # TODO can be done in _interpolator_access_to_stencils field.absolute_access == simd_gather + raise NotImplementedError() + else: + substitutions = {i: i.implementation_with_stencils() + for i in interpolation_accesses if not can_use_hw_interpolation(i)} + if isinstance(ast_node, AssignmentCollection): + ast_node = ast_node.subs(substitutions) + else: + ast_node.subs(substitutions) return ast_node diff --git a/pystencils_tests/test_interpolation.py b/pystencils_tests/test_interpolation.py index fbb3067e6..948c63748 100644 --- a/pystencils_tests/test_interpolation.py +++ b/pystencils_tests/test_interpolation.py @@ -112,51 +112,43 @@ def test_rotate_interpolation(address_mode): @pytest.mark.parametrize('dtype', (np.int32, np.float32, np.float64)) @pytest.mark.parametrize('address_mode', ('border', 'wrap', 'clamp', 'mirror')) -def test_rotate_interpolation_gpu(dtype, address_mode): +@pytest.mark.parametrize('use_textures', ('use_textures', False)) +def test_rotate_interpolation_gpu(dtype, address_mode, use_textures): rotation_angle = sympy.pi / 5 scale = 1 - 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) + 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) - pystencils.show_code(ast) - kernel = ast.compile() + 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]}") - 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 as e: # 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() + 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', 'wrap', 'clamp', 'mirror']) -def test_shift_interpolation_gpu(address_mode): +@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): rotation_angle = 0 # sympy.pi / 5 scale = 1 -- GitLab