Optimalizace až na Cost
Naučíme se optimalizovat funkce. Začneme od čisté implementace v Pythonu, zkusíme vyřešit problém v NumPy. Poté si ukážeme, jak postupně přejít k jazykům nižší úrovně, jako je Fortran nebo C, a také si představíme balík Numba.
%pylab inline --no-import-all
from IPython.display import Image
Základní koncepce optimalizace¶
Už víme, že Python je jazyk, ve kterém jde velice efektivně programovat. A díky balíkům, jako jsou NumPy, SciPy nebo SymPy jde velice rychle řešit různé vědecké úlohy. Je pochopitelné, že s takovouto mírou abstrakce může být výsledný program pomalejší, než kdyby byl dobře napsán v nějakém kompilovaném jazyce typu C/C++ nebo Fortran. Musíme si ovšem uvědomit, že efektivita programu se měří v celkovém čase stráveném na vývoji a běhu programu (pro daný soubor úloh). Schématicky můžeme znázornit závislost rychlosti běhu programu v závislosti na délce vývoje asi takto:
Image(filename='optimizing-what.png')
Všimněte si například, že typicky existuje nějaký offset ve vývojovém čase, tj. trvá nám déle v nízkoúrovňovém jazyce, než vůbec dostaneme první výsledek. Potřeba optimalizace tedy silně závisí na objemu výpočtů, které budeme s daným programem řešit.
Toto ovšem neznamená, že pokud je objem velký, máme hned začít programovat v C nebo Fortranu. Za chvíli si ukážeme, jak optimalizaci řešit chytřeji a postupně. Empirické pravidlo říká, že 90 % výpočetního času zabere 10 % zdrojového kódu. Jedná se o konkrétní příklad obecného Paretova 80 / 20 principu. Je tedy vhodné nejprve těchto 10 % najít a poté je teprve začít optimalizovat. Python nám k tomuto účelu poskytuje velice mocné nástroje.
Profilování¶
Profilování je nástroj, který nám umožní najít kritická místa v našem programu, oněch 10 %, které stojí za to optimalizovat. Zkusme si to ukázat na jednoduchém příkladu.
def heavy_calc(X):
Y = X.copy()
for i in range(10):
Y = Y**i
return Y
def heavy_loop(inputs):
res = []
for X in inputs:
res.append(heavy_calc(X))
return res
def code_setup():
from numpy.random import rand
N = 20
M = 1000
print("Will generate {} random arrays".format(N))
inputs = [rand(M, M) for n in range(N)]
print("Will calculate now")
result = heavy_loop(inputs)
print("Finished calculation")
Python obsahuje dva základní mofuly pro profilování - profile a cProfile, z nichž ten druhý je rychlejší. Pomocí funkce run pustíme výpočet pod dohledem cProfile, výsledky uložíme do souboru.
import cProfile
cProfile.run('code_setup()', 'pstats')
Dále budeme potřebovat modul pstats, který nám umožní s výsledky pracovat. Použije k tomu třídu Stats.
from pstats import Stats
p = Stats('pstats')
print_stats nám zobrazí prvních n záznamů.
p.print_stats(10)
Ty jsou ovšem nesetříděné. Následující výstup už je užitečnější, záznamy jsou totiž setříděné podle celkového času stráveného v dané funkci. Navíc strip_dirs odstraní adresáře ze jmen funkcí.
p.strip_dirs().sort_stats('cumulative').print_stats(10)
Takto vypadá výstup setříděný pomocí nekumulovaného času.
p.sort_stats('time').print_stats(10)
Jupyter nám může usnadnit práci pomocí %prun a %%prun. Např.
%prun -s cumulative -l 10 code_setup()
Vzorová úloha - vzdálenost množiny bodů ve vícerozměrném prostoru¶
(Tento příklad byl převzat z http://jakevdp.github.io/blog/2013/06/15/numba-vs-cython-take-2.)
Zadání je jednoduché: pro M bodů v N rozměrném prostoru spočítejte vzájemnou vzdálenost $d$, která je pro dva body $x,y$ definovaná jako $\sqrt {\sum_{i=1}^N {{{\left( {{x_i} - {y_i}} \right)}^2}} } $. Výslekem je tedy (symetrická) matice $M\times M$.
# toto nechť jsou naše vstupní data
M = 1000
N = 3
X = np.random.random((M, N))
Implementace v čistém Pythonu¶
Nemůžeme asi očekávat, že toto bude nejrychlejší a nejsnadnější verze našeho programu. Přesto stojí za to ji vyzkoušet, navíc ji budeme ještě potřebovat.
def pairwise_python(X):
M = X.shape[0]
N = X.shape[1]
D = np.empty((M, M), dtype=np.float)
for i in range(M):
for j in range(M):
d = 0.0
for k in range(N):
tmp = X[i, k] - X[j, k]
d += tmp * tmp
D[i, j] = np.sqrt(d)
return D
Tahle funkce nám bude pomáhat ukládat výsledné časy z %timeit.
Do pairwise_times si uložíme výsledné časy.
import collections
pairwise_times = collections.OrderedDict()
timings = %timeit -o pairwise_python(X)
pairwise_times['plain_python'] = timings
To samé pomocí NumPy¶
V případě NumPy můžeme v tomto případě využít broadcasting. Celá funkce tak zabere jeden rádek.
def pairwise_numpy(X):
return np.sqrt(((X[:, np.newaxis, :] - X) ** 2).sum(-1))
Zkusíme, jestli výsledky jsou stejné pomocí assert a numpy.all.
assert np.all(pairwise_numpy(X) == pairwise_python(X))
Nejsou - máme to špatně! Nebo možná jen trochu nepřesně? To ověříme pomocí numpy.allclose.
assert np.allclose(pairwise_numpy(X), pairwise_python(X), rtol=1e-10, atol=1e-15)
Výsledky jsou stejné až na velmi malé rozdíly - to je nebezpečí numerických výpočtů s konečnou přesností.
timings = %timeit -o pairwise_numpy(X)
pairwise_times['numpy'] = timings
Vidíme, že jsme zkrátili běh programu více než 100-krát. To není špatné, navíc je implementace daleko jednodušší.
Přichází Cython¶
Cython je nástroj, který z Python programu, obohaceného o nějaké Cython direktivy, vytvoří program v C (případně C++), který je možné zkompilovat a okamžitě použít jako modul v Pythonu. Typickým příkladem Cython direktiv jsou statické typy. Cython samozřejmě umožňuje používat funkce z binárních knihoven s C rozhraním.
Zkusíme optimalizovat naší funkci pairwise_python.
- Cython zdroják má koncovku .pyx (za začátku byl Pyrex).
- Cython dokáže přeložit jakýkoli Python. Výsledkem je ale minimální (nebo spíš žádná) optimalizace.
-
cimportje analogieimport, pracuje ale s Cython definicemi funkcí (.pxd soubory). - Cython dodává
numpy.pyx, obsahující dodatečné informace pro kompilace NumPy modulů. Proto volámecimport numpy. -
Podobně
libcje speciální modul Cythonu. -
Funkce se deklarují (moho deklarovat) se statickými typy vstupních parametrů. My použijeme
np.ndarray[np.float64_t, ndim=2]. - Proměnné se deklarují pomocí
cdef.
%%file cyfuncs.pyx
import numpy as np
# numpy pro Cython
cimport numpy as np
from libc.math cimport sqrt
# tohle je čistý Python
def pairwise0(X):
M = X.shape[0]
N = X.shape[1]
D = np.empty((M, M), dtype=np.float)
for i in range(M):
for j in range(M):
d = 0.0
for k in range(N):
tmp = X[i, k] - X[j, k]
d += tmp * tmp
D[i, j] = np.sqrt(d)
return D
# tady už začínáme optimalizovat, změny ale nejsou drastické
def pairwise1(np.ndarray[np.float64_t, ndim=2] X):
cdef int M = X.shape[0]
cdef int N = X.shape[1]
cdef double tmp, d
cdef np.ndarray D = np.empty((M, M), dtype=np.float64)
for i in range(M):
for j in range(M):
d = 0.0
for k in range(N):
tmp = X[i, k] - X[j, k]
d += tmp * tmp
D[i, j] = sqrt(d)
return D
%%file setup.py
from distutils.core import setup
from Cython.Build import cythonize
import numpy
setup(
name='cyfuncs',
include_dirs=[numpy.get_include()],
ext_modules=cythonize("cyfuncs.pyx"),
)
!python setup.py build_ext --inplace
Jak jsme již říkali, Cython vytvoří C zdroják, který se pak kompiluje pomocí běžného překladače (např. gcc). Pojďme se na tento soubor podívat.
from IPython.display import FileLink
FileLink('cyfuncs.c')
Ten soubor je dlouhý ... Obsahuje spoustu Python "balastu", na kterém vidíme, jak je vlastně samotný CPython naprogramován. Naštěstí tento soubor obsahuje i komentáře, které říkají, které řádce daný blok odpovídá. Např.
/* "cyfuncs.pyx":16
* tmp = X[i, k] - X[j, k]
* d += tmp * tmp
* D[i, j] = np.sqrt(d) # <<<<<<<<<<<<<<
* return D
*
*/import cyfuncs
print("cyfuncs obsahuje: " + ", ".join((d for d in dir(cyfuncs) if not d.startswith("_"))))
Podívejme se, jestli dostávám stále stejné výsledky.
assert np.allclose(pairwise_numpy(X), cyfuncs.pairwise1(X), rtol=1e-10, atol=1e-15)
No a jak jsme na tom s časem?
timings = %timeit -o cyfuncs.pairwise0(X)
pairwise_times['cython0'] = timings
timings = %timeit -o cyfuncs.pairwise1(X)
pairwise_times['cython1'] = timings
IPython %%cython magic¶
IPython, tak jako v mnoha jiných případech, nám práci s Cythonem usnadňuje pomocí triku %%cython. Zkusíme ho použít. Zároveň zkusíme ještě více náš kód optimalizovat, zatím je totiž pomalejší než numpy.
%load_ext Cython
%%cython -a
import numpy as np
cimport numpy as np
cimport cython
from libc.math cimport sqrt
@cython.boundscheck(False)
@cython.wraparound(False)
def pairwise_cython(double[:, ::1] X):
cdef int M = X.shape[0]
cdef int N = X.shape[1]
cdef double tmp, d
cdef double[:, ::1] D = np.empty((M, M), dtype=np.float64)
for i in range(M):
for j in range(M):
d = 0.0
for k in range(N):
tmp = X[i, k] - X[j, k]
d += tmp * tmp
D[i, j] = sqrt(d)
return np.asarray(D)
assert np.allclose(pairwise_numpy(X), pairwise_cython(X), rtol=1e-10, atol=1e-15)
timings = %timeit -o pairwise_cython(X)
pairwise_times['cython2'] = timings
Tohle je už konečně výrazné zrychlení oproti NumPy.
Cython toho nabízí mnoho¶
Podívejte se na http://docs.cython.org co všechno Cython nabízí -- není toho málo, např.
- použití C++
- šablony (templates)
- OpenMP (k tomu se možná ještě dostaneme)
- vytváření C-API
- třídy (cdef classes)
Z Fortranu (nebo C) do Pythonu pomocí F2PY¶
F2PY je nástroj, který byl v podstatě vytvořen pro NumPy a SciPy, protože, jak dobře víme, tyto moduly volají externí knihovny napsané ve Fortrane nebo C. Dokumentaci (trochu zastaralou) najdeme zde. Bylo tedy velice výhodné vytvořit nástroj, který toto usnadní. A tak se zrodilo F2PY. Ve zkratce, F2PY umožňuje velice jednoduše z Fortran nebo C funkcí vytvořit Python modul. Využívá navíc vlastností NumPy pro předávání vícerozměnrých polí.
Poďme chvilku programovat ve Fortranu :)
%%file pairwise_fort.f90
subroutine pairwise_fort(X, D, m, n)
integer :: n,m
double precision, intent(in) :: X(m, n)
double precision, intent(out) :: D(m, m)
integer :: i, j, k
double precision :: r
do j = 1,m
do i = 1,m
r = 0
do k = 1,n
r = r + (X(i,k) - X(j,k)) * (X(i,k) - X(j,k))
end do
D(i,j) = sqrt(r)
end do
end do
end subroutine pairwise_fort
Z čeho f2py bere informace o vytvoření modulu?
- Pole (double precision) převádí na numpy array.
-
intent(in)= vstupní argument. -
intent(out)= výstupní argument. - f2py schová explicitně zadané rozměry polí (m, n).
Pokud bychom programovali v C, je potřeba dodat f2py nějaké informace navíc, neboť např. intent v C neexistuje.
Tento soubor přeložíme pomocí f2py:
!f2py -c pairwise_fort.f90 -m pairwise_fort
-m pairwise_fort je požadované jméno modulu. Můžeme ho rovnou importovat, resp. jeho stejnojmennou funkci.
from pairwise_fort import pairwise_fort
print(pairwise_fort.__doc__)
Fortran a C používají jiné uspořádání paměti pro ukládání vícerozměrných polí. Fortran je "column-major" zatímco C je "row-major". NumPy dokáže pracovat s obojím a pro uživatele je to obvykle jedno. Pokud ovšem chceme předat vícerozměrné pole do Fortran funkce, je lepší mít prvky uložené v paměti jako to dělá Fortran. V takovém případě totiž f2py předá pouze referenci (ukazatel) na dané místo v paměti. V opačném případě f2py nejprve pole musí transponovat, tj. vytvořit kopii s jiným uspořádáním, což může být samozřejmě náročné na pamět a procesor.
Vytvoříme si proměnnou XF, která má Fortran uspořádání, pomocí numpy.asfortranarray (prozaický název :)
XF = np.asfortranarray(X)
Vyzkoušíme, jestli stále dostáváme stejné výsledky.
assert np.allclose(pairwise_numpy(X), pairwise_fort(X), rtol=1e-10, atol=1e-15)
No a konečně se můžeme podívat, jak je to s rychlostí ...
timings = %timeit -o pairwise_fort(X)
pairwise_times['fortran'] = timings
Představuje se numba¶
Numba kompiluje Python kód pomocí LLVM. Podporujme just-in-time kompilaci pomocí dekorátoru jit (http://numba.pydata.org/numba-doc/latest/user/jit.html).
@numba.jit(
signature=None,
nopython=False,
nogil=False,
cache=False,
forceobj=False,
parallel=False,
error_model='python',
fastmath=False, locals={}
)import numba
pairwise_numba = numba.jit(pairwise_python)
Tradiční kontrola. Po prvním spuštění navíc Numba funkci poprvé zkompiluje.
assert np.allclose(pairwise_numpy(X), pairwise_numba(X), rtol=1e-10, atol=1e-15)
Jaký čas od takto "jednoduché" optimalizace můžeme očekávat?
timings = %timeit -o pairwise_numba(X)
pairwise_times['numba'] = timings
Vidíme, že zrychlení je výborné - jsme na úrovni zatím nejlepšího výsledku! A navíc že jsme toho dosáhli jediným řádkem (kromě importů).
Ještě můžeme zkusit výledek vylepšit pomocí paralelizace, nopython a / nebo fastmath režimu. Pro nopython musíme vytvořit výsledný numpy objekt vně kompilované funkce. Paralelizace docílíme pomocí parallel=True a numba.prange. Všimněte si použití @jit jako dekorátoru.
@numba.jit(nopython=True, parallel=True, fastmath=True)
def _pairwise_nopython(X: numpy.ndarray, D: numpy.ndarray) -> numpy.ndarray:
M = X.shape[0]
N = X.shape[1]
for i in numba.prange(M):
for j in numba.prange(M):
d = 0.0
for k in range(N):
tmp = X[i, k] - X[j, k]
d += tmp * tmp
D[i, j] = numpy.sqrt(d)
return D
def pairwise_numba_fast_parallel(X: numpy.ndarray) -> numpy.ndarray:
D = numpy.empty((X.shape[0], X.shape[0]), dtype = numpy.float)
_pairwise_nopython(X, D)
return D
assert np.allclose(pairwise_numpy(X), pairwise_numba_fast_parallel(X), rtol=1e-10, atol=1e-15)
timings = %timeit -o pairwise_numba_fast_parallel(X)
pairwise_times['numba_fast_parallel'] = timings
Srovnání výsledků¶
Výsledky můžeme porovnat pomocí grafu.
pairwise_times
fig, ax = plt.subplots()
values = np.array([t.average for t in pairwise_times.values()])
x = range(len(pairwise_times))
ax.bar(x, values)
ax.set_xticks(x)
ax.set_xticklabels(tuple(pairwise_times.keys()), rotation='vertical')
ax.set_ylabel('time [ms]')
ax.set_yscale('log')
Další možnosti¶
- Vestavěný modul
ctypesdovoluje volat funkce z externích dynamických knihoven. Pro použití s NumPy viz Cookbook. - Alternativou k ctypes je cffi.
- CuPy využívá GPU.
- numexpr dokáže kompilovat Numpy výrazy.
- Theano se zaměřuje na strojové učení, také optimalizuje vektorové (maticové) operace, dovoluje je spouštět na CPU.
- Nuitka kompiluje Python, ale na rozdíl od Numba nespecializuje funkce na základě typů.
- SWIG jde použít pro propojení s mnoha jazyky.'
Cvičení¶
Obdobný postup aplikujte na výpočet kumulativního součtu, který je definovaný jako
$\displaystyle S_j = \sum\limits_{i = 1}^j x_i $
Výsledky a časování porovnejte s numpy.cumsum.
Nápověda: Ve vaší funkci vytvořte nejprve výsledné numpy pole pomocí numpy.empty_like.
Komentáře
Comments powered by Disqus