{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "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.\n", "" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2024-04-18T08:17:08.363327Z", "start_time": "2024-04-18T08:17:08.072665Z" }, "collapsed": false }, "outputs": [], "source": [ "%matplotlib inline\n", "from IPython.display import Image\n", "\n", "import matplotlib.pyplot as plt\n", "import numpy as np" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Základní koncepce optimalizace\n", "\n", "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:" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2024-04-18T08:17:10.193420Z", "start_time": "2024-04-18T08:17:10.188542Z" }, "collapsed": false }, "outputs": [], "source": [ "Image(filename='optimizing-what.png')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "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*.\n", "\n", "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](https://cs.wikipedia.org/wiki/Paret%C5%AFv_princip). 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." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Profilování\n", "\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." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2024-04-18T08:17:12.089139Z", "start_time": "2024-04-18T08:17:12.084944Z" }, "collapsed": false }, "outputs": [], "source": [ "def heavy_calc(X):\n", " Y = X.copy()\n", " for i in range(10):\n", " Y = Y**i\n", " return Y\n", "\n", "def heavy_loop(inputs):\n", " res = []\n", " for X in inputs:\n", " res.append(heavy_calc(X))\n", " return res\n", "\n", "def code_setup():\n", " from numpy.random import rand\n", " N = 20\n", " M = 1000\n", " print(\"Will generate {} random arrays\".format(N))\n", " inputs = [rand(M, M) for n in range(N)]\n", " print(\"Will calculate now\")\n", " result = heavy_loop(inputs)\n", " print(\"Finished calculation\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "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." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2024-04-18T08:17:16.723579Z", "start_time": "2024-04-18T08:17:15.663576Z" }, "collapsed": false }, "outputs": [], "source": [ "import cProfile\n", "cProfile.run('code_setup()', 'pstats')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Dále budeme potřebovat modul `pstats`, který nám umožní s výsledky pracovat. Použije k tomu třídu `Stats`." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2024-04-18T08:17:18.050402Z", "start_time": "2024-04-18T08:17:18.047259Z" }, "collapsed": false }, "outputs": [], "source": [ "from pstats import Stats\n", "p = Stats('pstats')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "`print_stats` nám zobrazí prvních n záznamů." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2024-04-18T08:17:20.961915Z", "start_time": "2024-04-18T08:17:20.956328Z" }, "collapsed": false }, "outputs": [], "source": [ "p.print_stats(10)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "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í." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2024-04-18T08:17:31.098424Z", "start_time": "2024-04-18T08:17:31.093655Z" }, "collapsed": false }, "outputs": [], "source": [ "p.strip_dirs().sort_stats('cumulative').print_stats(10)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Takto vypadá výstup setříděný pomocí nekumulovaného času." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2024-04-11T12:49:44.698688Z", "start_time": "2024-04-11T12:49:44.695290Z" }, "collapsed": false }, "outputs": [], "source": [ "p.sort_stats('time').print_stats(10)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Jupyter nám může usnadnit práci pomocí `%prun` a `%%prun`. Např." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2024-04-18T08:21:57.247796Z", "start_time": "2024-04-18T08:21:55.778256Z" } }, "outputs": [], "source": [ "%prun -s cumulative -l 10 code_setup()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Z obou výstupů celkem jasně vydíme, že naprostou většinu času trávíme ve funkci `heavy_calc`. Pokud se tedy chceme pustit do optimalizace, musíme se zaměřit právě na tuto část našeho programu.\n", "\n", "Výsledky můžete navíc spojit s nástroji pro vizualizaci, např.[SnakeViz](http://jiffyclub.github.io/snakeviz/) nebo [vprof](https://github.com/nvdv/vprof), popř. pokročilý editor jako [PyCharm](https://www.jetbrains.com/pycharm/)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Vzorová úloha - vzdálenost množiny bodů ve vícerozměrném prostoru\n", "\n", "(Tento příklad byl převzat z http://jakevdp.github.io/blog/2013/06/15/numba-vs-cython-take-2.)\n", "\n", "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$." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2024-04-18T08:22:00.402824Z", "start_time": "2024-04-18T08:22:00.399751Z" }, "collapsed": false }, "outputs": [], "source": [ "# toto nechť jsou naše vstupní data\n", "M = 1000\n", "N = 3\n", "X = np.random.random((M, N))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Implementace v čistém Pythonu\n", "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." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2024-04-18T08:22:01.347164Z", "start_time": "2024-04-18T08:22:01.342692Z" }, "collapsed": false }, "outputs": [], "source": [ "def pairwise_python(X):\n", " M = X.shape[0]\n", " N = X.shape[1]\n", " D = np.empty((M, M), dtype=float)\n", " for i in range(M):\n", " for j in range(M):\n", " d = 0.0\n", " for k in range(N):\n", " tmp = X[i, k] - X[j, k]\n", " d += tmp * tmp\n", " D[i, j] = np.sqrt(d)\n", " return D" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Tahle funkce nám bude pomáhat ukládat výsledné časy z `%timeit`." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Do `pairwise_times` si uložíme výsledné časy." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2024-04-18T08:22:02.796241Z", "start_time": "2024-04-18T08:22:02.792442Z" }, "collapsed": false }, "outputs": [], "source": [ "pairwise_times = {}" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2024-04-18T08:22:15.299522Z", "start_time": "2024-04-18T08:22:03.284684Z" }, "collapsed": false }, "outputs": [], "source": [ "timings = %timeit -o pairwise_python(X)\n", "pairwise_times['plain_python'] = timings" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### To samé pomocí NumPy\n", "V případě NumPy můžeme v tomto případě využít broadcasting. Celá funkce tak zabere jeden rádek." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2024-04-18T08:22:15.303807Z", "start_time": "2024-04-18T08:22:15.300664Z" }, "collapsed": false }, "outputs": [], "source": [ "def pairwise_numpy(X):\n", " return np.sqrt(((X[:, np.newaxis, :] - X) ** 2).sum(-1))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Zkusíme, jestli výsledky jsou stejné pomocí `assert` a `numpy.allclose`." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2024-04-18T08:22:19.469565Z", "start_time": "2024-04-18T08:22:17.419668Z" } }, "outputs": [], "source": [ "assert np.allclose(pairwise_numpy(X), pairwise_python(X), rtol=1e-10, atol=1e-15)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Výsledky jsou stejné až na velmi malé rozdíly - to je nebezpečí numerických výpočtů s konečnou přesností." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2024-04-18T08:22:23.531700Z", "start_time": "2024-04-18T08:22:21.409836Z" }, "collapsed": false }, "outputs": [], "source": [ "timings = %timeit -o pairwise_numpy(X)\n", "pairwise_times['numpy'] = timings" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Vidíme, že jsme zkrátili běh programu více než 100-krát. To není špatné, navíc je implementace daleko jednodušší." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Přichází Cython\n", "\n", "[Cython](https://cython.org/) 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.\n", "\n", "Zkusíme optimalizovat naší funkci `pairwise_python`.\n", "\n", "* Cython zdroják má koncovku .pyx (za začátku byl Pyrex).\n", "* Cython dokáže přeložit jakýkoli Python. Výsledkem je ale minimální (nebo spíš žádná) optimalizace.\n", "* `cimport` je analogie `import`, pracuje ale s Cython definicemi funkcí (.pxd soubory).\n", "* Cython dodává `numpy.pyx`, obsahující dodatečné informace pro kompilace NumPy modulů. Proto voláme `cimport numpy`.\n", "* Podobně `libc` je speciální modul Cythonu.\n", "\n", "* Funkce se deklarují (moho deklarovat) se statickými typy vstupních parametrů. My použijeme `np.ndarray[np.float64_t, ndim=2]`.\n", "* Proměnné se deklarují pomocí `cdef`.\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2024-04-18T08:22:30.642196Z", "start_time": "2024-04-18T08:22:30.639708Z" } }, "outputs": [], "source": [ "# Odkomentujte pro instalaci Cythonu\n", "# !pip install cython" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2024-04-18T08:22:32.130939Z", "start_time": "2024-04-18T08:22:32.127198Z" }, "collapsed": false }, "outputs": [], "source": [ "%%file cyfuncs.pyx\n", "\n", "language_level = \"3str\"\n", "\n", "import numpy as np\n", "# numpy pro Cython\n", "cimport numpy as np\n", "from libc.math cimport sqrt\n", "\n", "# tohle je čistý Python\n", "def pairwise0(X):\n", " M = X.shape[0]\n", " N = X.shape[1]\n", " D = np.empty((M, M), dtype=float)\n", " for i in range(M):\n", " for j in range(M):\n", " d = 0.0\n", " for k in range(N):\n", " tmp = X[i, k] - X[j, k]\n", " d += tmp * tmp\n", " D[i, j] = np.sqrt(d)\n", " return D\n", "\n", "# tady už začínáme optimalizovat, změny ale nejsou drastické\n", "def pairwise1(np.ndarray[np.float64_t, ndim=2] X):\n", " cdef int M = X.shape[0]\n", " cdef int N = X.shape[1]\n", " cdef double tmp, d\n", " cdef np.ndarray D = np.empty((M, M), dtype=np.float64)\n", " for i in range(M):\n", " for j in range(M):\n", " d = 0.0\n", " for k in range(N):\n", " tmp = X[i, k] - X[j, k]\n", " d += tmp * tmp\n", " D[i, j] = sqrt(d)\n", " return D" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2024-04-18T08:23:10.420625Z", "start_time": "2024-04-18T08:23:10.417392Z" }, "collapsed": false }, "outputs": [], "source": [ "%%file setup.py\n", "\n", "from distutils.core import setup\n", "from Cython.Build import cythonize\n", "import numpy\n", "\n", "setup(\n", " name='cyfuncs',\n", " include_dirs=[numpy.get_include()],\n", " ext_modules=cythonize(\"cyfuncs.pyx\"),\n", ")" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2024-04-18T08:23:13.073072Z", "start_time": "2024-04-18T08:23:10.984125Z" }, "collapsed": false }, "outputs": [], "source": [ "!python setup.py build_ext --inplace" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "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." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2024-04-18T08:23:17.343320Z", "start_time": "2024-04-18T08:23:17.339205Z" }, "collapsed": false }, "outputs": [], "source": [ "from IPython.display import FileLink\n", "FileLink('cyfuncs.c')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "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ř.\n", "\n", " /* \"cyfuncs.pyx\":16\n", " * tmp = X[i, k] - X[j, k]\n", " * d += tmp * tmp\n", " * D[i, j] = np.sqrt(d) # <<<<<<<<<<<<<<\n", " * return D\n", " * \n", " */\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2024-04-18T08:23:19.619108Z", "start_time": "2024-04-18T08:23:19.615667Z" }, "collapsed": false }, "outputs": [], "source": [ "import cyfuncs" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2024-04-18T08:23:20.331403Z", "start_time": "2024-04-18T08:23:20.328054Z" }, "collapsed": false }, "outputs": [], "source": [ "print(\"cyfuncs obsahuje: \" + \", \".join((d for d in dir(cyfuncs) if not d.startswith(\"_\"))))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Podívejme se, jestli dostávám stále stejné výsledky." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2024-04-18T08:23:28.253271Z", "start_time": "2024-04-18T08:23:28.160670Z" }, "collapsed": false }, "outputs": [], "source": [ "assert np.allclose(pairwise_numpy(X), cyfuncs.pairwise1(X), rtol=1e-10, atol=1e-15)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "No a jak jsme na tom s časem?" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2024-04-18T08:23:39.234005Z", "start_time": "2024-04-18T08:23:29.580814Z" }, "collapsed": false }, "outputs": [], "source": [ "timings = %timeit -o cyfuncs.pairwise0(X)\n", "pairwise_times['cython0'] = timings" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2024-04-18T08:23:44.674128Z", "start_time": "2024-04-18T08:23:39.235362Z" }, "collapsed": false }, "outputs": [], "source": [ "timings = %timeit -o cyfuncs.pairwise1(X)\n", "pairwise_times['cython1'] = timings" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### IPython `%%cython` magic\n", "\n", "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." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2024-04-18T08:23:45.671385Z", "start_time": "2024-04-18T08:23:45.405330Z" }, "collapsed": false }, "outputs": [], "source": [ "%load_ext Cython" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2024-04-18T08:23:49.316037Z", "start_time": "2024-04-18T08:23:45.830809Z" }, "collapsed": false }, "outputs": [], "source": [ "%%cython -a\n", "\n", "import numpy as np\n", "cimport numpy as np\n", "cimport cython\n", "from libc.math cimport sqrt\n", "\n", "@cython.boundscheck(False)\n", "@cython.wraparound(False)\n", "def pairwise_cython(double[:, ::1] X):\n", " cdef int M = X.shape[0]\n", " cdef int N = X.shape[1]\n", " cdef double tmp, d\n", " cdef double[:, ::1] D = np.empty((M, M), dtype=np.float64)\n", " for i in range(M):\n", " for j in range(M):\n", " d = 0.0\n", " for k in range(N):\n", " tmp = X[i, k] - X[j, k]\n", " d += tmp * tmp\n", " D[i, j] = sqrt(d)\n", " return np.asarray(D)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2024-04-18T08:24:09.786064Z", "start_time": "2024-04-18T08:24:09.723085Z" }, "collapsed": false }, "outputs": [], "source": [ "assert np.allclose(pairwise_numpy(X), pairwise_cython(X), rtol=1e-10, atol=1e-15)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2024-04-18T08:24:23.467435Z", "start_time": "2024-04-18T08:24:10.196820Z" }, "collapsed": false }, "outputs": [], "source": [ "timings = %timeit -o pairwise_cython(X)\n", "pairwise_times['cython2'] = timings" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Tohle je už konečně výrazné zrychlení oproti NumPy." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Cython toho nabízí mnoho\n", "Podívejte se na http://docs.cython.org co všechno Cython nabízí -- není toho málo, např.\n", "\n", "* použití C++\n", "* šablony (templates)\n", "* OpenMP (k tomu se možná ještě dostaneme)\n", "* vytváření C-API\n", "* třídy (cdef classes)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Z Fortranu (nebo C) do Pythonu pomocí F2PY\n", "\n", "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](http://cens.ioc.ee/projects/f2py2e/usersguide/index.html). 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í. \n", "\n", "Poďme chvilku programovat ve Fortranu :)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2024-04-18T08:25:02.527832Z", "start_time": "2024-04-18T08:25:02.523778Z" }, "collapsed": false }, "outputs": [], "source": [ "%%file pairwise_fort.f90\n", "\n", "subroutine pairwise_fort(X, D, m, n)\n", " integer :: n,m\n", " double precision, intent(in) :: X(m, n)\n", " double precision, intent(out) :: D(m, m)\n", " integer :: i, j, k\n", " double precision :: r\n", "\n", " do j = 1,m\n", " do i = 1,m\n", " r = 0\n", " do k = 1,n\n", " r = r + (X(i,k) - X(j,k)) * (X(i,k) - X(j,k))\n", " end do\n", " D(i,j) = sqrt(r)\n", " end do\n", " end do\n", "end subroutine pairwise_fort" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Z čeho f2py bere informace o vytvoření modulu?\n", "\n", "1. Pole (double precision) převádí na numpy array.\n", "2. `intent(in)` = vstupní argument.\n", "3. `intent(out)` = výstupní argument.\n", "4. f2py schová explicitně zadané rozměry polí (m, n).\n", "\n", "Pokud bychom programovali v C, je potřeba dodat f2py nějaké informace navíc, neboť např. intent v C neexistuje. \n", "\n", "Tento soubor přeložíme pomocí `f2py`:" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2024-04-18T08:25:05.699853Z", "start_time": "2024-04-18T08:25:04.807425Z" }, "collapsed": false }, "outputs": [], "source": [ "!f2py -c pairwise_fort.f90 -m pairwise_fort" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "`-m pairwise_fort` je požadované jméno modulu. Můžeme ho rovnou importovat, resp. jeho stejnojmennou funkci." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2024-04-18T08:25:11.219983Z", "start_time": "2024-04-18T08:25:11.216408Z" }, "collapsed": false }, "outputs": [], "source": [ "from pairwise_fort import pairwise_fort\n", "print(pairwise_fort.__doc__)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "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.\n", "\n", "Vytvoříme si proměnnou XF, která má Fortran uspořádání, pomocí `numpy.asfortranarray` (prozaický název :)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2024-04-18T08:25:15.329108Z", "start_time": "2024-04-18T08:25:15.326631Z" }, "collapsed": false }, "outputs": [], "source": [ "XF = np.asfortranarray(X)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Vyzkoušíme, jestli stále dostáváme stejné výsledky." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2024-04-18T08:25:17.031710Z", "start_time": "2024-04-18T08:25:16.997015Z" }, "collapsed": false }, "outputs": [], "source": [ "assert np.allclose(pairwise_numpy(X), pairwise_fort(X), rtol=1e-10, atol=1e-15)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "No a konečně se můžeme podívat, jak je to s rychlostí ..." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2024-04-18T08:25:34.844746Z", "start_time": "2024-04-18T08:25:20.063009Z" }, "collapsed": false }, "outputs": [], "source": [ "timings = %timeit -o pairwise_fort(X)\n", "pairwise_times['fortran'] = timings" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Představuje se `numba`\n", "\n", "Numba kompiluje Python kód pomocí [LLVM](http://llvm.org/). Podporujme just-in-time kompilaci pomocí dekorátoru `jit` (http://numba.pydata.org/numba-doc/latest/user/jit.html)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "```\n", "@numba.jit(\n", " signature=None, \n", " nopython=False, \n", " nogil=False, \n", " cache=False, \n", " forceobj=False, \n", " parallel=False, \n", " error_model='python', \n", " fastmath=False, locals={}\n", ")\n", "```" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2024-04-18T08:25:41.072952Z", "start_time": "2024-04-18T08:25:41.070431Z" } }, "outputs": [], "source": [ "# Odkomentujte pro instalaci balíku Numba\n", "# !pip install numba" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2024-04-18T08:25:59.325306Z", "start_time": "2024-04-18T08:25:59.202454Z" } }, "outputs": [], "source": [ "import numba" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2024-04-18T08:26:00.185673Z", "start_time": "2024-04-18T08:26:00.177083Z" } }, "outputs": [], "source": [ "pairwise_numba = numba.jit(pairwise_python)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Tradiční kontrola. Po prvním spuštění navíc Numba funkci poprvé zkompiluje." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2024-04-18T08:26:01.849016Z", "start_time": "2024-04-18T08:26:01.381482Z" }, "collapsed": false }, "outputs": [], "source": [ "assert np.allclose(pairwise_numpy(X), pairwise_numba(X), rtol=1e-10, atol=1e-15)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Jaký čas od takto \"jednoduché\" optimalizace můžeme očekávat?" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2024-04-18T08:26:05.311692Z", "start_time": "2024-04-18T08:26:03.651744Z" }, "collapsed": false }, "outputs": [], "source": [ "timings = %timeit -o pairwise_numba(X)\n", "pairwise_times['numba'] = timings" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "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ů)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "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`](http://numba.pydata.org/numba-doc/latest/user/parallel.html?highlight=prange). Všimněte si použití `@jit` jako dekorátoru." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2024-04-18T08:26:07.690637Z", "start_time": "2024-04-18T08:26:07.686997Z" } }, "outputs": [], "source": [ "@numba.jit(nopython=True, parallel=True, fastmath=True)\n", "def _pairwise_nopython(X: np.ndarray, D: np.ndarray) -> np.ndarray:\n", " M = X.shape[0]\n", " N = X.shape[1]\n", " for i in numba.prange(M):\n", " for j in numba.prange(M):\n", " d = 0.0\n", " for k in range(N):\n", " tmp = X[i, k] - X[j, k]\n", " d += tmp * tmp\n", " D[i, j] = np.sqrt(d)\n", " return D\n", "\n", "\n", "def pairwise_numba_fast_parallel(X: np.ndarray) -> np.ndarray:\n", " D = np.empty((X.shape[0], X.shape[0]), dtype = float)\n", " _pairwise_nopython(X, D)\n", " return D" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2024-04-18T08:26:08.828091Z", "start_time": "2024-04-18T08:26:08.498511Z" } }, "outputs": [], "source": [ "assert np.allclose(pairwise_numpy(X), pairwise_numba_fast_parallel(X), rtol=1e-10, atol=1e-15)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2024-04-18T08:26:13.446459Z", "start_time": "2024-04-18T08:26:08.927949Z" } }, "outputs": [], "source": [ "timings = %timeit -o pairwise_numba_fast_parallel(X)\n", "pairwise_times['numba_fast_parallel'] = timings" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Pojdme počítat na GPU i CPU: JAX\n", "\n", "_JAX: High-Performance Array Computing_ poskytuje plnou paletu nastrojů pro výpočty jak na CPU, tak na GPU. \n", "\n", "Mezi jeho hlavní výhody patří:\n", "\n", "* NumPy-like API - JAX podporuje většinu NumPy operací.\n", "* JIT kompilace - JAX umožňuje, podobně jako `numba`, kompilaci funkcí pomocí `jax.jit`.\n", "* Běží všude - stejný kód může běžet na CPU, GPU i TPU.\n", "* Automatická diferenciace - JAX umožňuje výpočet gradientů a Hessiánů.\n", "\n", " " ] }, { "cell_type": "code", "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2024-04-18T08:26:14.973233Z", "start_time": "2024-04-18T08:26:14.970674Z" } }, "outputs": [], "source": [ "# Odkomentujte pro instalaci balíku JAX\n", "# !pip install jax[cpu] jaxlib" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "JAX NumPy API je téměř identické s NumPy API. Jediný rozdíl je, že JAX funkce vrací vlastní datovy typ, které reprezentují pole na určitém zařízení (CPU, GPU, TPU).).\n", "Je proto vyhodné importovat JAX jako `import jax.numpy as jnp` společně s NumPy jako `import numpy as np`." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2024-04-18T08:26:16.447427Z", "start_time": "2024-04-18T08:26:15.990599Z" } }, "outputs": [], "source": [ "import jax.numpy as jnp\n", "import jax" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "V případě, že bychom měli dostupné GPU (AMD/NVIDIA), vyděli bychom jej v nasledujícím seznamu. " ] }, { "cell_type": "code", "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2024-04-18T08:26:18.876047Z", "start_time": "2024-04-18T08:26:18.859776Z" } }, "outputs": [], "source": [ "jax.devices()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Pro zpřístupnění vypočtů na GPU stačí nahradit `np` za `jnp` a JAX se postará o zbytek." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2024-04-18T08:26:21.002794Z", "start_time": "2024-04-18T08:26:20.999454Z" } }, "outputs": [], "source": [ "def pairwise_jax(X):\n", " return jnp.sqrt(jnp.power(X[:, jnp.newaxis, :] - X, 2).sum(-1))\n", " #return jnp.sqrt((jnp.power(X[:, jnp.newaxis, :] - X, 2)).sum(-1))\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2024-04-18T08:26:22.355890Z", "start_time": "2024-04-18T08:26:22.156958Z" } }, "outputs": [], "source": [ "assert np.allclose(pairwise_numpy(X), pairwise_jax(X), rtol=1e-10, atol=1e-15)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Ouha, co je špatně? Jelikož GPU pracuje s nižší přesností, JAX pracuje v zakladním režimu `float32` namísto `float64`." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2024-04-18T08:26:34.792567Z", "start_time": "2024-04-18T08:26:34.768471Z" } }, "outputs": [], "source": [ "pairwise_jax(X).dtype" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Pro vypočty s 64-bitovou přesností je potřeba povolit `jax_enable_x64`." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2024-04-18T08:26:36.962400Z", "start_time": "2024-04-18T08:26:36.959692Z" } }, "outputs": [], "source": [ "from jax import config\n", "config.update(\"jax_enable_x64\", True)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2024-04-18T08:26:40.204841Z", "start_time": "2024-04-18T08:26:40.108146Z" } }, "outputs": [], "source": [ "pairwise_jax(X).dtype" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2024-04-18T08:26:41.465030Z", "start_time": "2024-04-18T08:26:41.383496Z" } }, "outputs": [], "source": [ "assert np.allclose(pairwise_numpy(X), pairwise_jax(X), rtol=1e-10, atol=1e-15)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Proj JIT kompilaci stačí přidat dekorátor `jax.jit`. \n", "\n", "Nicméně pozor: \n", "* První volání funkce s JIT kompilací může byt pomalejší, než když bychom funkci volali bez kompilace.\n", "* JIT vyžaduje znát typy a tvar (shape) vstupních parametrů v době kompilace." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2024-04-18T08:26:43.629525Z", "start_time": "2024-04-18T08:26:43.626418Z" } }, "outputs": [], "source": [ "@jax.jit\n", "def pairwise_jax_jit(X):\n", " res = (jnp.sqrt((jnp.power(X[:, jnp.newaxis, :] - X, 2)))).sum(-1)\n", " return res\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2024-04-18T08:26:44.614037Z", "start_time": "2024-04-18T08:26:44.611318Z" } }, "outputs": [], "source": [ "X_jnp = np.asarray(X, dtype=jnp.float64)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2024-04-18T08:26:47.551431Z", "start_time": "2024-04-18T08:26:45.087129Z" } }, "outputs": [], "source": [ "timings = %timeit -o pairwise_jax(X_jnp)\n", "pairwise_times['jax'] = timings" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2024-04-18T08:26:55.285580Z", "start_time": "2024-04-18T08:26:50.089028Z" } }, "outputs": [], "source": [ "timings = %timeit -o pairwise_jax_jit(X_jnp)\n", "pairwise_times['jax_jit'] = timings" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Srovnání výsledků\n", "\n", "Výsledky můžeme porovnat pomocí grafu." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2024-04-18T08:26:59.171339Z", "start_time": "2024-04-18T08:26:59.167016Z" } }, "outputs": [], "source": [ "pairwise_times" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2024-04-18T08:27:01.532941Z", "start_time": "2024-04-18T08:27:01.208395Z" }, "collapsed": false }, "outputs": [], "source": [ "fig, ax = plt.subplots()\n", "values = np.array([t.average for t in pairwise_times.values()])\n", "x = range(len(pairwise_times))\n", "ax.bar(x, values)\n", "ax.set_xticks(x)\n", "ax.set_xticklabels(tuple(pairwise_times.keys()), rotation='vertical')\n", "ax.set_ylabel('time [ms]')\n", "ax.set_yscale('log')" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2024-04-18T08:31:41.312579Z", "start_time": "2024-04-18T08:31:41.309966Z" } }, "outputs": [], "source": [ "### Důkladnější srovnaní:\n", "pairwise_functions = {\n", " 'plain_python': pairwise_python,\n", " 'numpy': pairwise_numpy,\n", " 'cython0': cyfuncs.pairwise0,\n", " 'cython1': cyfuncs.pairwise1,\n", " 'cython2': pairwise_cython,\n", " 'fortran': pairwise_fort,\n", " 'numba': pairwise_numba,\n", " 'numba_fast_parallel': pairwise_numba_fast_parallel,\n", " 'jax': pairwise_jax,\n", " 'jax_jit': pairwise_jax_jit\n", "}" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2024-04-18T08:43:49.848782Z", "start_time": "2024-04-18T08:41:23.293147Z" } }, "outputs": [], "source": [ "\n", "Ms = np.logspace(1, 3.5, 6).astype(int)\n", "N = 3\n", "\n", "paiwise_times_all = {}\n", "for name, func in pairwise_functions.items():\n", " print(name)\n", " timings = []\n", " for M in Ms:\n", " X = np.random.random((M, N))\n", " print(M)\n", " # t = %timeit -o func(X)\n", " if \"jax\" in name:\n", " X_jnp = np.asarray(X, dtype=jnp.float64)\n", " # Přeskočíme první měření, protože probíha kompilace\n", " func(X_jnp)\n", " t = %timeit -o -r 2 func(X_jnp)\n", " else:\n", " t = %timeit -o -r 2 func(X)\n", "\n", " timings.append(t.average)\n", "\n", " paiwise_times_all[name] = timings\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "ExecuteTime": { "end_time": "2024-04-18T08:48:52.464845Z", "start_time": "2024-04-18T08:48:52.062914Z" } }, "outputs": [], "source": [ "fig, ax = plt.subplots()\n", "\n", "for name, times in paiwise_times_all.items():\n", " ax.plot(Ms, times, label=name)\n", " ax.set_xscale('log')\n", " ax.set_yscale('log')\n", "\n", "ax.set_xlabel('M')\n", "ax.set_ylabel('time [s]')\n", "\n", "ax.legend()\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Další možnosti\n", "\n", "* Vestavěný modul [`ctypes`](http://docs.python.org/2/library/ctypes.html) dovoluje volat funkce z externích dynamických knihoven. Pro použití s NumPy viz [Cookbook](http://wiki.scipy.org/Cookbook/Ctypes).\n", "* Alternativou k ctypes je [cffi](https://pypi.org/project/cffi/).\n", "* [CuPy](https://cupy.chainer.org/) využívá GPU.\n", "* [numexpr](https://github.com/pydata/numexpr) dokáže kompilovat Numpy výrazy.\n", "* [Theano](http://www.deeplearning.net/software/theano/index.html) se zaměřuje na strojové učení, také optimalizuje vektorové (maticové) operace, dovoluje je spouštět na CPU.\n", "* [Nuitka](http://nuitka.net/) kompiluje Python, ale na rozdíl od Numba nespecializuje funkce na základě typů.\n", "* [SWIG](http://www.swig.org/) jde použít pro propojení s mnoha jazyky.'" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Cvičení\n", "\n", "Obdobný postup aplikujte na výpočet kumulativního součtu, který je definovaný jako\n", "\n", "$\\displaystyle\n", "S_j = \\sum\\limits_{i = 1}^j x_i $\n", "\n", "Výsledky a časování porovnejte s `numpy.cumsum`.\n", "\n", "*Nápověda: Ve vaší funkci vytvořte nejprve výsledné numpy pole pomocí `numpy.empty_like`.*" ] } ], "metadata": { "anaconda-cloud": {}, "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.9.16" } }, "nbformat": 4, "nbformat_minor": 2 }