From be93795833fb37696e968378dc7653c7a8dce2c8 Mon Sep 17 00:00:00 2001
From: Stephan Seitz <stephan.seitz@fau.de>
Date: Wed, 22 Jan 2020 21:58:29 +0100
Subject: [PATCH] Cherry-pick changes from interpolation-refactoring

---
 pystencils_tests/test_interpolation.py | 150 ++++++++++++-------------
 1 file changed, 72 insertions(+), 78 deletions(-)

diff --git a/pystencils_tests/test_interpolation.py b/pystencils_tests/test_interpolation.py
index 0776a7de0..92f20a236 100644
--- a/pystencils_tests/test_interpolation.py
+++ b/pystencils_tests/test_interpolation.py
@@ -11,6 +11,8 @@ import itertools
 from os.path import dirname, join
 
 import numpy as np
+import pycuda.autoinit  # NOQA
+import pycuda.gpuarray as gpuarray
 import pytest
 import sympy
 
@@ -108,93 +110,79 @@ def test_rotate_interpolation(address_mode):
     pyconrad.imshow(out, "out " + address_mode)
 
 
-@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
+@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):
 
     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)
-
-            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:  # 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()
+    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]}")
+    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.importorskip('pycuda')
-    import pycuda.autoinit  # NOQA
-    import pycuda.gpuarray as gpuarray
+@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
     # shift = - sympy.Matrix([1.5, 1.5])
     shift = sympy.Matrix((0.0, 0.0))
 
-    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)
-            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 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]}")
+    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'])
@@ -222,7 +210,7 @@ def test_rotate_interpolation_size_change(address_mode):
 
 
 @pytest.mark.parametrize('address_mode, target',
-                         itertools.product(['border', 'wrap', 'clamp', 'mirror'], ['cpu', 'gpu']))
+                         itertools.product(['border', 'wrap', 'clamp', 'mirror'], ['cpu']))
 def test_field_interpolated(address_mode, target):
     x_f, y_f = pystencils.fields('x,y: float64 [2d]')
 
@@ -230,7 +218,7 @@ 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)
+    ast = pystencils.create_kernel(assignments, target=target)
     print(ast)
     print(pystencils.show_code(ast))
     kernel = ast.compile()
@@ -238,5 +226,11 @@ def test_field_interpolated(address_mode, target):
     out = np.zeros_like(lenna)
     kernel(x=lenna, y=out)
     pyconrad.imshow(out, "out " + address_mode)
-    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))
-- 
GitLab