WIP: Astnodes for interpolation

This PR needs maybe still needs some clean-up. However, it would be good to recieve already some feed-back.

What works:

• Using CUDA textures
• Using HW accelerated interpolation for float32 textures
• Implement linear interpolations either via software (CPU, GPU), texture accesses without HW-interpolation but HW boundary handling
• Adding transformed coordinate systems to fields

What does not work:

• HW boundary handling for CUDA textures for the boundary handling modes mirror and wrap (apparently they have been removed from CUDA's API but are still present in pycuda. Now there's only

cudaBoundaryModeZero = 0
    Zero boundary mode
cudaBoundaryModeClamp = 1
    Clamp boundary mode
cudaBoundaryModeTrap = 2
    Trap boundary mode

Wtf is trap boundary mode? Nothing is documented so we can only experiment.

What kind of works:

• B-Spline interpolation on GPU using this repo as a submodule (http://www.dannyruijters.nl/cubicinterpolation/), to lazy for tests. Don't know how to prove correctness
• Textures for dtypes with itemsize > 4. PyCUDA has helper header (https://github.com/inducer/pycuda/blob/master/pycuda/cuda/pycuda-helpers.hpp) that loads doubles by two int fetches. However, this hack seems to be only working if we add a 0.5 offset and make all functions in this header accept float.

pystencils/backends/cbackend.py

@@ -9,11 +9,12 @@ from sympy.printing.ccode import C89CodePrinter
 from pystencils.astnodes import KernelFunction, Node
 from pystencils.cpu.vectorization import vec_all, vec_any
 from pystencils.data_types import (
-    PointerType, VectorType, address_of, cast_func, create_type, get_type_of_expression, reinterpret_cast_func,
-    vector_memory_access)
+    PointerType, VectorType, address_of, cast_func, create_type, get_type_of_expression,
+    reinterpret_cast_func, vector_memory_access)
 from pystencils.fast_approximation import fast_division, fast_inv_sqrt, fast_sqrt
 from pystencils.integer_functions import (
-    bit_shift_left, bit_shift_right, bitwise_and, bitwise_or, bitwise_xor, int_div, int_power_of_2, modulo_ceil)
+    bit_shift_left, bit_shift_right, bitwise_and, bitwise_or, bitwise_xor,
+    int_div, int_power_of_2, modulo_ceil)
 from pystencils.kernelparameters import FieldPointerSymbol
 
 try:
@@ -124,6 +125,12 @@ class CustomCodeNode(Node):
     def undefined_symbols(self):
         return self._symbols_read - self._symbols_defined
 
+    def __eq___(self, other):
+        return self._code == other._code
+
+    def __hash__(self):
+        return hash(self._code)
+
 
 class PrintNode(CustomCodeNode):
     # noinspection SpellCheckingInspection
@@ -168,7 +175,8 @@ class CBackend:
         raise NotImplementedError(self.__class__ + " does not support node of type " + str(type(node)))
 
     def _print_KernelFunction(self, node):
-        function_arguments = ["%s %s" % (str(s.symbol.dtype), s.symbol.name) for s in node.get_parameters()]
+        parameters = node.get_parameters()
+        function_arguments = ["%s %s" % (str(s.symbol.dtype), s.symbol.name) for s in parameters]
         launch_bounds = ""
         if self.__class__ == 'cuda':
             max_threads = node.indexing.max_threads_per_block()
@@ -247,6 +255,12 @@ class CBackend:
     def _print_CustomCodeNode(self, node):
         return node.get_code(self._dialect, self._vector_instruction_set)
 
+    def _print_SourceCodeComment(self, node):
+        return "/* " + node.text + " */"
+
+    def _print_EmptyLine(self, node):
+        return ""
+
     def _print_Conditional(self, node):
         cond_type = get_type_of_expression(node.condition_expr)
         if isinstance(cond_type, VectorType):
@@ -369,10 +383,61 @@ class CustomSympyPrinter(CCodePrinter):
                     res += ".0f"
                 else:
                     res += "f"
+            else:
+                if '.' not in res:
+                    res += "."
             return res
         else:
             return res
 
+    def _print_Sum(self, expr):
+        template = jinja2.Template(
+            """[&]() {
+    {{dtype}} sum = ({{dtype}}) 0;
+    for ( {{iterator_dtype}} {{var}} = {{start}}; {{condition}}; {{var}} += {{increment}} ) {
+        sum += {{expr}};
+    }
+    return sum;
+}()""")
+        var = expr.limits[0][0]
+        start = expr.limits[0][1]
+        end = expr.limits[0][2]
+        code = template.render(
+            dtype=get_type_of_expression(expr.args[0]),
+            iterator_dtype='int',
+            var=self._print(var),
+            start=self._print(start),
+            end=self._print(end),
+            expr=self._print(expr.function),
+            increment=str(1),
+            condition=self._print(var) + ' <= ' + self._print(end)  # if start < end else '>='
+        )
+        return code
+
+    def _print_Product(self, expr):
+        template = jinja2.Template(
+            """[&]() {
+    {{dtype}} product = ({{dtype}}) 1;
+    for ( {{iterator_dtype}} {{var}} = {{start}}; {{condition}}; {{var}} += {{increment}} ) {
+        product *= {{expr}};
+    }
+    return product;
+}()""")
+        var = expr.limits[0][0]
+        start = expr.limits[0][1]
+        end = expr.limits[0][2]
+        code = template.render(
+            dtype=get_type_of_expression(expr.args[0]),
+            iterator_dtype='int',
+            var=self._print(var),
+            start=self._print(start),
+            end=self._print(end),
+            expr=self._print(expr.function),
+            increment=str(1),
+            condition=self._print(var) + ' <= ' + self._print(end)  # if start < end else '>='
+        )
+        return code
+
     _print_Max = C89CodePrinter._print_Max
     _print_Min = C89CodePrinter._print_Min