*pystencils* is a package that can speed up computations on *numpy* arrays. All computations are carried out fully parallel on CPUs (single node with OpenMP, multiple nodes with MPI) or on GPUs.
*pystencils* is a package that can speed up computations on *numpy* arrays. All computations are carried out fully parallel on CPUs (single node with OpenMP, multiple nodes with MPI) or on GPUs.
It is suited for applications that run the same operation on *numpy* arrays multiple times. It can be used to accelerate computations on images or voxel fields. Its main application, however, are numerical simulations using finite differences, finite volumes, or lattice Boltzmann methods.
It is suited for applications that run the same operation on *numpy* arrays multiple times. It can be used to accelerate computations on images or voxel fields. Its main application, however, are numerical simulations using finite differences, finite volumes, or lattice Boltzmann methods.
There already exist a variety of packages to speed up numeric Python code. One could use pure numpy or solutions that compile your code, like *Cython* and *numba*. See [this page](demo_benchmark.ipynb) for a comparison of these tools.
There already exist a variety of packages to speed up numeric Python code. One could use pure numpy or solutions that compile your code, like *Cython* and *numba*. See [this page](demo_benchmark.ipynb) for a comparison of these tools.
As the name suggests, *pystencils* was developed for **stencil codes**, i.e. operations that update array elements using only a local neighborhood.
As the name suggests, *pystencils* was developed for **stencil codes**, i.e. operations that update array elements using only a local neighborhood.
It generates C code, compiles it behind the scenes, and lets you call the compiled C function as if it was a native Python function.
It generates C code, compiles it behind the scenes, and lets you call the compiled C function as if it was a native Python function.
But lets not dive too deep into the concepts of *pystencils* here, they are covered in detail in the following tutorials. Let's instead look at a simple example, that computes the average neighbor values of a *numpy* array. Therefor we first create two rather large arrays for input and output:
But lets not dive too deep into the concepts of *pystencils* here, they are covered in detail in the following tutorials. Let's instead look at a simple example, that computes the average neighbor values of a *numpy* array. Therefor we first create two rather large arrays for input and output:
%% Cell type:code id: tags:
%% Cell type:code id: tags:
``` python
``` python
input_arr=np.random.rand(1024,1024)
input_arr=np.random.rand(1024,1024)
output_arr=np.zeros_like(input_arr)
output_arr=np.zeros_like(input_arr)
```
```
%% Cell type:markdown id: tags:
%% Cell type:markdown id: tags:
We first implement a version of this algorithm using pure numpy and benchmark it.
We first implement a version of this algorithm using pure numpy and benchmark it.
Here we first have created a symbolic notation of the stencil itself. This representation is built on top of *sympy* and is explained in detail in the next section.
Here we first have created a symbolic notation of the stencil itself. This representation is built on top of *sympy* and is explained in detail in the next section.
This description is then compiled and loaded as a Python function.
This description is then compiled and loaded as a Python function.
This whole process might seem overly complicated. We have already spent more lines of code than we needed for the *numpy* implementation and don't have anything running yet! However, this multi-stage process of formulating the algorithm symbolically, and just in the end actually running it, is what makes *pystencils* faster and more flexible than other approaches.
This whole process might seem overly complicated. We have already spent more lines of code than we needed for the *numpy* implementation and don't have anything running yet! However, this multi-stage process of formulating the algorithm symbolically, and just in the end actually running it, is what makes *pystencils* faster and more flexible than other approaches.
Now finally lets benchmark the *pystencils* kernel.
Now finally lets benchmark the *pystencils* kernel.
%% Cell type:code id: tags:
%% Cell type:code id: tags:
``` python
``` python
defpystencils_kernel():
defpystencils_kernel():
kernel(src=input_arr,dst=output_arr)
kernel(src=input_arr,dst=output_arr)
```
```
%% Cell type:code id: tags:
%% Cell type:code id: tags:
``` python
``` python
%%timeit
%%timeit
pystencils_kernel()
pystencils_kernel()
```
```
%%%% Output: stream
%%%% Output: stream
639 µs ± 35 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
643 µs ± 8.66 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
%% Cell type:markdown id: tags:
%% Cell type:markdown id: tags:
This benchmark shows that *pystencils* is a lot faster than pure *numpy*, especially for large arrays.
This benchmark shows that *pystencils* is a lot faster than pure *numpy*, especially for large arrays.
If you are interested in performance details and comparison to other packages like Cython, have a look at [this page](demo_benchmark.ipynb).
If you are interested in performance details and comparison to other packages like Cython, have a look at [this page](demo_benchmark.ipynb).
%% Cell type:markdown id: tags:
%% Cell type:markdown id: tags:
## Short *sympy* introduction
## Short *sympy* introduction
In this tutorial we continue with a short *sympy* introduction, since the symbolic kernel definition is built on top of this package. If you already know *sympy* you can skip this section.
In this tutorial we continue with a short *sympy* introduction, since the symbolic kernel definition is built on top of this package. If you already know *sympy* you can skip this section.
You can also read the full [sympy documentation here](http://docs.sympy.org/latest/index.html).
You can also read the full [sympy documentation here](http://docs.sympy.org/latest/index.html).
%% Cell type:code id: tags:
%% Cell type:code id: tags:
``` python
``` python
importsympyassp
importsympyassp
sp.init_printing()# enable nice LaTeX output
sp.init_printing()# enable nice LaTeX output
```
```
%% Cell type:markdown id: tags:
%% Cell type:markdown id: tags:
*sympy* is a package for symbolic calculation. So first we need some symbols:
*sympy* is a package for symbolic calculation. So first we need some symbols:
%% Cell type:code id: tags:
%% Cell type:code id: tags:
``` python
``` python
x=sp.Symbol("x")
x=sp.Symbol("x")
y=sp.Symbol("y")
y=sp.Symbol("y")
type(x)
type(x)
```
```
%%%% Output: execute_result
%%%% Output: execute_result
sympy.core.symbol.Symbol
sympy.core.symbol.Symbol
%% Cell type:markdown id: tags:
%% Cell type:markdown id: tags:
The usual mathematical operations are defined for symbols:
The usual mathematical operations are defined for symbols:
