7-optimizacion.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "\n",
    "# Acelerando Python\n",
    "\n",
    "(Notebook tomada de Introducción a Python para ciencias e ingenierías (clase 6) - Ing. Martín Gaitán)\n",
    "\n",
    "Algunas cosas que hay que saber:\n",
    "\n",
    "- Python es un lenguaje de programación estandarizado que tiene **múltiples implementaciones**\n",
    "- La implementación típica y oficial se llama estrictamente **CPython**, porque está programada en **C**\n",
    "- Existen **otras implementaciones** del estándar: IronPython (.Net), Jython (Java) y hasta una implementación en una versión reducida Python llama **PyPy** (que, de paso, es rapidísima)\n",
    "- Cada versión tiene sus propias ventajas, pero el hecho de que CPython sea C **permite extender con código muy eficiente hecho en C** (o C++)\n",
    "- Pero **programar** en lenguajes de bajo nivel como **C es un *perno* **. \n",
    "\n",
    "\n",
    "## Repaso express de f2py\n",
    "\n",
    "\n",
    "Ya vimos la herramienta F2py que aprovecha que Fortran y CPython se basan en C, y analizando el código fortran genera el \"pegamento\" (cabecesas y código C) para comunicar Python con Fortran. El resultado es un módulo (una extensión) importable desde python, donde las subrutinas fortran se convierten en funciones Python. \n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%file suma.f90 \n",
    "\n",
    "subroutine suma2(n, v)\n",
    "    ! Compute the dot product between u and v (length n) and put the result in product\n",
    "     \n",
    "    integer, intent(in) :: n\n",
    "    integer, intent(out) :: v\n",
    "    v = n + 2\n",
    "    \n",
    "end subroutine"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Si queres evitar el output podés guardar en una variable\n",
    "_ = !f2py -m suma -c suma.f90   \n",
    "#!/home/tin/.virtualenvs/curso/bin/f2py3.4 -m suma -c suma.f90 "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import suma\n",
    "suma.suma2(10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%install_ext https://raw.github.com/mgaitan/fortran_magic/master/fortranmagic.py"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%load_ext fortranmagic"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%fortran\n",
    "\n",
    "subroutine f1(x, y, z)\n",
    "    real, intent(in) :: x,y\n",
    "    real, intent(out) :: z\n",
    "\n",
    "    z = sin(x+y)\n",
    "\n",
    "end subroutine f1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f1(4, 2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Cython\n",
    "\n",
    "Cython es conceptualmente similar: **es un lenguaje** (con su correspondiente compilador), que extiende la sintaxis de Python (python con sabor a C). El **compilador Cython produce código C estándar** listo para ser compilado y usado como una extensión.\n",
    "\n",
    "* código python válido es código cython válido\n",
    "* Se pueden **llamar funciones en C**, o funciones/métodos de C++, **directamente** desde el código en Cython.   \n",
    "* Es posible usar declarar tipos explícitamente (enteros, flotantes, o cualquier tipo de dato).\n",
    "\n",
    "### Instalación\n",
    "\n",
    "\n",
    "    conda install cython\n",
    "\n",
    "o \n",
    "\n",
    "    pip install cython\n",
    "\n",
    "\n",
    "### Hola mundo (versión larga)\n",
    "\n",
    "Por convención, los archivos Cython se guardan con extensión .pyx  (viene de Pyrex, que es el paquete precedente)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "%%file helloworld.pyx\n",
    "\n",
    "print('Hola Curso, soy Cython. #NoFueMagia')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Podemos compilar este código Cython para que genere el código C necesario para crear nuestra extensión Python"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "!cython -a helloworld.pyx"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "!cat helloworld.c"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from IPython.core.display import HTML\n",
    "HTML(filename='helloworld.html')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Con el código C necesario, podemos compilar nuestra extensión"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!gcc -shared -pthread -fPIC -fwrapv -O2 -Wall -fno-strict-aliasing -I /usr/include/python3.5m -o helloworld.so helloworld.c"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import helloworld"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Hola mundo, the magic way\n",
    "\n",
    "Pero si lo usamos en el notebook y lo modificamos seguido, es mejor usa el comándo magic **cython** "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# %load_ext cythonmagic    en versiones anteriores a 0.21\n",
    "%load_ext cython"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%cython \n",
    "\n",
    "def hola():\n",
    "    print('hola, soy una funcion cython')\n",
    "\n",
    "print('Hola, soy cython, a través de un magic.')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hola()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "### Paréntesis: funciones generadoras"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def fib(n):\n",
    "    \"\"\"cálcula la serie hasta n y devuelve una lista\"\"\"\n",
    "    serie = []\n",
    "    a, b = 0, 1\n",
    "    while b < n:\n",
    "        serie.append(b)\n",
    "        a, b = b, a + b\n",
    "    return serie\n",
    "\n",
    "fib(120)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%timeit fib(1000000)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### El mismo generador en Cython"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%%cython\n",
    "\n",
    "def fib_cython(n):\n",
    "    \"\"\"Print the Fibonacci series up to n.\"\"\"\n",
    "    a, b = 0, 1\n",
    "    while b < n:\n",
    "        yield b\n",
    "        a, b = b, a + b"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%timeit list(fib_cython(100000))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "### Cython on steroids\n",
    "\n",
    "Pero veamos un código Cython que no sea Python puro. Supongamos que queremos calcular $\\sum_{i=0}^n \\alpha^i$ para $\\alpha, n$ dados. \n",
    "\n",
    "Obviamente, antes de optimizar un algoritmo podriamos aplicar la formula de la progresión geométrica \n",
    "\n",
    "$$\n",
    "\\sum_{i=0}^n \\alpha^i = \\frac{1 - \\alpha^{n+1}}{1 - \\alpha}\n",
    "$$\n",
    "\n",
    "pero para \"medir músculos\" (?) lo haremos iterativamente"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def geo_prog(alpha, n):\n",
    "    current = 1.0\n",
    "    sum = current\n",
    "    for i in range(n):\n",
    "        current = current * alpha\n",
    "        sum = sum + current\n",
    "    return sum"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%timeit geo_prog(-4.0, 1000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%%cython \n",
    "def geo_prog_cython(double alpha, int n):\n",
    "    cdef double current = 1.0\n",
    "    cdef double sum = current          # declaracion estática de tipo!\n",
    "    cdef int i\n",
    "    for i in range(n):\n",
    "        current = current * alpha\n",
    "        sum = sum + current\n",
    "    return sum"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%timeit geo_prog_cython(-4.0, 1000)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Numba\n",
    "\n",
    "Numba es otra manera de optimizar código Python a través de un proceso de compilación. A diferencia de Cython, el proceso de compilación no produce código C válido para ser compilado como una extensión, sino que directamente se compila a \"código maquina\", pero *justo cuando se lo necesita* (ver [Just in Time compilation](http://en.wikipedia.org/wiki/Just-in-time_compilation))\n",
    "\n",
    "Numba utiliza como tecnología subyacente la infraestructura de compilación [LLVM](http://en.wikipedia.org/wiki/LLVM) que es la encargada de realizar estas optimizaciones al vuelo \n",
    "\n",
    "### Instalación\n",
    "\n",
    "Usando conda\n",
    "\n",
    "    conda install numba\n",
    "    \n",
    "Via pip/virtualenv (linux)\n",
    "\n",
    "    sudo apt-get install libllvm build-essential libllvm-dev libedit-dev\n",
    "    pip install llvmlite\n",
    "    pip install numba\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def sum2d(arr):\n",
    "    M, N = arr.shape\n",
    "    result = 0.0\n",
    "    for i in range(M):\n",
    "        for j in range(N):\n",
    "            result += arr[i,j]\n",
    "    return result"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "x = np.ones((10, 20))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sum2d(x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "m_random = np.random.rand(1000, 1000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sum2d(m_random)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%timeit sum2d(m_random)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Así se aplica la magia de Numba"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from numba import jit\n",
    "\n",
    "sum2d_numba = jit(sum2d)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Para quienes desconfian, comprobemos que la funcion \"numbizada\" produce el mismo resultado"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sum2d_numba(m_random)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Pero anda ~10$^2$ más rápido"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%timeit sum2d_numba(m_random)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Una forma forma más elegante de usar usar `jit` es como un **decorador**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "@jit\n",
    "def sum2d_(arr):\n",
    "    M, N = arr.shape\n",
    "    result = 0.0\n",
    "    for i in range(M):\n",
    "        for j in range(N):\n",
    "            result += arr[i,j]\n",
    "    return result"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%timeit sum2d_(m_random)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "¡Más de 100 veces más rápido con una línea de código!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Los límites de la magia de Numba y como maximizarla\n",
    "\n",
    "Cuando no se definen tipos en la signatura, Numba hace una compilación *\"lazy\"* (vaga, que no se produce hasta el momento de se invoncada la función), e intentará operar sobre cualquier tipo de datos que reciba"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from numba import jit\n",
    "\n",
    "@jit\n",
    "def plus(x, y):\n",
    "    # A somewhat trivial example\n",
    "    return x + y"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plus(3.4, 3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plus('a', 'b')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Cuando especificamos la signatura, la optimización es *\"eager\"*, y obviamente esta especialización aumenta la performance sacrificando versatilidad"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "@jit(int32(int32, int32))\n",
    "def plus32(x, y):\n",
    "    return x + y"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Numba tiene **dos estrategias de compilación**: una optimización directa y de máxima perfomance denominado **modo nopython**, que genera código maquina que \"no requiere consultar a la API C de Python\" logrando inferir los tipos de datos con los que está trabajando. Cuando esta inferencia falla (cuando Numba no logró inferir cuál es el tipo de dato \"dinámico\" de alguna operación), automáticamente conmuta al **modo object** que sí utiliza la API C de Python. En este modo las únicas optimizaciones que se producen son con los *loops*. \n",
    "\n",
    "\n",
    "\n",
    "Pero la técnica de optimización es **factorizar el código \"nopython\"**. Para evitar que implicitamente se pase a modo object, se puede forzar el modo con `nopython=True`\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "@jit(nopython=True)\n",
    "def suma_nopy(x, y):\n",
    "    return x + y"
   ]
  }
 ],
 "metadata": {
  "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.6.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}