.gitlab-ci.yml

 stages:
   - pretest
   - test
+  - nightly
+  - docs
   - deploy
 
 
+# --------------------------  Templates ------------------------------------------------------------------------------------
+
+# Base configuration for jobs meant to run at every commit
+.every-commit:
+  rules:
+    - if: $CI_PIPELINE_SOURCE != "schedule"
+
+# Configuration for jobs meant to run on each commit to pycodegen/pystencils/master
+.every-commit-master:
+  rules:
+    - if: '$CI_PIPELINE_SOURCE != "schedule" && $CI_PROJECT_PATH == "pycodegen/pystencils" && $CI_COMMIT_BRANCH == "master"'
+
+# Base configuration for jobs meant to run at a schedule
+.scheduled:
+  rules:
+    - if: $CI_PIPELINE_SOURCE == "schedule"
+
 # --------------------------  Tests ------------------------------------------------------------------------------------
 
 # Normal test - runs on every commit all but "long run" tests
 tests-and-coverage:
   stage: pretest
-  except:
-    variables:
-      - $ENABLE_NIGHTLY_BUILDS
+  extends: .every-commit
   image: i10git.cs.fau.de:5005/pycodegen/pycodegen/full
   before_script:
     - pip install -e .
@@ -65,9 +82,7 @@ tests-and-coverage-with-longrun:
 # pipeline with latest python version
 latest-python:
   stage: test
-  except:
-    variables:
-      - $ENABLE_NIGHTLY_BUILDS
+  extends: .every-commit
   image: i10git.cs.fau.de:5005/pycodegen/pycodegen/latest_python
   before_script:
     - pip install -e .
@@ -92,9 +107,6 @@ latest-python:
 # Minimal tests in windows environment
 #minimal-windows:
 #  stage: test
-#  except:
-#    variables:
-#      - $ENABLE_NIGHTLY_BUILDS
 #  tags:
 #    - win
 #  script:
@@ -108,9 +120,7 @@ latest-python:
 
 ubuntu:
   stage: test
-  except:
-    variables:
-      - $ENABLE_NIGHTLY_BUILDS
+  extends: .every-commit
   image: i10git.cs.fau.de:5005/pycodegen/pycodegen/ubuntu
   before_script:
     - ln -s /usr/include/locale.h /usr/include/xlocale.h
@@ -134,9 +144,7 @@ ubuntu:
 
 .multiarch_template:
   stage: test
-  except:
-    variables:
-      - $ENABLE_NIGHTLY_BUILDS
+  extends: .every-commit
   before_script: &multiarch_before_script
     # - pip3 install -v .
     - export PYTHONPATH=src
@@ -187,7 +195,7 @@ arm64v9:
   image: i10git.cs.fau.de:5005/pycodegen/pycodegen/arm64
   before_script:
     - *multiarch_before_script
-    - sed -i s/march=native/march=armv8-a+sve+sme/g ~/.config/pystencils/config.json
+    - sed -i s/march=native/march=armv9-a+sve2+sme/g ~/.config/pystencils/config.json
     - sed -i s/g\+\+/clang++/g ~/.config/pystencils/config.json
 
 riscv64:
@@ -205,9 +213,7 @@ riscv64:
 
 minimal-conda:
   stage: pretest
-  except:
-    variables:
-      - $ENABLE_NIGHTLY_BUILDS
+  extends: .every-commit
   image: i10git.cs.fau.de:5005/pycodegen/pycodegen/minimal_conda
   before_script:
     - pip install -e .
@@ -220,9 +226,7 @@ minimal-conda:
 
 minimal-sympy-master:
   stage: test
-  except:
-    variables:
-      - $ENABLE_NIGHTLY_BUILDS
+  extends: .every-commit
   image: i10git.cs.fau.de:5005/pycodegen/pycodegen/minimal_conda
   before_script:
     - pip install -e .
@@ -279,14 +283,42 @@ pycodegen-integration:
     reports:
       junit: pycodegen/*/report.xml
 
+
+# -------------------- Scheduled Tasks --------------------------------------------------------------------------
+
+
+# Nightly test against the latest (pre-release) version of SymPy published on PyPI
+nightly-sympy:
+  stage: nightly
+  needs: []
+  extends: .scheduled
+  image: i10git.cs.fau.de:5005/pycodegen/pycodegen/latest_python
+  before_script:
+    - pip install -e .
+    - pip install --upgrade --pre sympy
+  script:
+    - env
+    - pip list
+    - export NUM_CORES=$(nproc --all)
+    - mkdir -p ~/.config/matplotlib
+    - echo "backend:template" > ~/.config/matplotlib/matplotlibrc
+    - mkdir public
+    - pytest -v -n $NUM_CORES -m "not longrun" --junitxml=report.xml
+  tags:
+    - docker
+    - AVX
+    - cuda
+  artifacts:
+    when: always
+    reports:
+      junit: report.xml
+
 # -------------------- Linter & Documentation --------------------------------------------------------------------------
 
 
 flake8-lint:
   stage: pretest
-  except:
-    variables:
-      - $ENABLE_NIGHTLY_BUILDS
+  extends: .every-commit
   image: i10git.cs.fau.de:5005/pycodegen/pycodegen/full
   script:
     - flake8 src/pystencils
@@ -295,8 +327,10 @@ flake8-lint:
 
 
 build-documentation:
-  stage: test
+  stage: docs
+  extends: .every-commit
   image: i10git.cs.fau.de:5005/pycodegen/pycodegen/documentation
+  needs: []
   before_script:
     - pip install -e .
   script:
@@ -312,7 +346,9 @@ build-documentation:
 
 pages:
   image: i10git.cs.fau.de:5005/pycodegen/pycodegen/full
+  extends: .every-commit-master
   stage: deploy
+  needs: ["tests-and-coverage", "build-documentation"]
   script:
     - ls -l
     - mv coverage_report html_doc
@@ -322,5 +358,3 @@ pages:
       - public
   tags:
     - docker
-  only:
-    - master@pycodegen/pystencils

doc/notebooks/02_tutorial_basic_kernels.ipynb

 %% Cell type:code id: tags:
 
 ``` python
 from pystencils.session import *
 ```
 
 %% Cell type:markdown id: tags:
 
 # Tutorial 02: Basic Kernel generation with *pystencils*
 
 Now that you have an [overview of pystencils](01_tutorial_getting_started.ipynb),
 this tutorial shows in more detail how to formulate, optimize and run stencil kernels.
 
 ## 1) Kernel Definition
 
 ### a) Defining kernels with assignment lists and the  `kernel` decorator
 
 *pystencils* gets a symbolic formulation of the kernel. This can be either an `Assignment` or a sequence of `Assignment`s that follow a set of restrictions.
 
 Lets first create a kernel that consists of multiple assignments:
 
 %% Cell type:code id: tags:
 
 ``` python
 src_arr = np.zeros([20, 30])
 dst_arr = np.zeros_like(src_arr)
 
 dst, src = ps.fields(dst=dst_arr, src=src_arr)
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 grad_x, grad_y = sp.symbols("grad_x, grad_y")
 
 symbolic_description = [
     ps.Assignment(grad_x, (src[1, 0] - src[-1, 0]) / 2),
     ps.Assignment(grad_y, (src[0, 1] - src[0, -1]) / 2),
     ps.Assignment(dst[0, 0], grad_x + grad_y),
 ]
 kernel = ps.create_kernel(symbolic_description)
 symbolic_description
 ```
 
 %% Output
 
-
-    $\displaystyle \left[ grad_{x} \leftarrow \frac{{src}_{(1,0)}}{2} - \frac{{src}_{(-1,0)}}{2}, \  grad_{y} \leftarrow \frac{{src}_{(0,1)}}{2} - \frac{{src}_{(0,-1)}}{2}, \  {dst}_{(0,0)} \leftarrow grad_{x} + grad_{y}\right]$
+
+    $\displaystyle \left[ grad_{x} \leftarrow_{} \frac{{src}_{(1,0)}}{2} - \frac{{src}_{(-1,0)}}{2}, \  grad_{y} \leftarrow_{} \frac{{src}_{(0,1)}}{2} - \frac{{src}_{(0,-1)}}{2}, \  {dst}_{(0,0)} \leftarrow_{} grad_{x} + grad_{y}\right]$
     ⎡         src_E   src_W            src_N   src_S                         ⎤
     ⎢gradₓ := ───── - ─────, grad_y := ───── - ─────, dst_C := gradₓ + grad_y⎥
     ⎣           2       2                2       2                           ⎦
 
 %% Cell type:markdown id: tags:
 
 We created subexpressions, using standard sympy symbols on the left hand side, to split the kernel into multiple assignments. Defining a kernel using a list of `Assignment`s is quite tedious and hard to read.
 To simplify the formulation of a kernel, *pystencils* offers the `kernel` decorator, that transforms a normal Python function with `@=` assignments into an assignment list that can be passed to `create_kernel`.
 
 %% Cell type:code id: tags:
 
 ``` python
 @ps.kernel
 def symbolic_description_using_function():
     grad_x @= (src[1, 0] - src[-1, 0]) / 2
     grad_y @= (src[0, 1] - src[0, -1]) / 2
     dst[0, 0] @= grad_x + grad_y
 symbolic_description_using_function
 ```
 
 %% Output
 
-
-    $\displaystyle \left[ grad_{x} \leftarrow \frac{{src}_{(1,0)}}{2} - \frac{{src}_{(-1,0)}}{2}, \  grad_{y} \leftarrow \frac{{src}_{(0,1)}}{2} - \frac{{src}_{(0,-1)}}{2}, \  {dst}_{(0,0)} \leftarrow grad_{x} + grad_{y}\right]$
+
+    $\displaystyle \left[ grad_{x} \leftarrow_{} \frac{{src}_{(1,0)}}{2} - \frac{{src}_{(-1,0)}}{2}, \  grad_{y} \leftarrow_{} \frac{{src}_{(0,1)}}{2} - \frac{{src}_{(0,-1)}}{2}, \  {dst}_{(0,0)} \leftarrow_{} grad_{x} + grad_{y}\right]$
     ⎡         src_E   src_W            src_N   src_S                         ⎤
     ⎢gradₓ := ───── - ─────, grad_y := ───── - ─────, dst_C := gradₓ + grad_y⎥
     ⎣           2       2                2       2                           ⎦
 
 %% Cell type:markdown id: tags:
 
 The decorated function can contain any Python code, only the `@=` operator, and the ternary inline `if-else` operator have different meaning.
 
 ### b) Ternary 'if' with `Piecewise`
 
 The ternary operator maps to `sympy.Piecewise` functions, that can be used to introduce branching into the kernel. Piecewise defined functions must give a value for every input, i.e. there must be a 'otherwise' clause in the end that is indicated by the condition `True`. Piecewise objects are standard sympy terms that can be integrated into bigger expressions:
 
 %% Cell type:code id: tags:
 
 ``` python
 sp.Piecewise((1.0, src[0,1] > 0), (0.0, True)) + src[1, 0]
 ```
 
 %% Output
 
     $\displaystyle {src}_{(1,0)} + \begin{cases} 1.0 & \text{for}\: {src}_{(0,1)} > 0 \\0.0 & \text{otherwise} \end{cases}$
             ⎛⎧1.0  for src_N > 0⎞
     src_E + ⎜⎨                  ⎟
             ⎝⎩0.0    otherwise  ⎠
 
 %% Cell type:markdown id: tags:
 
 Piecewise objects are created by the `kernel` decorator for ternary if-else statements.
 
 %% Cell type:code id: tags:
 
 ``` python
 @ps.kernel
 def kernel_with_piecewise():
     grad_x @= (src[1, 0] - src[-1, 0]) / 2 if src[-1, 0] > 0 else 0.0
 kernel_with_piecewise
 ```
 
 %% Output
 
-
-    $\displaystyle \left[ grad_{x} \leftarrow \begin{cases} \frac{{src}_{(1,0)}}{2} - \frac{{src}_{(-1,0)}}{2} & \text{for}\: {src}_{(-1,0)} > 0 \\0.0 & \text{otherwise} \end{cases}\right]$
+
+    $\displaystyle \left[ grad_{x} \leftarrow_{} \begin{cases} \frac{{src}_{(1,0)}}{2} - \frac{{src}_{(-1,0)}}{2} & \text{for}\: {src}_{(-1,0)} > 0 \\0.0 & \text{otherwise} \end{cases}\right]$
     ⎡         ⎧src_E   src_W               ⎤
     ⎢         ⎪───── - ─────  for src_W > 0⎥
     ⎢gradₓ := ⎨  2       2                 ⎥
     ⎢         ⎪                            ⎥
     ⎣         ⎩     0.0         otherwise  ⎦
 
 %% Cell type:markdown id: tags:
 
 ### c) Assignment level optimizations using `AssignmentCollection`
 
 When the kernels get larger and more complex, it is helpful to organize the list of assignment into a more structured way. The `AssignmentCollection` offers optimizating transformation on a list of assignments. It holds two assignment lists, one for subexpressions and one for the main assignments. Main assignments are typically those that write to an array.
 
 %% Cell type:code id: tags:
 
 ``` python
 @ps.kernel
 def somewhat_longer_dummy_kernel(s):
     s.a @= src[0, 1] + src[-1, 0]
     s.b @= 2 * src[1, 0] + src[0, -1]
     s.c @= src[0, 1] + 2 * src[1, 0] + src[-1, 0] + src[0, -1] - src[0,0]
     dst[0, 0] @= s.a + s.b + s.c
 
 ac = ps.AssignmentCollection(main_assignments=somewhat_longer_dummy_kernel[-1:],
                              subexpressions=somewhat_longer_dummy_kernel[:-1])
 ac
 ```
 
 %% Output
 
-    AssignmentCollection: dst_C, <- f(src_N, src_E, src_W, src_C, src_S)
+    AssignmentCollection: dst_C, <- f(src_C, src_W, src_S, src_N, src_E)
 
 %% Cell type:code id: tags:
 
 ``` python
 ac.operation_count
 ```
 
 %% Output
 
     {'adds': 8,
      'muls': 2,
      'divs': 0,
      'sqrts': 0,
      'fast_sqrts': 0,
      'fast_inv_sqrts': 0,
      'fast_div': 0}
 
 %% Cell type:markdown id: tags:
 
 The `pystencils.simp` submodule offers several functions to optimize a collection of assignments.
 It also offers functionality to group optimization into strategies and evaluate them.
 In this example we reduce the number of operations by reusing existing subexpressions to get rid of two unnecessary floating point additions. For more information about assignment collections and simplifications see the [demo notebook](demo_assignment_collection.ipynb).
 
 %% Cell type:code id: tags:
 
 ``` python
 opt_ac = ps.simp.subexpression_substitution_in_existing_subexpressions(ac)
 opt_ac
 ```
 
 %% Output
 
-    AssignmentCollection: dst_C, <- f(src_N, src_E, src_W, src_C, src_S)
+    AssignmentCollection: dst_C, <- f(src_C, src_W, src_S, src_N, src_E)
 
 %% Cell type:code id: tags:
 
 ``` python
 opt_ac.operation_count
 ```
 
 %% Output
 
     {'adds': 6,
      'muls': 1,
      'divs': 0,
      'sqrts': 0,
      'fast_sqrts': 0,
      'fast_inv_sqrts': 0,
      'fast_div': 0}
 
 %% Cell type:markdown id: tags:
 
 ### d) Ghost layers and iteration region
 
 When creating a kernel with neighbor accesses, *pystencils* automatically restricts the iteration region, such that all accesses are safe.
 
 %% Cell type:code id: tags:
 
 ``` python
 kernel = ps.create_kernel(ps.Assignment(dst[0,0], src[2, 0] + src[-1, 0]))
 ps.show_code(kernel)
 ```
 
 %% Output
 
 
 
 %% Cell type:markdown id: tags:
 
 When no additional ghost layer information is given, *pystencils* looks at all neighboring field accesses and introduces the required number of ghost layers **for all directions**. In the example above the largest neighbor accesses was ``src[2, 0]``, so theoretically we would need 2 ghost layers only the the end of the x coordinate.
 By default *pystencils* introduces 2 ghost layers at all borders of the domain. The next cell shows how to change this behavior. Be careful with manual ghost layer specification, wrong values may lead to SEGFAULTs.
 
 %% Cell type:code id: tags:
 
 ``` python
 gl_spec = [(0, 2),   # 0 ghost layers at the left, 2 at the right border
            (1, 0)]   # 1 ghost layer at the lower y, one at the upper y coordinate
 kernel = ps.create_kernel(ps.Assignment(dst[0,0], src[2, 0] + src[-1, 0]), ghost_layers=gl_spec)
 ps.show_code(kernel)
 ```
 
 %% Output
 
 
 
 %% Cell type:markdown id: tags:
 
 ## 2 ) Restrictions
 
 
 ### a) Independence Restriction
 
 *pystencils* only works for kernels where each array element can be updated independently from all other elements. This restriction ensures that the kernels can be easily parallelized and also be run on the GPU. Trying to define kernels where the results depends on the iteration order, leads to a ValueError.
 
 %% Cell type:code id: tags:
 
 ``` python
 invalid_description = [
     ps.Assignment(dst[1, 0], src[1, 0] + src[-1, 0]),
     ps.Assignment(dst[0, 0], src[1, 0] - src[-1, 0]),
 ]
 try:
     invalid_kernel = ps.create_kernel(invalid_description)
     assert False, "Should never be executed"
 except ValueError as e:
     print(e)
 ```
 
 %% Output
 
     Field dst is written at two different locations
 
 %% Cell type:markdown id: tags:
 
 The independence restriction makes sure that the kernel can be safely parallelized by checking the following conditions: If a field is modified inside the kernel, it may only be modified at a single spatial position. In that case the field may also only be read at this position. Fields that are not modified may be read at multiple neighboring positions.
 
 Specifically, this rule allows for in-place updates that don't access neighbors.
 
 %% Cell type:code id: tags:
 
 ``` python
 valid_kernel = ps.create_kernel(ps.Assignment(src[0,0], 2*src[0,0] + 42))
 ```
 
 %% Cell type:markdown id: tags:
 
 If a field stores multiple values per cell, as in the next example, this restriction only applies for accesses with the same index.
 
 %% Cell type:code id: tags:
 
 ``` python
 v = ps.fields("v(2): double[2D]")
 valid_kernel = ps.create_kernel([ps.Assignment(v[0,0](1), 2*v[0,0](1) + 42),
-                                 ps.Assignment(v[0,1](0), 2*v[1,0](0) + 42)])
+                                 ps.Assignment(v[0,1](0), 2*v[0,1](0) + 42)])
 ```
 
 %% Cell type:markdown id: tags:
 
 ### b) Static Single Assignment Form
 
 All assignments that don't write to a field must be in SSA form
 1. Each sympy symbol may only occur once as a left-hand-side (fields can be written multiple times)
 2. A symbol has to be defined before it is used. If it is never defined it is introduced as function parameter
 
 The next cell demonstrates the first SSA restriction:
 
 %% Cell type:code id: tags:
 
 ``` python
 @ps.kernel
 def not_allowed():
     a, b = sp.symbols("a b")
     a @= src[0, 0]
     b @= a + 3
     a @= src[-1, 0]
     dst[0, 0] @= a + b
 try:
     ps.create_kernel(not_allowed)
     assert False
 except ValueError as e:
     print(e)
 ```
 
 %% Output
 
     Assignments not in SSA form, multiple assignments to a
 
 %% Cell type:markdown id: tags:
 
 Also it is not allowed to write a field at the same location
 
 %% Cell type:code id: tags:
 
 ``` python
 @ps.kernel
 def not_allowed():
     dst[0, 0] @= src[0, 1] + src[1, 0]
     dst[0, 0] @= 2 * dst[0, 0]
 
 try:
     ps.create_kernel(not_allowed)
     assert False
 except ValueError as e:
     print(e)
 ```
 
 %% Output
 
     Field dst is written twice at the same location
 
 %% Cell type:markdown id: tags:
 
 This situation should be resolved by introducing temporary variables
 
 %% Cell type:code id: tags:
 
 ``` python
 tmp_var = sp.Symbol("a")
 
 @ps.kernel
 def allowed():
     tmp_var @= src[0, 1] + src[1, 0]
     dst[0, 0] @= 2 * tmp_var
 
 
 ast = ps.create_kernel(allowed)
 ps.show_code(ast)
 ```
 
 %% Output
 
 

doc/notebooks/demo_assignment_collection.ipynb

 %% Cell type:code id: tags:
 
 ``` python
 from pystencils.session import *
 ```
 
 %% Cell type:markdown id: tags:
 
 # Demo: Assignment collections and simplification
 
 
 ## Assignment collections
 
 The assignment collection class helps to formulate and simplify assignments for numerical kernels.
 
 An ``AssignmentCollection`` is an ordered collection of assignments, together with an optional ordered collection of subexpressions, that are required to evaluate the main assignments. There are various simplification rules available that operate on ``AssignmentCollection``s.
 
 %% Cell type:markdown id: tags:
 
 We start by defining some stencil update rule. Here we also use the *pystencils* ``Field``, note however that the assignment collection module works purely on the *sympy* level.
 
 %% Cell type:code id: tags:
 
 ``` python
 a,b,c = sp.symbols("a b c")
 f = ps.fields("f(2) : [2D]")
+g = ps.fields("g(2) : [2D]")
 
-a1 = ps.Assignment(f[0,0](1), (a**2 +b) * f[0,1] + \
+a1 = ps.Assignment(g[0,0](1), (a**2 +b) * f[0,1] + \
                   (a**2 - c) * f[1,0] + \
                   (a**2 - 2*c) * f[-1,0] + \
                   (a**2) * f[0, -1])
 
-a2 = ps.Assignment(f[0,0](0), (c**2 +b) * f[0,1] + \
+a2 = ps.Assignment(g[0,0](0), (c**2 +b) * f[0,1] + \
                   (c**2 - c) * f[1,0] + \
                   (c**2 - 2*c) * f[-1,0] + \
                   (c**2 - a**2) * f[0, -1])
 
 
 ac = ps.AssignmentCollection([a1, a2], subexpressions=[])
 ac
 ```
 
 %% Output
 
-    Equation Collection for f_C^1,f_C^0
+    AssignmentCollection: g_C^0, g_C^1 <- f(f_N^0, b, f_S^0, f_E^0, a, f_W^0, c)
 
 %% Cell type:markdown id: tags:
 
 *sympy* operations can be applied on an assignment collection: In this example we first expand the collection, then look for common subexpressions.
 
 %% Cell type:code id: tags:
 
 ``` python
 expand_all = ps.simp.apply_to_all_assignments(sp.expand)
 expandedEc = expand_all(ac)
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 ac_cse = ps.simp.sympy_cse(expandedEc)
 ac_cse
 ```
 
 %% Output
 
-    Equation Collection for f_C^1,f_C^0
+    AssignmentCollection: g_C^0, g_C^1 <- f(f_N^0, b, f_S^0, f_E^0, a, f_W^0, c)
 
 %% Cell type:markdown id: tags:
 
 Symbols occuring in assignment collections are classified into 3 categories:
 - ``free_symbols``: symbols that occur in right-hand-sides but never on left-hand-sides
 - ``bound_symbols``: symbols that occur on left-hand-sides
 - ``defined_symbols``: symbols that occur on left-hand-sides of a main assignment
 
 %% Cell type:code id: tags:
 
 ``` python
 ac_cse.free_symbols
 ```
 
 %% Output
 
-
-    $$\left\{{{f}_{E}^{0}}, {{f}_{N}^{0}}, {{f}_{S}^{0}}, {{f}_{W}^{0}}, a, b, c\right\}$$
-    set([f_E__0, f_N__0, f_S__0, f_W__0, a, b, c])
+
+    $\displaystyle \left\{{f}_{(1,0)}^{0}, {f}_{(0,1)}^{0}, {f}_{(0,-1)}^{0}, {f}_{(-1,0)}^{0}, a, b, c\right\}$
+    {f_E__0, f_N__0, f_S__0, f_W__0, a, b, c}
 
 %% Cell type:code id: tags:
 
 ``` python
 ac_cse.bound_symbols
 ```
 
 %% Output
 
-
-    $$\left\{{{f}_{C}^{0}}, {{f}_{C}^{1}}, \xi_{0}, \xi_{1}, \xi_{2}, \xi_{3}\right\}$$
-    set([f_C__0, f_C__1, ξ₀, ξ₁, ξ₂, ξ₃])
+
+    $\displaystyle \left\{{g}_{(0,0)}^{0}, {g}_{(0,0)}^{1}, \xi_{0}, \xi_{1}, \xi_{2}, \xi_{3}\right\}$
+    {g_C__0, g_C__1, ξ₀, ξ₁, ξ₂, ξ₃}
 
 %% Cell type:code id: tags:
 
 ``` python
 ac_cse.defined_symbols
 ```
 
 %% Output
 
-
-    $$\left\{{{f}_{C}^{0}}, {{f}_{C}^{1}}\right\}$$
-    set([f_C__0, f_C__1])
+
+    $\displaystyle \left\{{g}_{(0,0)}^{0}, {g}_{(0,0)}^{1}\right\}$
+    {g_C__0, g_C__1}
 
 %% Cell type:markdown id: tags:
 
 Assignment collections can be splitted up, and merged together. For splitting, a list of symbols that occur on the left-hand-side in the main assignments has to be passed. The returned assignment collection only contains these main assignments together with all necessary subexpressions.
 
 %% Cell type:code id: tags:
 
 ``` python
-ac_f0 = ac_cse.new_filtered([f(0)])
-ac_f1 = ac_cse.new_filtered([f(1)])
+ac_f0 = ac_cse.new_filtered([g(0)])
+ac_f1 = ac_cse.new_filtered([g(1)])
 ac_f1
 ```
 
 %% Output
 
-    Equation Collection for f_C^1
+    AssignmentCollection: g_C^1, <- f(f_N^0, b, f_S^0, f_E^0, a, f_W^0, c)
 
 %% Cell type:markdown id: tags:
 
 Note here that $\xi_4$ is no longer part of the subexpressions, since it is not used in the main assignment of $f_C^1$.
 
 If we merge both collections together, we end up with the original collection.
 
 %% Cell type:code id: tags:
 
 ``` python
 ac_f0.new_merged(ac_f1)
 ```
 
 %% Output
 
-    Equation Collection for f_C^0,f_C^1
+    AssignmentCollection: g_C^0, g_C^1 <- f(f_N^0, b, f_S^0, f_E^0, a, f_W^0, c)
 
 %% Cell type:markdown id: tags:
 
 There is also a method that inserts all subexpressions into the main assignments. This is the inverse operation of common subexpression elimination.
 
 %% Cell type:code id: tags:
 
 ``` python
 assert sp.simplify(ac_f0.new_without_subexpressions().main_assignments[0].rhs - a2.rhs) == 0
 ac_f0.new_without_subexpressions()
 ```
 
 %% Output
 
-    Equation Collection for f_C^0
+    AssignmentCollection: g_C^0, <- f(f_N^0, b, f_S^0, f_E^0, a, f_W^0, c)
 
 %% Cell type:markdown id: tags:
 
 To evaluate an assignment collection, use the ``lambdify`` method. It is very similar to *sympy*s ``lambdify`` function.
 
 %% Cell type:code id: tags:
 
 ``` python
 evalFct = ac_cse.lambdify([f[0,1], f[1,0]],  # new parameters of returned function
                           fixed_symbols={a:1, b:2, c:3, f[0,-1]: 4, f[-1,0]: 5}) # fix values of other symbols
 evalFct(2,1)
 ```
 
 %% Output
 
-
-    $$\left \{ {{f}_{C}^{0}} : 75, \quad {{f}_{C}^{1}} : -17\right \}$$
-    {f_C__0: 75, f_C__1: -17}
+
+    $\displaystyle \left\{ {g}_{(0,0)}^{0} : 75, \  {g}_{(0,0)}^{1} : -17\right\}$
+    {g_C__0: 75, g_C__1: -17}
 
 %% Cell type:markdown id: tags:
 
 lambdify is rather slow for evaluation. The intended way to evaluate an assignment collection is *pystencils* i.e. create a fast kernel, that applies the update at every site of a structured grid. The collection can be directly passed to the `create_kernel` function.
 
 %% Cell type:code id: tags:
 
 ``` python
 func = ps.create_kernel(ac_cse).compile()
 ```
 
 %% Cell type:markdown id: tags:
 
 ## Simplification Strategies
 
 In above examples, we already applied simplification rules to assignment collections. Simplification rules are functions that take, as a single argument, an assignment collection and return an modified/simplified copy of it. The ``SimplificationStrategy`` class holds a list of simplification rules and can apply all of them in the specified order. Additionally it provides useful printing and reporting functions.
 
 We start by creating a simplification strategy, consisting of the expand and CSE simplifications we have already applied above:
 
 %% Cell type:code id: tags:
 
 ``` python
 strategy = ps.simp.SimplificationStrategy()
 strategy.add(ps.simp.apply_to_all_assignments(sp.expand))
 strategy.add(ps.simp.sympy_cse)
 ```
 
 %% Cell type:markdown id: tags:
 
 This strategy can be applied to any assignment collection:
 
 %% Cell type:code id: tags:
 
 ``` python
 strategy(ac)
 ```
 
 %% Output
 
-    Equation Collection for f_C^1,f_C^0
+    AssignmentCollection: g_C^0, g_C^1 <- f(f_N^0, b, f_S^0, f_E^0, a, f_W^0, c)
 
 %% Cell type:markdown id: tags:
 
 The strategy can also print the simplification results at each stage.
 The report contains information about the number of operations after each simplification as well as the runtime of each simplification routine.
 
 %% Cell type:code id: tags:
 
 ``` python
 strategy.create_simplification_report(ac)
 ```
 
 %% Output
 
-    <pystencils.simp.simplificationstrategy.SimplificationStrategy.create_simplification_report.<locals>.Report at 0x7f9be404fda0>
+    <pystencils.simp.simplificationstrategy.SimplificationStrategy.create_simplification_report.<locals>.Report at 0x147de3e90>
 
 %% Cell type:markdown id: tags:
 
 The strategy can also print the full collection after each simplification...
 
 %% Cell type:code id: tags:
 
 ``` python
 strategy.show_intermediate_results(ac)
 ```
 
 %% Output
 
-    <pystencils.simp.simplificationstrategy.SimplificationStrategy.show_intermediate_results.<locals>.IntermediateResults at 0x7f9bad688dd8>
+    <pystencils.simp.simplificationstrategy.SimplificationStrategy.show_intermediate_results.<locals>.IntermediateResults at 0x147e09c90>
 
 %% Cell type:markdown id: tags:
 
 ... or only specific assignments for better readability
 
 %% Cell type:code id: tags:
 
 ``` python
-strategy.show_intermediate_results(ac, symbols=[f(1)])
+strategy.show_intermediate_results(ac, symbols=[g(1)])
 ```
 
 %% Output
 
-    <pystencils.simp.simplificationstrategy.SimplificationStrategy.show_intermediate_results.<locals>.IntermediateResults at 0x7f9bad688b00>
+    <pystencils.simp.simplificationstrategy.SimplificationStrategy.show_intermediate_results.<locals>.IntermediateResults at 0x1265a1b90>
+
+%% Cell type:code id: tags:
+
+``` python
+```

pyproject.toml

@@ -12,7 +12,7 @@ authors = [
 ]
 license = { file = "COPYING.txt" }
 requires-python = ">=3.10"
-dependencies = ["sympy>=1.6,<=1.11.1", "numpy>=1.8.0", "appdirs", "joblib", "pyyaml"]
+dependencies = ["sympy>=1.9,<=1.12.1", "numpy>=1.8.0", "appdirs", "joblib", "pyyaml"]
 classifiers = [
     "Development Status :: 4 - Beta",
     "Framework :: Jupyter",

src/pystencils/backends/arm_instruction_sets.py

@@ -18,8 +18,11 @@ def get_argument_string(function_shortcut, first=''):
 def get_vector_instruction_set_arm(data_type='double', instruction_set='neon'):
     if instruction_set not in ['neon', 'sme'] and not instruction_set.startswith('sve'):
         raise NotImplementedError(instruction_set)
-    if instruction_set in ['sve', 'sme']:
+    if instruction_set in ['sve', 'sve2', 'sme']:
         cmp = 'cmp'
+    elif instruction_set.startswith('sve2') and instruction_set not in ('sve256', 'sve2048'):
+        cmp = 'cmp'
+        bitwidth = int(instruction_set[4:])
     elif instruction_set.startswith('sve'):
         cmp = 'cmp'
         bitwidth = int(instruction_set[3:])
@@ -52,7 +55,7 @@ def get_vector_instruction_set_arm(data_type='double', instruction_set='neon'):
 
     result = dict()
 
-    if instruction_set in ['sve', 'sme']:
+    if instruction_set in ['sve', 'sve2', 'sme']:
         width = 'svcntd()' if data_type == 'double' else 'svcntw()'
         intwidth = 'svcntw()'
         result['bytes'] = 'svcntb()'
@@ -61,13 +64,14 @@ def get_vector_instruction_set_arm(data_type='double', instruction_set='neon'):
         intwidth = bitwidth // bits['int']
         result['bytes'] = bitwidth // 8
     if instruction_set.startswith('sve') or instruction_set == 'sme':
+        base_names['stream'] = 'stnt1[0, 1]'
         prefix = 'sv'
         suffix = f'_f{bits[data_type]}' 
     elif instruction_set == 'neon':
         prefix = 'v'
         suffix = f'q_f{bits[data_type]}' 
 
-    if instruction_set in ['sve', 'sme']:
+    if instruction_set in ['sve', 'sve2', 'sme']:
         predicate = f'{prefix}whilelt_b{bits[data_type]}_u64({{loop_counter}}, {{loop_stop}})'
         int_predicate = f'{prefix}whilelt_b{bits["int"]}_u64({{loop_counter}}, {{loop_stop}})'
     else:
@@ -86,7 +90,7 @@ def get_vector_instruction_set_arm(data_type='double', instruction_set='neon'):
 
         result[intrinsic_id] = prefix + name + suffix + undef + arg_string
 
-    if instruction_set in ['sve', 'sme']:
+    if instruction_set in ['sve', 'sve2', 'sme']:
         from pystencils.backends.cbackend import CFunction
         result['width'] = CFunction(width, "int")
         result['intwidth'] = CFunction(intwidth, "int")
@@ -105,15 +109,18 @@ def get_vector_instruction_set_arm(data_type='double', instruction_set='neon'):
                                vindex.format("{2}") + ', {1})'
             result['loadS'] = f'svld1_gather_u{bits[data_type]}index_f{bits[data_type]}({predicate}, {{0}}, ' + \
                               vindex.format("{1}") + ')'
+        if instruction_set.startswith('sve2') and instruction_set not in ('sve256', 'sve2048'):
+            result['streamS'] = f'svstnt1_scatter_u{bits[data_type]}offset_f{bits[data_type]}({predicate}, {{0}}, ' + \
+                                vindex.format(f"{{2}}*{bits[data_type]//8}") + ', {1})'
 
         result['+int'] = f"svadd_s{bits['int']}_x({int_predicate}, " + "{0}, {1})"
 
-        result['float'] = f'svfloat{bits["float"]}_{"s" if instruction_set not in ["sve", "sme"] else ""}t'
-        result['double'] = f'svfloat{bits["double"]}_{"s" if instruction_set not in ["sve", "sme"] else ""}t'
-        result['int'] = f'svint{bits["int"]}_{"s" if instruction_set not in ["sve", "sme"] else ""}t'
-        result['bool'] = f'svbool_{"s" if instruction_set not in ["sve", "sme"] else ""}t'
+        result['float'] = f'svfloat{bits["float"]}_{"s" if instruction_set not in ["sve", "sve2", "sme"] else ""}t'
+        result['double'] = f'svfloat{bits["double"]}_{"s" if instruction_set not in ["sve", "sve2", "sme"] else ""}t'
+        result['int'] = f'svint{bits["int"]}_{"s" if instruction_set not in ["sve", "sve2", "sme"] else ""}t'
+        result['bool'] = f'svbool_{"s" if instruction_set not in ["sve", "sve2", "sme"] else ""}t'
 
-        result['headers'] = ['<arm_sve.h>', '"arm_neon_helpers.h"']
+        result['headers'] = ['<arm_sve.h>', '<arm_acle.h>', '"arm_neon_helpers.h"']
 
         result['&'] = f'svand_b_z({predicate},' + ' {0}, {1})'
         result['|'] = f'svorr_b_z({predicate},' + ' {0}, {1})'
@@ -122,12 +129,17 @@ def get_vector_instruction_set_arm(data_type='double', instruction_set='neon'):
         result['all'] = f'svcntp_b{bits[data_type]}({predicate}, {{0}}) == {width}'
 
         result['maskStoreU'] = result['storeU'].replace(predicate, '{2}')
+        result['maskStream'] = result['stream'].replace(predicate, '{2}')
         if instruction_set != 'sme':
             result['maskStoreS'] = result['storeS'].replace(predicate, '{3}')
+            if instruction_set.startswith('sve2') and instruction_set not in ('sve256', 'sve2048'):
+                result['maskStreamS'] = result['streamS'].replace(predicate, '{3}')
+        
+        result['streamFence'] = '__dmb(15)'
 
         if instruction_set == 'sme':
             result['function_prefix'] = '__attribute__((arm_locally_streaming))'
-        elif instruction_set not in ['sve', 'sme']:
+        elif instruction_set not in ['sve', 'sve2', 'sme']:
             result['compile_flags'] = [f'-msve-vector-bits={bitwidth}']
     else:
         result['makeVecConst'] = f'vdupq_n_f{bits[data_type]}' + '({0})'
@@ -152,7 +164,9 @@ def get_vector_instruction_set_arm(data_type='double', instruction_set='neon'):
         result['any'] = f'vaddlvq_u8(vreinterpretq_u8_u{bits[data_type]}({{0}})) > 0'
         result['all'] = f'vaddlvq_u8(vreinterpretq_u8_u{bits[data_type]}({{0}})) == 16*0xff'
 
+        # SVE has real nontemporal stores, so we only need to zero cachlines on Neon
+        result['cachelineZero'] = 'cachelineZero((void*) {0})'
+
     result['cachelineSize'] = 'cachelineSize()'
-    result['cachelineZero'] = 'cachelineZero((void*) {0})'
 
     return result

src/pystencils/backends/cbackend.py

@@ -280,14 +280,25 @@ class CBackend:
             if type(lhs_type) is VectorType and isinstance(node.lhs, CastFunc):
                 arg, data_type, aligned, nontemporal, mask, stride = node.lhs.args
                 instr = 'storeU'
-                if aligned:
+                if nontemporal and 'storeA' not in self._vector_instruction_set and \
+                        'stream' in self._vector_instruction_set:
+                    instr = 'stream'
+                elif aligned:
                     instr = 'stream' if nontemporal and 'stream' in self._vector_instruction_set else 'storeA'
                 if mask != True:  # NOQA
-                    instr = 'maskStoreA' if aligned else 'maskStoreU'
+                    instr = 'maskStream' if nontemporal and 'maskStream' in self._vector_instruction_set else \
+                            'maskStoreA' if aligned else 'maskStoreU'
                     if instr not in self._vector_instruction_set:
-                        self._vector_instruction_set[instr] = self._vector_instruction_set['store' + instr[-1]].format(
+                        if instr == 'maskStream' and 'stream' in self._vector_instruction_set:
+                            store, load = 'stream', 'loadA'
+                        elif (instr in ('maskStream', 'maskStoreA')) and 'storeA' in self._vector_instruction_set:
+                            store, load = 'storeA', 'loadA'
+                        else:
+                            store, load = 'storeU', 'loadU'
+                        load = load if load in self._vector_instruction_set else 'loadU'
+                        self._vector_instruction_set[instr] = self._vector_instruction_set[store].format(
                             '{0}', self._vector_instruction_set['blendv'].format(
-                                self._vector_instruction_set['load' + instr[-1]].format('{0}', **self._kwargs),
+                                self._vector_instruction_set[load].format('{0}', **self._kwargs),
                                 '{1}', '{2}', **self._kwargs), **self._kwargs)
                     printed_mask = self.sympy_printer.doprint(mask)
                     if data_type.base_type.c_name == 'double':
@@ -312,12 +323,14 @@ class CBackend:
                 ptr = "&" + self.sympy_printer.doprint(node.lhs.args[0])
 
                 if stride != 1:
-                    instr = 'maskStoreS' if mask != True else 'storeS'  # NOQA
+                    instr = ('maskStreamS' if nontemporal and 'maskStreamS' in self._vector_instruction_set else
+                             'maskStoreS') if mask != True else \
+                            ('streamS' if nontemporal and 'streamS' in self._vector_instruction_set else 'storeS')  # NOQA
                     return self._vector_instruction_set[instr].format(ptr, self.sympy_printer.doprint(rhs),
                                                                       stride, printed_mask, **self._kwargs) + ';'
 
                 pre_code = ''
-                if nontemporal and 'cachelineZero' in self._vector_instruction_set:
+                if nontemporal and 'cachelineZero' in self._vector_instruction_set and mask == True:  # NOQA
                     first_cond = f"((uintptr_t) {ptr} & {CachelineSize.mask_symbol}) == 0"
                     offset = sp.Add(*[sp.Symbol(LoopOverCoordinate.get_loop_counter_name(i))
                                       * node.lhs.args[0].field.spatial_strides[i] for i in
@@ -337,15 +350,22 @@ class CBackend:
                     code2 = self._vector_instruction_set['flushCacheline'].format(
                         ptr, self.sympy_printer.doprint(rhs), **self._kwargs) + ';'
                     code = f"{code}\nif ({flushcond}) {{\n\t{code2}\n}}"
-                elif nontemporal and 'storeAAndFlushCacheline' in self._vector_instruction_set:
+                elif aligned and nontemporal and 'storeAAndFlushCacheline' in self._vector_instruction_set:
                     lhs_hash = hashlib.sha1(self.sympy_printer.doprint(node.lhs).encode('ascii')).hexdigest()[:8]
                     rhs_hash = hashlib.sha1(self.sympy_printer.doprint(rhs).encode('ascii')).hexdigest()[:8]
                     tmpvar = f'_tmp_{lhs_hash}_{rhs_hash}'
                     code = 'const ' + self._print(node.lhs.dtype).replace(' const', '') + ' ' + tmpvar + ' = ' \
                         + self.sympy_printer.doprint(rhs) + ';'
                     code1 = self._vector_instruction_set[instr].format(ptr, tmpvar, printed_mask, **self._kwargs) + ';'
-                    code2 = self._vector_instruction_set['storeAAndFlushCacheline'].format(ptr, tmpvar, printed_mask,
-                                                                                           **self._kwargs) + ';'
+                    maskStore, store, load = 'maskStoreAAndFlushCacheline', 'storeAAndFlushCacheline', 'loadA'
+                    instr2 = maskStore if mask != True else store  # NOQA
+                    if instr2 not in self._vector_instruction_set:
+                        self._vector_instruction_set[maskStore] = self._vector_instruction_set[store].format(
+                            '{0}', self._vector_instruction_set['blendv'].format(
+                                self._vector_instruction_set[load].format('{0}', **self._kwargs),
+                                '{1}', '{2}', **self._kwargs),
+                            **self._kwargs)
+                    code2 = self._vector_instruction_set[instr2].format(ptr, tmpvar, printed_mask, **self._kwargs) + ';'
                     code += f"\nif ({flushcond}) {{\n\t{code2}\n}} else {{\n\t{code1}\n}}"
                 return pre_code + code
             else:

src/pystencils/backends/riscv_instruction_sets.py

@@ -34,7 +34,7 @@ def get_vector_instruction_set_riscv(data_type='double', instruction_set='rvv'):
         'maskStoreU': f'se{bits[data_type]}_v[2, 0, 1]',
         'loadS': f'lse{bits[data_type]}_v[0, 1]',
         'storeS': f'sse{bits[data_type]}_v[0, 2, 1]',
-        'maskStoreS': f'sse{bits[data_type]}_v[2, 0, 3, 1]',
+        'maskStoreS': f'sse{bits[data_type]}_v[3, 0, 2, 1]',
 
         'abs': 'fabs_v[0]',
         '==': 'mfeq_vv[0, 1]',
@@ -89,7 +89,7 @@ def get_vector_instruction_set_riscv(data_type='double', instruction_set='rvv'):
 
     result['storeS'] = result['storeS'].replace('{2}', f'{{2}}*{bits[data_type]//8}')
     result['loadS'] = result['loadS'].replace('{1}', f'{{1}}*{bits[data_type]//8}')
-    result['maskStoreS'] = result['maskStoreS'].replace('{3}', f'{{3}}*{bits[data_type]//8}')
+    result['maskStoreS'] = result['maskStoreS'].replace('{2}', f'{{2}}*{bits[data_type]//8}')
 
     result['+int'] = f"vadd_vv_i{bits['int']}m1({{0}}, {{1}}, {int_vl})"
 

src/pystencils/backends/simd_instruction_sets.py

@@ -45,15 +45,19 @@ def get_supported_instruction_sets():
         result = ['neon']  # Neon is mandatory on 64-bit ARM
         libc = CDLL('libc.so.6')
         hwcap = libc.getauxval(16)  # AT_HWCAP
+        hwcap2 = libc.getauxval(26)  # AT_HWCAP2
         if hwcap & (1 << 22):  # HWCAP_SVE
+            if hwcap2 & (1 << 1):  # HWCAP2_SVE2
+                name = 'sve2'
+            else:
+                name = 'sve'
             length = 8 * libc.prctl(51, 0, 0, 0, 0)  # PR_SVE_GET_VL
             if length < 0:
                 raise OSError("SVE length query failed")
             while length >= 128:
-                result.append(f"sve{length}")
+                result.append(f"{name}{length}")
                 length //= 2
-            result.append("sve")
-        hwcap2 = libc.getauxval(26)  # AT_HWCAP2
+            result.append(name)
         if hwcap2 & (1 << 23):  # HWCAP2_SME
             result.append("sme")
         return result

src/pystencils/boundaries/boundaryhandling.py

@@ -35,11 +35,11 @@ class FlagInterface:
         >>> dh = create_data_handling((4, 5))
         >>> fi = FlagInterface(dh, 'flag_field', np.uint8)
         >>> assert dh.has_data('flag_field')
-        >>> fi.reserve_next_flag()
+        >>> int(fi.reserve_next_flag())
         2
-        >>> fi.reserve_flag(4)
+        >>> int(fi.reserve_flag(4))
         4
-        >>> fi.reserve_next_flag()
+        >>> int(fi.reserve_next_flag())
         8
     """
 
@@ -450,5 +450,6 @@ def create_boundary_kernel(field, index_field, stencil, boundary_functor, target
     dir_symbol = TypedSymbol("dir", np.int32)
     elements += [SympyAssignment(dir_symbol, index_field[0]('dir'))]
     elements += boundary_functor(field, direction_symbol=dir_symbol, index_field=index_field)
-    config = CreateKernelConfig(index_fields=[index_field], target=target, **kernel_creation_args)
+    config = CreateKernelConfig(index_fields=[index_field], target=target, skip_independence_check=True,
+                                **kernel_creation_args)
     return create_kernel(elements, config=config)

src/pystencils/config.py

@@ -135,8 +135,9 @@ class CreateKernelConfig:
     """
     skip_independence_check: bool = False
     """
-    Don't check that loop iterations are independent. This is needed e.g. for
-    periodicity kernel, that access the field outside the iteration bounds. Use with care!
+    By default the assignment list is checked for read/write independence. This means fields are only written at
+    locations where they are read. Doing so guarantees thread safety. In some cases e.g. for
+    periodicity kernel, this can not be assured and does the check needs to be deactivated. Use with care!
     """
 
     class DataTypeFactory:

src/pystencils/gpu/indexing.py

@@ -224,6 +224,9 @@ class BlockIndexing(AbstractIndexing):
             assert len(self._iteration_space) == len(arr_shape), "Iteration space must be equal to the array shape"
             numeric_iteration_slice = _get_numeric_iteration_slice(self._iteration_space, arr_shape)
         end = [s.stop if s.stop != 0 else 1 for s in numeric_iteration_slice]
+        for i, s in enumerate(numeric_iteration_slice):
+            if s.step and s.step != 1:
+                end[i] = div_ceil(s.stop - s.start, s.step) + s.start
 
         if self._dim < 4:
             conditions = [c < e for c, e in zip(self.coordinates, end)]

src/pystencils/include/PyStencilsField.h

-#pragma once
-
-extern "C++" {
-#ifdef __CUDA_ARCH__
-template <typename DTYPE_T, std::size_t DIMENSION> struct PyStencilsField {
-  DTYPE_T *data;
-  DTYPE_T shape[DIMENSION];
-  DTYPE_T stride[DIMENSION];
-};
-#else
-#include <array>
-
-template <typename DTYPE_T, std::size_t DIMENSION> struct PyStencilsField {
-  DTYPE_T *data;
-  std::array<DTYPE_T, DIMENSION> shape;
-  std::array<DTYPE_T, DIMENSION> stride;
-};
-#endif
-}

src/pystencils/include/aesni_rand.h

+/*
+Copyright 2010-2011, D. E. Shaw Research. All rights reserved.
+Copyright 2019-2023, Michael Kuron.
+
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions are
+met:
+
+* Redistributions of source code must retain the above copyright
+  notice, this list of conditions, and the following disclaimer.
+
+* Redistributions in binary form must reproduce the above copyright
+  notice, this list of conditions, and the following disclaimer in the
+  documentation and/or other materials provided with the distribution.
+
+* Neither the name of of the copyright holder nor the names of its
+  contributors may be used to endorse or promote products derived from
+  this software without specific prior written permission.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
+OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
+SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
+LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
+DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
+THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+*/
+
 #include <emmintrin.h> // SSE2
 #include <wmmintrin.h> // AES
 #ifdef __AVX__

src/pystencils/include/arm_neon_helpers.h

+/*
+Copyright 2021-2023, Michael Kuron.
+
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions are
+met:
+
+* Redistributions of source code must retain the above copyright
+  notice, this list of conditions, and the following disclaimer.
+
+* Redistributions in binary form must reproduce the above copyright
+  notice, this list of conditions, and the following disclaimer in the
+  documentation and/or other materials provided with the distribution.
+
+* Neither the name of of the copyright holder nor the names of its
+  contributors may be used to endorse or promote products derived from
+  this software without specific prior written permission.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
+OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
+SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
+LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
+DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
+THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+*/
+
 #if defined(_MSC_VER)
 #define __ARM_NEON
 #endif

src/pystencils/include/gpu_defines.h

+/*
+Copyright 2023, Markus Holzer.
+Copyright 2023, Michael Kuron.
+
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions are
+met:
+
+* Redistributions of source code must retain the above copyright
+  notice, this list of conditions, and the following disclaimer.
+
+* Redistributions in binary form must reproduce the above copyright
+  notice, this list of conditions, and the following disclaimer in the
+  documentation and/or other materials provided with the distribution.
+
+* Neither the name of of the copyright holder nor the names of its
+  contributors may be used to endorse or promote products derived from
+  this software without specific prior written permission.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
+OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
+SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
+LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
+DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
+THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+*/
+
 #pragma once
 
 #define POS_INFINITY __int_as_float(0x7f800000)

src/pystencils/include/half_precision.h

+/*
+Copyright 2023, Markus Holzer.
+
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions are
+met:
+
+* Redistributions of source code must retain the above copyright
+  notice, this list of conditions, and the following disclaimer.
+
+* Redistributions in binary form must reproduce the above copyright
+  notice, this list of conditions, and the following disclaimer in the
+  documentation and/or other materials provided with the distribution.
+
+* Neither the name of of the copyright holder nor the names of its
+  contributors may be used to endorse or promote products derived from
+  this software without specific prior written permission.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
+OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
+SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
+LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
+DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
+THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+*/
+
+
 /// Half precision support. Experimental. Use carefully.
 ///
 /// This feature is experimental, since it strictly depends on the underlying architecture and compiler support.

src/pystencils/include/myintrin.h

+/*
+Copyright 2019-2023, Michael Kuron.
+
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions are
+met:
+
+* Redistributions of source code must retain the above copyright
+  notice, this list of conditions, and the following disclaimer.
+
+* Redistributions in binary form must reproduce the above copyright
+  notice, this list of conditions, and the following disclaimer in the
+  documentation and/or other materials provided with the distribution.
+
+* Neither the name of of the copyright holder nor the names of its
+  contributors may be used to endorse or promote products derived from
+  this software without specific prior written permission.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
+OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
+SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
+LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
+DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
+THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+*/
+
 #pragma once
 
 #if defined(__SSE2__) || (defined(_MSC_VER) && !defined(_M_ARM64))

src/pystencils/include/philox_rand.h

+/*
+Copyright 2010-2011, D. E. Shaw Research. All rights reserved.
+Copyright 2019-2024, Michael Kuron.
+
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions are
+met:
+
+* Redistributions of source code must retain the above copyright
+  notice, this list of conditions, and the following disclaimer.
+
+* Redistributions in binary form must reproduce the above copyright
+  notice, this list of conditions, and the following disclaimer in the
+  documentation and/or other materials provided with the distribution.
+
+* Neither the name of of the copyright holder nor the names of its
+  contributors may be used to endorse or promote products derived from
+  this software without specific prior written permission.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
+OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
+SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
+LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
+DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
+THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+*/
+
 #if !defined(__OPENCL_VERSION__) && !defined(__HIPCC_RTC__)
 #if defined(__SSE2__) || (defined(_MSC_VER) && !defined(_M_ARM64))
 #include <emmintrin.h> // SSE2

src/pystencils/include/ppc_altivec_helpers.h

+/*
+Copyright 2021, Michael Kuron.
+
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions are
+met:
+
+* Redistributions of source code must retain the above copyright
+  notice, this list of conditions, and the following disclaimer.
+
+* Redistributions in binary form must reproduce the above copyright
+  notice, this list of conditions, and the following disclaimer in the
+  documentation and/or other materials provided with the distribution.
+
+* Neither the name of of the copyright holder nor the names of its
+  contributors may be used to endorse or promote products derived from
+  this software without specific prior written permission.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
+OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
+SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
+LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
+DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
+THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+*/
+
 #include <altivec.h>
 #undef vector
 #undef bool