Python erlaubt das schnelle Entwickeln und Implementieren von Prototypen und ist dabei sehr geeignet um einen Algorithmus auf seine Eignung zu testen. Für wissenschaftliche Berechnungen steht man häufig vor zwei Problemen: Python kann für Berechnungen sehr langsam sein. Zum anderen soll vorhandener Code und externe Bibliotheken verwendet werden. Dieser Artikel gibt einen Überblick über verschiedene Techniken.
Im Folgenden wird anhand eines einfachen Beispiels auf das Cython, das Scipy module Weave, Ctypes und die Python C-API eingegangen. Der Artikel soll einen Eindruck der jeweiligen Technik vermitteln, kann und soll aber nicht die jeweilige Dokumentation ersetzen.
Im folgenden Beispiel soll der rechenintensivste Teil eines Monte Carlo Algorithmus zur Erzeugung von Photonen optimiert werden. Bei diesem Monte Carlo Experiment wird in jedem Zeitschritt zwischen Fluoreszenz (=0) und FRET (=1) bzw. keines der Beiden gewürfelt werden. Die Wahrscheinlichkeit für Fluoreszenz ist Zeitunabhängig während die FRET Wahrscheinlichkeit Zeitabhängig ist. So sieht die Implementierung in “Pure Python” aus (python.py):
import random
def getPhoton(p0,pvar,start,seed):
"""
p0: probability for donor fluorescence
pvar: probability trajectory for FRET + probability for donor fluorescence
start: start position on trajectory
seed: random number generator seed
"""
random.seed(seed)
for curndx in range(start,len(pvar)):
rnd=random.random()
if rnd < p0:
return 0
if rnd < pvar[curndx]:
return 1
return -1
Die an die Funktion übergebenen Variablen sind die Fluoreszenzwahrscheinlichkeit p0, die FRET Wahrscheinlichkeit pvar (Liste, Tupel oder Numpy Array), der Startzeitpunkt und die Seed für die Initialisierung des Zufallzahlengenerators.
Wir möchten nun die Performance des Codes verbessern. Cython konvertiert Python Code in C, welcher dann kompiliert und als Python Modul verwendet werden kann. Hier erstmal der unmodifizierte Python code (cython.pyx):
import random
def getPhoton(p0,pvar,start,seed):
random.seed(seed)
for curndx in range(start,len(pvar)):
rnd=random.random()
if rnd < p0:
return 0
if rnd < pvar[curndx]:
return 1
return -1
Übersetzt wird der Code in C durch den Aufruf von
cython cython.pyx
um dann anschließend das C-Modul zu übersetzten.
gcc -shared -O3 -fPIC -I/usr/include/python2.7 -o cython.so cython.c
Mit Cython übersetzter unmodifizierter Code läuft häufig nicht viel schneller als das Original. Ein Grund ist die
dynamische Typisierung von Python. Cython erlaubt daher statische Typdeklaration als Erweiterung des Python Syntax.
Somit kann man durch einfaches Hinzufügen von Typendeklarationen oft einen deutlichen Geschwindigkeitszuwachs
bekommen (ctypes-statictypedec.pyx):
cimport numpy as np
from numpy import random
def getPhoton(double p0,np.ndarray pvar,long start,long seed):
cdef double rnd
cdef long end,curndx
end = len(pvar)
random.seed(seed)
for curndx in range(start,end):
rnd=random.random()
if rnd < p0:
return 0
if rnd < pvar[curndx]:
return 1
return -1
Standard Datentypen wie int, long etc. können wie in C deklariert werden. Für spezielle Datentypen wie Numpy Arrays hält cython ein angepasstes Numpy modul vor, welches wir mit cimport für Cython importieren.
Cython ist auch relativ flexibel, wenn es um das Einbinden von externem Code geht. Im folgenden Fall möchten wir einen SIMD optimierten Mersenne Twister Zufallszahlengenerator einbinden. Die Signaturen aus der SFMT.h Header kopieren wir und formen sie wie folgt um (cython-sfmt.pyx):
cdef extern from "SFMT.h":
void init_gen_rand(int seed)
double genrand_res53()
cimport numpy as np
def getPhoton(double p0,np.ndarray pvar,long start,long seed):
cdef double rnd
cdef long end,curndx
end = len(pvar)
init_gen_rand(seed)
for curndx in range(start,end):
rnd=genrand_res53()
if rnd < p0:
return 0
if rnd < pvar[curndx]:
return 1
return -1
Der cdef der ersten 3 Zeilen ermöglicht das verwenden der Funktionen im “Python Code”. Wir können z.B. die beiden Objektdateien von SFMT und Cython wie folgt zu einem Python C-Modul linken:
gcc -shared -fPIC -o cython-sfmt.so SFMT.o cython-sfmt.o
Zusammenfassend ist Cython besonders gut geeignet um bestehenden Python Code zu Optimieren und um externen C-Code bzw Bibliotheken einzubinden.
Was, wenn ein Codeschnipsel in C aus Python heraus aufgerufen werden soll? Hier bietet Weave mit dem Typenconverter Blitz eine einfache Möglichkeit (weave.py):
from scipy.weave import inline
from scipy.weave import converters
def getPhoton(pconst,u,startpos,rngseed):
code = """
double rnd;
int retval,i;
retval=-1;
srand(rngseed);
for (i=startpos;i<Nu[0];i++)
{
rnd = rand()/((float)RAND_MAX);
if (rnd < pconst) {
retval=0;
break;
}
if (rnd < u(i)){
retval=1;
break;
}
}
return_val = retval;
"""
return inline(code, ['pconst','u','startpos','rngseed'],
type_converters=converters.blitz,compiler='gcc',
headers=['<stdlib.h>',],)
Wie man sieht ist besteht die getPhoton Funktion ausschließlich aus einer Variablenzuweisung eines multiline Strings zu code welcher in der letzten Zeile mittels der inline Funktion “on-the-fly” übersetzt wird. In der Inline funktion werden der Codestring und die Liste der Parameter spezifiziert. Außerdem werden Typconverter, Compiler und Header spezifiziert. Die Übersetzung erfolgt zur Laufzeit und übersetzte Module werden gecached, so daß diese nur einmalig anfällt.
Meiner Einschätzung nach ist Weave v.a. für kurze einfache Codeabschnitte geeignet. Zum einen ist das einbinden externer Objektdateien und Bibliotheken etwas hakelig, zum anderen kann der C-Syntax üblicherweise nicht hervorgehoben werden, da dieser einen Python String darstellt.
Nehmen wir nun an, daß wir den Code für den Photonengenerator in einer externen shared Library verfügbar haben und darauf zugreifen möchten. Hierfür bietet sich ctypes aus der Python Standardbibliothek an. Ctypes ermöglicht das direkte Aufrufen der Funktion und entsprechende Konvertierung von Python Objekten in Objekte die die C Typen abbilden. Was hier (ctypes.py) folgt ist ein sehr einfaches Beispiel. Ctypes ermöglicht auch die Abbildung von Strukturen, das spezifizieren von Funktionsinput und -output Variablen u.v.m..
from ctypes import cdll,c_double,byref,c_long,POINTER
import ctypes
myclib = cdll.LoadLibrary("libfretext.so")
def getPhoton(p0,pvar,start,seed):
nvar = len(pvar)
doublepointer = ctypes.POINTER(ctypes.c_double)
c_pconst = c_double(p0)
c_pvarlen = c_long(nvar)
c_pvarp = pvar.ctypes.data_as(doublepointer)
c_startpos = c_long(start)
c_rngseed = c_long(seed)
retval = myclib.tryGetPhoton(c_pconst,c_pvarlen, c_pvarp,c_startpos, c_rngseed)
return retval
Zuerst “importieren” wir die Bibliothek “libfretext.so” über die LoadLibrary Zeile als myclib. Die Folgende Funktion getPhoton hat als einzige Aufgabe die Konvertierung von Python Objekten. Für das Array erzeugen wir uns zuerst den Typ “doublepointer”, welchen wir später der ctypes.data_as Methode des Numpy Arrays übergeben. Die so konvertierten Python Objekte können anschließend beim Funktionsaufruf übergeben werden.
Ctypes ist vorallem zum Aufruf von Funktionen in bestehenden Bibliotheken geeignet. Es ermöglicht umfangreiche Kapselung von bestehenden Funktionen in Bibliotheken ohne deren Modifikation.
Die umfangreichsten Möglichkeiten bietet die C-API von Python. Diese kann verwendet werden um Module in C zu schreiben oder aber um Python Interpreter in bestehende C-Anwendungen einzubetten.
#include <Python.h>
#include <numpy/arrayobject.h>
#include "cext.h"
#include "SFMT.h"
static PyObject *getPhoton(PyObject *self, PyObject *args)
{
PyArrayObject *pvar;
long startpos,rngseed;
double pconst,*pvarp,rnd;
npy_intp *pvarlen;
int retval,i,j;
if (!PyArg_ParseTuple(args,"dOll",pconst,pvar,startpos,rngseed)) return NULL;
/*Checking if array is C-Contiguous, double and then initializing pointers*/
if (!PyArray_ISCONTIGUOUS(pvar) || (sizeof(double) != PyArray_ITEMSIZE(pvar))) return NULL;
pvarp=(double*)PyArray_DATA(pvar);
pvarlen = PyArray_DIMS(pvar);
retval=-1;
init_gen_rand(rngseed);
for (i=startpos;i<pvarlen[0];i++){
rnd = genrand_res53();
if (rnd < pconst) {
retval=0;
break;
}
if (rnd < pvarp[i]){
retval=1;
break;
}
}
return Py_BuildValue("i",retval);
}
static PyMethodDef fretmethods[] = {
{"getPhoton", getPhoton,METH_VARARGS,"Generation attempt of a photon in plain C-Code"},
{NULL,NULL,0,NULL}
};
#ifndef PyMODINIT_FUNC /* declarations for DLL import/export */
#define PyMODINIT_FUNC void
#endif
PyMODINIT_FUNC initcext(void)
{
(void) Py_InitModule("cext",fretmethods);
import_array();
printf("-> FRET numpy extension initialized using SF Mersenne Twister\n");
}
Der gezeigte Code besteht im Wesentlichen aus 3 Blöcken: Die eigentliche Funktion getPhoton, welche Python Objekte entgegennimmt und zurückgibt. Dem PyMethodDef Block, in dem die im Modul enthaltenen Funktionen aufgelistet werden und der initcext Funktion, die bei der Modulinitialisierung aufgerufen wird. In obigem Beispiel wird zusätzlich die Numpy C-API verwendet.
Übersetzt und gelinkt wird mit
gcc -I/usr/lib64/python2.7/site-packages/numpy/core/include/ -I/usr/include/python2.7 -c -O3 -fPIC -D_FORTIFY_SOURCE=0 -DMEXP=607 -msse2 -DHAVE_SSE2 SFMT.c
gcc -I/usr/lib64/python2.7/site-packages/numpy/core/include/ -I/usr/include/python2.7 -c -O3 -fPIC -D_FORTIFY_SOURCE=0 -DMEXP=607 -msse2 -DHAVE_SSE2 cext.c
gcc -shared -fPIC -Wl,-soname,cext.so -o cext.so SFMT.o cext.o
Wer das Maximum an Kontrolle, Flexibilität und Geschwindikeit erhalten möchte, sollte sich mit der C-API vertraut machen. Sie erfordert jedoch auch eine umfangreiche Beschäftigung mit den Python Interna und deren Implementierung in C, als Beispiel sei hier nur die Inkrementierung und Dekrementierung von Pythonobjektreferenzzählern genannt.