In [1]:
from IPython.display import Image
Image("https://raw.github.com/jrjohansson/scientific-python-lectures/master/images/optimizing-what.png")
Out[1]:
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.
In [2]:
# olakšava rad u Cythonu u IPythonu
%load_ext Cython
In [3]:
a=10
b=20
%%cython_inline magična funkcija koristi Cython.inline da bi se kompajlirao Cython izraz. return služi za slanje izlaza.
In [4]:
%%cython_inline
return a+b
Out[4]:
%%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.
In [5]:
%%cython_pyximport foo
def f(x):
return 4.0*x
In [6]:
f(10)
Out[6]:
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.
In [7]:
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
In [8]:
approx=integral4pi(10**7)
print ('pi={}'.format(approx))
In [9]:
%timeit integral4pi(10**7)
In [10]:
%%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
In [11]:
cy_approx = cy_integral4pi(10**7)
print ('pi={}'.format(cy_approx))
In [12]:
%timeit cy_integral4pi(10**7)
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
In [13]:
%%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
In [14]:
black_scholes(100.0, 100.0, 1.0, 0.3, 0.03, 0.0, -1)
Out[14]:
In [15]:
%timeit black_scholes(100.0, 100.0, 1.0, 0.3, 0.03, 0.0, -1)
Cython omogućava linkanje dodatnih biblioteka pri kompajliranju, u ovom slučaju standardne matematičke biblioteke
In [33]:
%%cython -lm
from libc.math cimport sin
print ('sin(1)=', sin(1))
In [34]:
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
In [35]:
a = np.random.rand(100000)
b = np.empty_like(a)
In [36]:
%%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
In [37]:
%timeit cy_dcumsum2(a,b)
In [38]:
%timeit py_dcumsum(a)
In [39]:
# 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.
In [40]:
# @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
In [41]:
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
In [42]:
a = arange(9).reshape(3,3)
In [43]:
%timeit sum2d(a)
In [44]:
%timeit py_sum2d(a)
In [45]:
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)
In [46]:
image = numpy.random.random((100, 100))
filt = numpy.random.random((10, 10))
In [47]:
%timeit fastfilter_2d(image, filt)
In [48]:
%timeit filter2d(image, filt)
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.
In [49]:
from verzije import *
from IPython.display import HTML
HTML(print_sysinfo()+info_packages('cython,numpy,numba'))
Out[49]: