Wir backen uns ein schnelles Mandelbrötchen

Die Visualisierung der Mandelbrot-Menge, die im Volksmund oft Apfelmännchen oder nach ihrem Schöpfer, dem französisch-amerikanischen Mathematiker Benoît Mandelbrot, auch Mandelbrötchen genannt wird, ist, obwohl sie durch die einfache Formel z^2 + c erzeugt wird (z und c sind komplexe Zahlen), eine sehr rechenintensive Operation. Zwar geht es mit den heutigen Rechnern schon recht flott, wie man sieht, wenn man dieses IPython-Programm laufen läßt:

import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline

size = 400
iterations = 100

def mandelbrot(m, size, iterations):
    for i in range(size):
        for j in range(size):
            c = -2 + 3./size*j + 1j*(1.5 - 3./size*i)
            z = 0
            for n in range(iterations):
                if np.abs(z) <= 10:
                    z = z*z + c
                    m[i, j] = n
                else:
                    break

m = np.zeros((size, size))
mandelbrot(m, size, iterations)

plt.imshow(np.log(m), cmap = plt.cm.hot)
plt.xticks([])
plt.yticks([])
plt.show()

Als ich in den späten 1980er Jahren diese Menge auf meinem Atari Mega ST 2 (in freundlichem schwarz-weiß) berechnen ließ, dauerte es Stunden, bis sie komplett auf dem Monitor erschien. Da ich nun auch die Zeit messen wollte, die zur Berechnung benötigt wurde, habe ich dieses Programm (und das folgende) ausnahmsweise in einem Jupyter-Notebook laufen lassen, da nur dort die zellbasierte Zeitmessung mit dem magischen Kommande %%timeit durchgeführt werden kann:

%%timeit m = np.zeros((size, size))
mandelbrot(m, size, iterations)

Das Ergebnis war mit

1 loop, best of 3: 8.83 s per loop

schon recht flott, aber es geht noch viel schneller. Dazu wird Numba benötigt, ein Paket, das mit Hilfe der LLVM-Architektur rechenintensive Operationen just in time in Maschinencode überführt, ohne daß der Python-Interpreter verlassen werden muß. Um Numba zu installieren, wird der Anaconda-eigene Paketmanager conda bemüht:

conda install numba

Dann wird mit

import numba
from numba import jit, complex128

das Paket und der Just-in-Time-Compiler jit und der Datentyp complex128 importiert. Zwar versucht in der Regel jit selber, die Datentypen zu erkennen, aber manchmal (zum Beispiel bei komplexen Zahlen) kann es sicherer sein, ihm den Datentyp mitzugeben.

Die Funktion, die die Mandelbrot-Menge berechnet, ist identisch (nur das ich sie zur Unterscheidung mandelbrot_numba genannt habe) aber ich habe sie mit einem Decorator geschmückt, der Numba anweist, diese Funktion on the fly in Maschinencode zu übersetzen:

@jit(locals = dict(c = complex128, z = complex128))
def mandelbrot_numba(m, size, iterations):
    for i in range(size):
        for j in range(size):
            c = -2 + 3./size*j + 1j*(1.5 - 3./size*i)
            z = 0
            for n in range(iterations):
                if np.abs(z) <= 10:
                    z = z*z + c
                    m[i, j] = n
                else:
                    break

m = np.zeros((size, size))

mandelbrot_numba(m, size, iterations)
plt.imshow(np.log(m), cmap = plt.cm.hot)

Ich habe dieses Mal auch die Ticks zeichnen lassen, um bei meinen Experimenten (es lief nicht alles von Anfang an so reibungslos, wie es hier aussieht) den schnellen Plot von dem langsameren Plot unterscheiden zu können.

Nun also auch hier die Zeitmessung mit %%timeit:

%%timeit m = np.zeros((size, size))
mandelbrot_numba(m, size, iterations)

Das Ergebnis war:

10 loops, best of 3: 20.3 ms per loop

Wow, 20,3 ms gegenüber 8,3 s! Das ist doch eine Beschleunigung, die sich sehen lassen kann. Vor allem, wenn man solche Berechnung für die Details der Mandelbrot-Menge mit größerer Iterationstiefe und häufig hintereinander ausführen möchte, lohnen sich die Installation und der Import von Numba auf jeden Fall.

---

(Kommentieren)