from IPython.display import Image
Image("https://raw.github.com/jrjohansson/scientific-python-lectures/master/images/optimizing-what.png")
Katkada korištenje Numpy-a ne daje dovoljno ubrzanje, ili je teško/nespretno vektorizirati kod. Tada postoji više opcija
%prun
) možemo napisati npr. u C-u. U Pythonu je taj proces prilično bezbolan.# olakšava rad u Cythonu u IPythonu
%load_ext Cython
a=10
b=20
%%cython_inline
magična funkcija koristi Cython.inline
da bi se kompajlirao Cython izraz. return
služi za slanje izlaza.
%%cython_inline
return a+b
Compiling /projects/95cb524c-f1ac-4a19-908e-5f47c10f47c1/.cython/inline/_cython_inline_6941fceee5936d163c6c3563c2edbe1e.pyx because it changed. [1/1] Cythonizing /projects/95cb524c-f1ac-4a19-908e-5f47c10f47c1/.cython/inline/_cython_inline_6941fceee5936d163c6c3563c2edbe1e.pyx
warning: /projects/95cb524c-f1ac-4a19-908e-5f47c10f47c1/.cython/inline/_cython_inline_6941fceee5936d163c6c3563c2edbe1e.pyx:6:4: Unreachable code
30
%%cython_pyximport
magična funkcija služi da se unese proizvoljan Cython kod u IPython notebook ćeliju. Taj kod se sprema u .pyx
u radnom direktoriju i onda importira koristeći pyximport
. Moramo specificirati ime modula (u donjem slučaju foo
). Svi objekti modula se automatski importiraju.
%%cython_pyximport foo
def f(x):
return 4.0*x
f(10)
40.0
U dosadašnjim primjerima Cython kod se nije razlikovao od Python koda. U sljedećem primjeru ćemo vidjeti neka od proširenja Pythona koja nudi Cython.
Magična funkcija %%cython
je slična funkciji %%cython_pyximport
, ali ne traži ime modula te sve datoteke sprema u privremeni direktorij.
Jedan od načina računanja broja $\pi$ je korištenjem formule
$$ \frac{\pi}{4} = \int_0^1 \sqrt{1-x^2}\,\mathrm{d}x.$$
A integral možemo aproksimirati pomoću trapezne formule, što nam daje
$$\frac{\pi}{4}\approx \frac{1}{n}\left(\frac{1}{2}+\sum_{i=1}^n\sqrt{1-\left(\frac{i}{n}\right)^2} \right) $$
Krenimo od običnog Pythona.
from math import sqrt
def funkcija(x):
return sqrt(1-x**2)
def integral4pi(n):
korak = 1.0/n
rez = (funkcija(0)+funkcija(1))/2
for i in range(n):
rez += funkcija(i*korak)
return 4*rez*korak
approx=integral4pi(10**7)
print ('pi={}'.format(approx))
pi=3.1415930535527115
%timeit integral4pi(10**7)
1 loop, best of 3: 4.17 s per loop
Ovdje je
cimport
: ekvivalent importa, ali možemo učitavati i iz C ili C++ bibliotekacdef
: ekvivalent za def, gdje definiramo (kao u C-u) tipove podataka; također služi za deklariranje tipa varijabli%%cython
cimport cython
from libc.math cimport sqrt
cdef double cy_funkcija(double x):
return sqrt(1-x**2)
def cy_integral4pi(int n):
cdef double korak, rez
cdef int i
korak = 1.0/n
rez = (cy_funkcija(0)+cy_funkcija(1))/2
for i in range(n):
rez += cy_funkcija(i*korak)
return 4*rez*korak
cy_approx = cy_integral4pi(10**7)
print ('pi={}'.format(cy_approx))
pi=3.1415930535527115
%timeit cy_integral4pi(10**7)
10 loops, best of 3: 72 ms per loop
Ovdje je
@
: Pythonova sintaksa za tzv. dekoratorenogil
: GIL je skraćenica od global interpreter lock, koji spriječava simultano izvršavanje koda u više niti (eng. threads), ovdje s nogil
deklariramo da je sigurno pozvati funkciju bez GIL-awith nogil
: odgovarajući dio koda ne koristi GIL%%cython
cimport cython
from libc.math cimport exp, sqrt, pow, log, erf
@cython.cdivision(True)
cdef double std_norm_cdf(double x) nogil:
return 0.5*(1+erf(x/sqrt(2.0)))
@cython.cdivision(True)
def black_scholes(double s, double k, double t, double v, double rf, double div, double cp):
"""s : početna cijena dionice
k : fiksirana cijena (strike price)
t : vrijeme
v : volatilnost
rf : bezrizična kamata
div : dividenda
cp : put/call paritet
"""
cdef double d1, d2, optprice
with nogil:
d1 = (log(s/k)+(rf-div+0.5*pow(v,2))*t)/(v*sqrt(t))
d2 = d1 - v*sqrt(t)
optprice = cp*s*exp(-div*t)*std_norm_cdf(cp*d1) - cp*k*exp(-rf*t)*std_norm_cdf(cp*d2)
return optprice
black_scholes(100.0, 100.0, 1.0, 0.3, 0.03, 0.0, -1)
10.327861752731728
%timeit black_scholes(100.0, 100.0, 1.0, 0.3, 0.03, 0.0, -1)
The slowest run took 53.24 times longer than the fastest. This could mean that an intermediate result is being cached. 1000000 loops, best of 3: 328 ns per loop
Cython omogućava linkanje dodatnih biblioteka pri kompajliranju, u ovom slučaju standardne matematičke biblioteke
%%cython -lm
from libc.math cimport sin
print ('sin(1)=', sin(1))
sin(1)= 0.8414709848078965
Računamo kumulativnu sumu niza
import numpy as np
def py_dcumsum(a):
b = np.empty_like(a)
b[0] = a[0]
for n in range(1,len(a)):
b[n] = b[n-1]+a[n]
return b
a = np.random.rand(100000)
b = np.empty_like(a)
%%cython
cimport numpy
def cy_dcumsum2(numpy.ndarray[numpy.float64_t, ndim=1] a, numpy.ndarray[numpy.float64_t, ndim=1] b):
cdef int i, n = len(a)
b[0] = a[0]
for i from 1 <= i < n:
b[i] = b[i-1] + a[i]
return b
%timeit cy_dcumsum2(a,b)
1000 loops, best of 3: 328 µs per loop
%timeit py_dcumsum(a)
10 loops, best of 3: 32.4 ms per loop
Numba je just-in-time kompajler koji korisiti LLVM. Korištenje numbe je zapravo dosta jednostavno.
# iz numbe učitavamo jit kompajler
from numba import jit
from numpy import arange
Funkcija sum2d(arr)
će biti optimizirana, a cijela procedura se svela na korištenje jednog dekoratora.
# @jit je dekorator
@jit
def sum2d(arr):
M, N = arr.shape
result = 0.0
for i in range(M):
for j in range(N):
result += arr[i,j]
return result
def py_sum2d(arr):
M, N = arr.shape
result = 0.0
for i in range(M):
for j in range(N):
result += arr[i,j]
return result
a = arange(9).reshape(3,3)
%timeit sum2d(a)
The slowest run took 214864.54 times longer than the fastest. This could mean that an intermediate result is being cached. 1000000 loops, best of 3: 344 ns per loop
%timeit py_sum2d(a)
The slowest run took 10.92 times longer than the fastest. This could mean that an intermediate result is being cached. 100000 loops, best of 3: 5.66 µs per loop
import numpy
def filter2d(image, filt):
M, N = image.shape
Mf, Nf = filt.shape
Mf2 = Mf // 2
Nf2 = Nf // 2
result = numpy.zeros_like(image)
for i in range(Mf2, M - Mf2):
for j in range(Nf2, N - Nf2):
num = 0.0
for ii in range(Mf):
for jj in range(Nf):
num += (filt[Mf-1-ii, Nf-1-jj] * image[i-Mf2+ii, j-Nf2+jj])
result[i, j] = num
return result
from numba import double
fastfilter_2d = jit(double[:,:](double[:,:], double[:,:]))(filter2d)
image = numpy.random.random((100, 100))
filt = numpy.random.random((10, 10))
%timeit fastfilter_2d(image, filt)
1000 loops, best of 3: 1.7 ms per loop
%timeit filter2d(image, filt)
1 loop, best of 3: 499 ms per loop
Drugi način ubrzavanja izvršavanja programa je pomoću paralelizacije. To je opsežna tema sama po sebi te u nju nećemo ulaziti. Popularni način paralelizacije je npr. pomoću Apache Sparka ili pomoću Dask biblioteke.
from verzije import *
from IPython.display import HTML
HTML(print_sysinfo()+info_packages('cython,numpy,numba'))
Python verzija | 3.5.3 |
kompajler | GCC 4.8.2 20140120 (Red Hat 4.8.2-15) |
sustav | Linux |
broj CPU-a | 8 |
interpreter | 64bit |
cython verzija | 0.25.2 |
numpy verzija | 1.11.3 |
numba verzija | 0.31.0 |