# Using Futhark with PyOpenCL

*Code updated in May 2019.*

Python is a language with many qualities, but few would claim that
performance is among them. While libraries such as Numpy can be used,
they are not as flexible as being able to write code directly in a
high-performance language. Unfortunately, writing the
performance-critical parts of a Python program in (say) C is not
always a good experience, and the interfacing between the Python code
and the C code can be awkward and inelegant (although to be fair, it
is still nicer in Python than in many other languages). It would be
more convenient if we could compile a high-performance language
directly to a Python module that we could then `import`

like any
other piece of Python code. Of course, this entire exercise is only
worthwhile if the code in the resulting Python module executes much
faster than manually written Python. Fortunately, when most of the
computation can be offloaded to the GPU via OpenCL, the Futhark
compiler is capable of this feat.

OpenCL works by having an ordinary program running on the CPU that
transmits code and data to the GPU (or any other *accelerator*, but
we’ll stick to GPUs). In the ideal case, the CPU-code is mostly glue
that performs bookkeeping and making API calls - in other words, not
resource-intensive, and exactly what Python is good at. No matter the
language the CPU code is written in, the GPU code will be written in
OpenCL C and translated at program initialisation to whatever machine
code is needed by the concrete GPU.

This is what is exploited by the PyOpenCL code generation backend in the Futhark compiler. Certainly, the CPU-level code is written in pure Python and quite slow, but all it does is use the PyOpenCL library to offload work to the GPU. The fact that this offloading takes place is hidden from the user of the generated code, who is provided a module with functions that accept and produce ordinary Numpy arrays.

To demonstrate this capability, we will implement a simple Mandelbrot
set visualisation
program in Python. We will use Futhark to compute the set and produce
the pixel data, and make use of Pythons png module to write an image
file with the result. Our program will consist of two files:
`mandelbrot.fut`

and `mandelbrot-run.py`

.

## mandelbrot.fut

(Do not worry too much about the specifics of how the Mandelbrot set is defined - it is not the point of this post.)

Since Futhark does not have built-in support for complex numbers, we have to define our own. Futhark does not yet support proper user-defined types, so we decide to simply represent complex numbers as pairs of `f32`s. We need three operations: dot product, multiplication, and addition:

```
let dot (r: f32, i: f32): f32 =
r * r + i * i
let multComplex (a: f32, b: f32) (c: f32, d: f32): (f32,f32) =
(a*c - b * d,
a*d + b * c)
let addComplex (a: f32, b: f32) (c: f32, d: f32): (f32,f32) =
(a + c,
b + d)
```

We can now define the core function that determines whether a given
point on the complex plane is part of the Mandelbrot set. We do this
by defining a function `divergence`

that returns the iteration at
which the loop diverges (or the limit, `depth`

, if it does not):

```
let divergence (depth: i32) (c0: (f32,f32)): i32 =
let (_, i) =
loop (c, i) = (c0, 0) while (i < depth) && (dot c < 4.0) do
(addComplex c0 (multComplex c c), i + 1)
in i
```

The `mandelbrot`

function returns the divergence point for the
complex number corresponding to a pixel in a given view of the complex
plane:

```
let mandelbrot (screenX: i32) (screenY: i32) (depth: i32)
(xmin: f32) (ymin: f32) (xmax: f32) (ymax: f32): [screenY][screenX]i32 =
let sizex = xmax - xmin
let sizey = ymax - ymin
in map (\(y: i32): [screenX]i32 ->
map (\(x: i32): i32 ->
let c0 = (xmin + (r32 x * sizex) / r32 screenX,
ymin + (r32 y * sizey) / r32 screenY)
in (divergence depth c0))
(iota screenX))
(iota screenY)
```

Given the point of divergence for a pixel, we can decide on a colour, which is encoded as RGB within a 32-bit integer (the alpha channel is not used):

```
let escapeToColour (depth: i32) (divergence: i32): i32 =
if depth == divergence
then 0
else let r = 3 * divergence
let g = 5 * divergence
let b = 7 * divergence
in (r<<16 | g<<8 | b)
```

Finally we tie it all together - the `main`

function computes the
point of divergence for each pixel, then colours them:

```
let main (screenX: i32) (screenY: i32) (depth: i32)
(xmin: f32) (ymin: f32) (xmax: f32) (ymax: f32): [screenY][screenX]i32 =
let escapes = mandelbrot screenX screenY depth xmin ymin xmax ymax
in map (\(row: []i32): [screenX]i32 ->
map (escapeToColour depth) row)
escapes
```

We can test our code by compiling it to a standalone program:

```
$ futhark pyopencl mandelbrot.fut
$ echo 3 2 255 -2.23 -1.15 0.83 1.15 | ./mandelbrot
[[0i32, 395790i32, 593685i32], [0i32, 0i32, 0i32]]
```

Of course, it is not very satisfying to look at fractals as arrays of
numerically encoded pixel values. Hence, we pass `--library`

to
`futhark pyopencl`

:

`$ futhark pyopencl --library mandelbrot.fut`

This produces a file `mandelbrot.py`

defining a single Python class
`mandelbrot`

, which we can access from ordinary Python code, as
shown below.

## mandelbrot-visualise.py

We will need to import a PNG encoder, Numpy, and of course the module
produced by `futhark pyopencl`

:

```
import png
import numpy
from mandelbrot import mandelbrot
```

Then we create an instance of the class `mandelbrot`

:

`m = mandelbrot()`

The constructor may take additional arguments specifying which OpenCL
platform and device to use, as well as other configuration parameters.
The class defines a single method, `main`

, corresponding to the main
function of the Futhark program. We define a handful of constants
which we pass to the method:

```
filename='mandelbrot.png'
width=800
height=600
limit=255
minx=-2.23
miny=-1.15
maxx=0.83
maxy=1.15
# The .get() is to obtain a Numpy array instead of a PyOpenCL array.
fut_image=m.main(width, height, limit, minx, miny, maxx, maxy).get()
```

The result value is stored in the variable `fut_image`

. Since we
declared the return type of `main`

to be
`[screenY][screenX]i32`

, the returned value will be a
two-dimensional Numpy array of shape `(screenY,screenX)`

. We cannot
pass this directly to the `png`

library, as it expects a
three-dimensional array explicitly encoding the different colour
channels. Fortunately, this array transformation is easy to do with
Numpy:

```
image=numpy.empty((height,width,3))
image[:,:,0] = (fut_image & 0xFF0000) >> 16
image[:,:,1] = (fut_image & 0xFF00) >> 8
image[:,:,2] = (fut_image & 0xFF)
```

And now we can simply invoke the `png`

library:

```
w = png.Writer(width, height, greyscale=False, alpha=False, bitdepth=8)
with open(filename, 'wb') as f:
w.write(f, numpy.reshape(image, (height, width*3)))
```

The result is this moderately attractive fractal in the file
`mandelbrot.png`

:

A slightly more elaborate Python program, which supports command-line parameters and reports timing, can be found here.

## Entry Points

Every *entry point* in the Futhark program becomes a method in the
generated class. An entry point is any function named `text`

, as
well as any function defined using the keyword `entry`

instead of
`let`

. In most cases, the type of the Futhark function maps easily
to the Python world. For example, a Futhark function accepting three
parameters of types `[][]f64`

, `[]i32`

and `bool`

will be
translated into a Python method accepting a two-dimensional Numpy
array of `numpy.double`

s, a one-dimensional array of `numpy.i32`

s,
and a single `numpy.bool`

. And if the Futhark function returns
`([]i32, f64)`

, the Python method will return a tuple of two values:
a Numpy array of integers and a Numpy double-precision float.

Things are more complicated when the entry point accepts or returns
types that do not correspond easily to Numpy types. Actually, the
reason is that the generated code makes use of Futhark’s internal
value representation, but I’m happy to blame Numpy instead. For
example, a function that accepts an array of pairs
(e.g. `[](i32,f32)`

) will be turned into a method that accepts two
arrays: one of integers and one of floats. Similarly, all tuples are
flattened. This not only means that a Futhark function returning
`(i32, (f32, f32))`

will be turned into a Python method returning a
tuple with three elements. It also means that a Futhark function
taking an argument of type `(f32,f32)`

will be turned into a Python
method accepting *two* arguments, each being a float.

The best workaround is to only use simple types in entry point
functions: return only flat tuples, and accept neither tuples nor
arrays of tuples. You can still use tuples and arrays of tuples in
your function bodies and internal functions, it is only the entry
points that are problematic. The `zip`

and `unzip`

operations are
entirely free in Futhark, so `zip`

ing two passed-in arrays into a
single array of pairs carries no overhead.

## More Examples

We have an implementation of Game of Life that uses Pygame to render the ongoing simulation. It supports several variants of the game rules, some of which look rather interesting when visualised. It is also an example of a program that uses multiple entry points.

We also have an interactive Mandelbrot explorer, the Futhark core of which is very similar to the one described above, but where we have written a Pygame interface that allows interactive scrolling, zooming, etc.

## Write your own!

We are quite interested in developing more interesting use cases for Python-Futhark interop. The best use cases are those that perform a good bit of work on the GPU, to amortise the relatively inefficient host-level Python (not to mention copying back and forth between system memory and the GPU). If you can think of something, or even want to try your hand at implementing it, please contribute!