04-numba.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# numba\n",
    "\n",
    "numba es un compilador de Python a LLVM, que funciona para un subconjunto de la funcionalidad del lenguaje orientado a código científico en el que predomina la manipulación de arrays. Este notebook está basado en nuestro artículo [cómo acelerar código Python con numba](http://pybonacci.org/2015/03/13/como-acelerar-tu-codigo-python-con-numba/)\n",
    "\n",
    "En la web puedes ver la [evolución del rendimiento](http://numba.pydata.org/numba-benchmark/) de algunos benchmarks gracias a [asv](http://github.com/spacetelescope/asv/)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## ¿Cómo funciona?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "http://numba.pydata.org/numba-doc/0.18.2/developer/architecture.html"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Guion\n",
    "\n",
    "1. Ejemplo trivial: numba.jit\n",
    "2. Especificando tipos\n",
    "3. Código que no se acelera: numba.njit\n",
    "4. Placa plana"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Ejemplo trivial"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import numba"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "36.0\n"
     ]
    }
   ],
   "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\n",
    "\n",
    "a = np.arange(9).reshape(3,3)\n",
    "print(sum2d(a))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "big_arr = np.random.randn(1000, 1000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1 loops, best of 3: 232 ms per loop\n"
     ]
    }
   ],
   "source": [
    "%timeit sum2d(big_arr)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "sum2d_jit = numba.jit(sum2d)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "36.0"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sum2d_jit(a)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The slowest run took 69.00 times longer than the fastest. This could mean that an intermediate result is being cached \n",
      "1000 loops, best of 3: 1.15 ms per loop\n"
     ]
    }
   ],
   "source": [
    "%timeit sum2d_jit(big_arr)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Especificando tipos"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Al utilizar `numba.jit` estamos pidiendo a numba que *infiera* los tipos de las variables en el momento de la primera llamada.\n",
    "\n",
    "http://numba.pydata.org/numba-doc/0.18.2/reference/compilation.html#numba.jit\n",
    "\n",
    "Podemos ver esto utilizando el método `inspect_types` de la función resultante:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "sum2d (array(float64, 2d, C, nonconst),)\n",
      "--------------------------------------------------------------------------------\n",
      "# File: <ipython-input-2-884e457c72df>\n",
      "# --- LINE 1 --- \n",
      "\n",
      "def sum2d(arr):\n",
      "\n",
      "    # --- LINE 2 --- \n",
      "    # label 0\n",
      "    #   arr.1 = arr  :: array(float64, 2d, C, nonconst)\n",
      "    #   del arr\n",
      "    #   $0.2 = getattr(attr=shape, value=arr.1)  :: (int64 x 2)\n",
      "    #   $0.5 = exhaust_iter(count=2, value=$0.2)  :: (int64 x 2)\n",
      "    #   del $0.2\n",
      "    #   $0.3 = static_getitem(value=$0.5, index=0)  :: int64\n",
      "    #   $0.4 = static_getitem(value=$0.5, index=1)  :: int64\n",
      "    #   del $0.5\n",
      "    #   M = $0.3  :: int64\n",
      "    #   del $0.3\n",
      "    #   N = $0.4  :: int64\n",
      "    #   del $0.4\n",
      "\n",
      "    M, N = arr.shape\n",
      "\n",
      "    # --- LINE 3 --- \n",
      "    #   $const0.6 = const(float, 0.0)  :: float64\n",
      "    #   result = $const0.6  :: float64\n",
      "    #   del $const0.6\n",
      "\n",
      "    result = 0.0\n",
      "\n",
      "    # --- LINE 4 --- \n",
      "    #   jump 21\n",
      "    # label 34\n",
      "    #   $34.2 = iternext(value=$phi34.1)  :: pair<int64, bool>\n",
      "    #   $34.3 = pair_first(value=$34.2)  :: int64\n",
      "    #   $34.4 = pair_second(value=$34.2)  :: bool\n",
      "    #   del $34.2\n",
      "    #   $phi37.1 = $34.3  :: int64\n",
      "    #   del $34.3\n",
      "    #   branch $34.4, 37, 86\n",
      "    # label 21\n",
      "    #   $21.1 = global(range: <class 'range'>)  :: range\n",
      "    #   $21.3 = call $21.1(M, )  :: (int64,) -> range_state64\n",
      "    #   del M\n",
      "    #   del $21.1\n",
      "    #   $21.4 = getiter(value=$21.3)  :: range_iter64\n",
      "    #   del $21.3\n",
      "    #   $phi34.1 = $21.4  :: range_iter64\n",
      "    #   del $21.4\n",
      "    #   jump 34\n",
      "    # label 37\n",
      "    #   i = $phi37.1  :: int64\n",
      "    #   del $phi37.1\n",
      "\n",
      "    for i in range(M):\n",
      "\n",
      "        # --- LINE 5 --- \n",
      "        # label 56\n",
      "        #   j = $phi56.1  :: int64\n",
      "        #   del $phi56.1\n",
      "        # label 40\n",
      "        #   $40.1 = global(range: <class 'range'>)  :: range\n",
      "        #   $40.3 = call $40.1(N, )  :: (int64,) -> range_state64\n",
      "        #   del $40.1\n",
      "        #   $40.4 = getiter(value=$40.3)  :: range_iter64\n",
      "        #   del $40.3\n",
      "        #   $phi53.1 = $40.4  :: range_iter64\n",
      "        #   del $40.4\n",
      "        #   jump 53\n",
      "        # label 53\n",
      "        #   $53.2 = iternext(value=$phi53.1)  :: pair<int64, bool>\n",
      "        #   $53.3 = pair_first(value=$53.2)  :: int64\n",
      "        #   $53.4 = pair_second(value=$53.2)  :: bool\n",
      "        #   del $53.2\n",
      "        #   $phi56.1 = $53.3  :: int64\n",
      "        #   del $53.3\n",
      "        #   branch $53.4, 56, 82\n",
      "        #   jump 40\n",
      "\n",
      "        for j in range(N):\n",
      "\n",
      "            # --- LINE 6 --- \n",
      "            #   $56.6 = build_tuple(items=[Var(i, <ipython-input-2-884e457c72df> (4)), Var(j, <ipython-input-2-884e457c72df> (5))])  :: (int64 x 2)\n",
      "            #   del j\n",
      "            #   $56.7 = getitem(value=arr.1, index=$56.6)  :: float64\n",
      "            #   del $56.6\n",
      "            #   $56.8 = inplace_binop(rhs=$56.7, lhs=result, fn=+)  :: float64\n",
      "            #   del $56.7\n",
      "            #   result = $56.8  :: float64\n",
      "            #   del $56.8\n",
      "            #   jump 53\n",
      "            # label 83\n",
      "            #   jump 34\n",
      "            # label 82\n",
      "            #   del i\n",
      "            #   del $phi56.1\n",
      "            #   del $phi53.1\n",
      "            #   del $53.4\n",
      "            #   jump 83\n",
      "            # label 86\n",
      "            #   del arr.1\n",
      "            #   del N\n",
      "            #   del $phi37.1\n",
      "            #   del $phi34.1\n",
      "            #   del $34.4\n",
      "\n",
      "            result += arr[i,j]\n",
      "\n",
      "    # --- LINE 7 --- \n",
      "    #   jump 87\n",
      "    # label 87\n",
      "    #   $87.2 = cast(value=result)  :: float64\n",
      "    #   del result\n",
      "    #   return $87.2\n",
      "\n",
      "    return result\n",
      "\n",
      "\n",
      "================================================================================\n",
      "sum2d (array(int64, 2d, C, nonconst),)\n",
      "--------------------------------------------------------------------------------\n",
      "# File: <ipython-input-2-884e457c72df>\n",
      "# --- LINE 1 --- \n",
      "\n",
      "def sum2d(arr):\n",
      "\n",
      "    # --- LINE 2 --- \n",
      "    # label 0\n",
      "    #   arr.1 = arr  :: array(int64, 2d, C, nonconst)\n",
      "    #   del arr\n",
      "    #   $0.2 = getattr(attr=shape, value=arr.1)  :: (int64 x 2)\n",
      "    #   $0.5 = exhaust_iter(count=2, value=$0.2)  :: (int64 x 2)\n",
      "    #   del $0.2\n",
      "    #   $0.3 = static_getitem(value=$0.5, index=0)  :: int64\n",
      "    #   $0.4 = static_getitem(value=$0.5, index=1)  :: int64\n",
      "    #   del $0.5\n",
      "    #   M = $0.3  :: int64\n",
      "    #   del $0.3\n",
      "    #   N = $0.4  :: int64\n",
      "    #   del $0.4\n",
      "\n",
      "    M, N = arr.shape\n",
      "\n",
      "    # --- LINE 3 --- \n",
      "    #   $const0.6 = const(float, 0.0)  :: float64\n",
      "    #   result = $const0.6  :: float64\n",
      "    #   del $const0.6\n",
      "\n",
      "    result = 0.0\n",
      "\n",
      "    # --- LINE 4 --- \n",
      "    #   jump 21\n",
      "    # label 34\n",
      "    #   $34.2 = iternext(value=$phi34.1)  :: pair<int64, bool>\n",
      "    #   $34.3 = pair_first(value=$34.2)  :: int64\n",
      "    #   $34.4 = pair_second(value=$34.2)  :: bool\n",
      "    #   del $34.2\n",
      "    #   $phi37.1 = $34.3  :: int64\n",
      "    #   del $34.3\n",
      "    #   branch $34.4, 37, 86\n",
      "    # label 21\n",
      "    #   $21.1 = global(range: <class 'range'>)  :: range\n",
      "    #   $21.3 = call $21.1(M, )  :: (int64,) -> range_state64\n",
      "    #   del M\n",
      "    #   del $21.1\n",
      "    #   $21.4 = getiter(value=$21.3)  :: range_iter64\n",
      "    #   del $21.3\n",
      "    #   $phi34.1 = $21.4  :: range_iter64\n",
      "    #   del $21.4\n",
      "    #   jump 34\n",
      "    # label 37\n",
      "    #   i = $phi37.1  :: int64\n",
      "    #   del $phi37.1\n",
      "\n",
      "    for i in range(M):\n",
      "\n",
      "        # --- LINE 5 --- \n",
      "        # label 56\n",
      "        #   j = $phi56.1  :: int64\n",
      "        #   del $phi56.1\n",
      "        # label 40\n",
      "        #   $40.1 = global(range: <class 'range'>)  :: range\n",
      "        #   $40.3 = call $40.1(N, )  :: (int64,) -> range_state64\n",
      "        #   del $40.1\n",
      "        #   $40.4 = getiter(value=$40.3)  :: range_iter64\n",
      "        #   del $40.3\n",
      "        #   $phi53.1 = $40.4  :: range_iter64\n",
      "        #   del $40.4\n",
      "        #   jump 53\n",
      "        # label 53\n",
      "        #   $53.2 = iternext(value=$phi53.1)  :: pair<int64, bool>\n",
      "        #   $53.3 = pair_first(value=$53.2)  :: int64\n",
      "        #   $53.4 = pair_second(value=$53.2)  :: bool\n",
      "        #   del $53.2\n",
      "        #   $phi56.1 = $53.3  :: int64\n",
      "        #   del $53.3\n",
      "        #   branch $53.4, 56, 82\n",
      "        #   jump 40\n",
      "\n",
      "        for j in range(N):\n",
      "\n",
      "            # --- LINE 6 --- \n",
      "            #   $56.6 = build_tuple(items=[Var(i, <ipython-input-2-884e457c72df> (4)), Var(j, <ipython-input-2-884e457c72df> (5))])  :: (int64 x 2)\n",
      "            #   del j\n",
      "            #   $56.7 = getitem(value=arr.1, index=$56.6)  :: int64\n",
      "            #   del $56.6\n",
      "            #   $56.8 = inplace_binop(rhs=$56.7, lhs=result, fn=+)  :: float64\n",
      "            #   del $56.7\n",
      "            #   result = $56.8  :: float64\n",
      "            #   del $56.8\n",
      "            #   jump 53\n",
      "            # label 83\n",
      "            #   jump 34\n",
      "            # label 82\n",
      "            #   del i\n",
      "            #   del $phi56.1\n",
      "            #   del $phi53.1\n",
      "            #   del $53.4\n",
      "            #   jump 83\n",
      "            # label 86\n",
      "            #   del arr.1\n",
      "            #   del N\n",
      "            #   del $phi37.1\n",
      "            #   del $phi34.1\n",
      "            #   del $34.4\n",
      "\n",
      "            result += arr[i,j]\n",
      "\n",
      "    # --- LINE 7 --- \n",
      "    #   jump 87\n",
      "    # label 87\n",
      "    #   $87.2 = cast(value=result)  :: float64\n",
      "    #   del result\n",
      "    #   return $87.2\n",
      "\n",
      "    return result\n",
      "\n",
      "\n",
      "================================================================================\n"
     ]
    }
   ],
   "source": [
    "sum2d_jit.inspect_types()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Como hemos llamado a la función con arrays de enteros y de flotantes, vemos que se han generado las dos especificaciones. No obstante, nosotros podemos forzar los tipos de entrada:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "@numba.jit('float64(float64[:, :])')\n",
    "def sum2d_jit_float(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": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "36.0"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sum2d_jit_float(a.astype(float))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Algunas limitaciones"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Hay muchas partes de Python que no están implementadas en las versiones actuales de numba: en algunas de ellas se está trabajando activamente, mientras que otras es complicado que vean la luz en un futuro cercano.\n",
    "\n",
    "Una de estas limitaciones es que no podemos crear arrays nuevos dentro de las funciones. Si lo hacemos, numba no será capaz de optimizar todo el código e irá mucho más lento de lo que esperamos."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "@numba.jit\n",
    "def sum2d_jit_big():\n",
    "    arr = np.random.randn(1000, 1000)\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": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The slowest run took 4.81 times longer than the fastest. This could mean that an intermediate result is being cached \n",
      "1 loops, best of 3: 56 ms per loop\n"
     ]
    }
   ],
   "source": [
    "%timeit sum2d_jit_big()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1 loops, best of 3: 290 ms per loop\n"
     ]
    }
   ],
   "source": [
    "%timeit sum2d_jit_big.py_func()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "sum2d_jit_big ()\n",
      "--------------------------------------------------------------------------------\n",
      "# File: <ipython-input-11-ca6a1125999c>\n",
      "# --- LINE 1 --- \n",
      "\n",
      "@numba.jit\n",
      "\n",
      "# --- LINE 2 --- \n",
      "\n",
      "def sum2d_jit_big():\n",
      "\n",
      "    # --- LINE 3 --- \n",
      "    # label 0\n",
      "    #   $0.1 = global(np: <module 'numpy' from '/home/juanlu/.miniconda3/envs/py34/lib/python3.4/site-packages/numpy/__init__.py'>)  :: pyobject\n",
      "    #   $0.2 = getattr(attr=random, value=$0.1)  :: pyobject\n",
      "    #   del $0.1\n",
      "    #   $0.3 = getattr(attr=randn, value=$0.2)  :: pyobject\n",
      "    #   del $0.2\n",
      "    #   $const0.4 = const(int, 1000)  :: pyobject\n",
      "    #   $const0.5 = const(int, 1000)  :: pyobject\n",
      "    #   $0.6 = call $0.3($const0.4, $const0.5, )  :: pyobject\n",
      "    #   del $const0.5\n",
      "    #   del $const0.4\n",
      "    #   del $0.3\n",
      "    #   arr = $0.6  :: pyobject\n",
      "    #   del $0.6\n",
      "\n",
      "    arr = np.random.randn(1000, 1000)\n",
      "\n",
      "    # --- LINE 4 --- \n",
      "    #   $0.8 = getattr(attr=shape, value=arr)  :: pyobject\n",
      "    #   $0.11 = exhaust_iter(count=2, value=$0.8)  :: pyobject\n",
      "    #   del $0.8\n",
      "    #   $0.9 = static_getitem(value=$0.11, index=0)  :: pyobject\n",
      "    #   $0.10 = static_getitem(value=$0.11, index=1)  :: pyobject\n",
      "    #   del $0.11\n",
      "    #   M = $0.9  :: pyobject\n",
      "    #   del $0.9\n",
      "    #   N = $0.10  :: pyobject\n",
      "    #   del $0.10\n",
      "\n",
      "    M, N = arr.shape\n",
      "\n",
      "    # --- LINE 5 --- \n",
      "    #   $const0.12 = const(float, 0.0)  :: pyobject\n",
      "    #   result = $const0.12  :: pyobject\n",
      "    #   del $const0.12\n",
      "\n",
      "    result = 0.0\n",
      "\n",
      "    # --- LINE 6 --- \n",
      "    #   jump 45.1\n",
      "    # label 45.1\n",
      "    #   $const45.1.1 = const(LiftedLoop, LiftedLoop(<function sum2d_jit_big at 0x7fc3a06b8950>))  :: pyobject\n",
      "    #   $45.1.6 = call $const45.1.1(arr, M, result, N, )  :: pyobject\n",
      "    #   del result\n",
      "    #   del arr\n",
      "    #   del N\n",
      "    #   del M\n",
      "    #   del $const45.1.1\n",
      "    #   $45.1.8 = exhaust_iter(count=1, value=$45.1.6)  :: pyobject\n",
      "    #   del $45.1.6\n",
      "    #   $45.1.7 = static_getitem(value=$45.1.8, index=0)  :: pyobject\n",
      "    #   del $45.1.8\n",
      "    #   result.1 = $45.1.7  :: pyobject\n",
      "    #   del $45.1.7\n",
      "    #   jump 108\n",
      "\n",
      "    for i in range(M):\n",
      "\n",
      "        # --- LINE 7 --- \n",
      "\n",
      "        for j in range(N):\n",
      "\n",
      "            # --- LINE 8 --- \n",
      "\n",
      "            result += arr[i,j]\n",
      "\n",
      "    # --- LINE 9 --- \n",
      "    # label 108\n",
      "    #   $108.2 = cast(value=result.1)  :: pyobject\n",
      "    #   del result.1\n",
      "    #   return $108.2\n",
      "\n",
      "    return result\n",
      "\n",
      "# The function contains lifted loops\n",
      "# Loop at line 6\n",
      "# Has 1 overloads\n",
      "# File: <ipython-input-11-ca6a1125999c>\n",
      "# --- LINE 1 --- \n",
      "\n",
      "@numba.jit\n",
      "\n",
      "# --- LINE 2 --- \n",
      "\n",
      "def sum2d_jit_big():\n",
      "\n",
      "    # --- LINE 3 --- \n",
      "\n",
      "    arr = np.random.randn(1000, 1000)\n",
      "\n",
      "    # --- LINE 4 --- \n",
      "\n",
      "    M, N = arr.shape\n",
      "\n",
      "    # --- LINE 5 --- \n",
      "\n",
      "    result = 0.0\n",
      "\n",
      "    # --- LINE 6 --- \n",
      "    # label 58\n",
      "    #   i = $phi58.1  :: int64\n",
      "    #   del $phi58.1\n",
      "    # label 55\n",
      "    #   $55.2 = iternext(value=$phi55.1)  :: pair<int64, bool>\n",
      "    #   $55.3 = pair_first(value=$55.2)  :: int64\n",
      "    #   $55.4 = pair_second(value=$55.2)  :: bool\n",
      "    #   del $55.2\n",
      "    #   $phi58.1 = $55.3  :: int64\n",
      "    #   del $55.3\n",
      "    #   branch $55.4, 58, 107\n",
      "    # label 42\n",
      "    #   arr.1 = arr  :: array(float64, 2d, C, nonconst)\n",
      "    #   del arr\n",
      "    #   M.1 = M  :: int64\n",
      "    #   del M\n",
      "    #   result.1 = result  :: float64\n",
      "    #   del result\n",
      "    #   N.1 = N  :: int64\n",
      "    #   del N\n",
      "    #   $42.1 = global(range: <class 'range'>)  :: range\n",
      "    #   $42.3 = call $42.1(M.1, )  :: (int64,) -> range_state64\n",
      "    #   del M.1\n",
      "    #   del $42.1\n",
      "    #   $42.4 = getiter(value=$42.3)  :: range_iter64\n",
      "    #   del $42.3\n",
      "    #   $phi55.1 = $42.4  :: range_iter64\n",
      "    #   del $42.4\n",
      "    #   jump 55\n",
      "\n",
      "    for i in range(M):\n",
      "\n",
      "        # --- LINE 7 --- \n",
      "        #   jump 61\n",
      "        # label 74\n",
      "        #   $74.2 = iternext(value=$phi74.1)  :: pair<int64, bool>\n",
      "        #   $74.3 = pair_first(value=$74.2)  :: int64\n",
      "        #   $74.4 = pair_second(value=$74.2)  :: bool\n",
      "        #   del $74.2\n",
      "        #   $phi77.1 = $74.3  :: int64\n",
      "        #   del $74.3\n",
      "        #   branch $74.4, 77, 103\n",
      "        # label 61\n",
      "        #   $61.1 = global(range: <class 'range'>)  :: range\n",
      "        #   $61.3 = call $61.1(N.1, )  :: (int64,) -> range_state64\n",
      "        #   del $61.1\n",
      "        #   $61.4 = getiter(value=$61.3)  :: range_iter64\n",
      "        #   del $61.3\n",
      "        #   $phi74.1 = $61.4  :: range_iter64\n",
      "        #   del $61.4\n",
      "        #   jump 74\n",
      "        # label 77\n",
      "        #   j = $phi77.1  :: int64\n",
      "        #   del $phi77.1\n",
      "\n",
      "        for j in range(N):\n",
      "\n",
      "            # --- LINE 8 --- \n",
      "            # label 107\n",
      "            #   del arr.1\n",
      "            #   del N.1\n",
      "            #   del $phi58.1\n",
      "            #   del $phi55.1\n",
      "            #   del $55.4\n",
      "            #   jump 108\n",
      "            # label 104\n",
      "            #   jump 55\n",
      "            # label 103\n",
      "            #   del i\n",
      "            #   del $phi77.1\n",
      "            #   del $phi74.1\n",
      "            #   del $74.4\n",
      "            #   jump 104\n",
      "            # label 108\n",
      "            #   $108.2 = build_tuple(items=[Var(result.1, <ipython-input-11-ca6a1125999c> (6))])  :: (float64 x 1)\n",
      "            #   del result.1\n",
      "            #   $108.3 = cast(value=$108.2)  :: (float64 x 1)\n",
      "            #   del $108.2\n",
      "            #   return $108.3\n",
      "            #   $77.6 = build_tuple(items=[Var(i, <ipython-input-11-ca6a1125999c> (6)), Var(j, <ipython-input-11-ca6a1125999c> (7))])  :: (int64 x 2)\n",
      "            #   del j\n",
      "            #   $77.7 = getitem(value=arr.1, index=$77.6)  :: float64\n",
      "            #   del $77.6\n",
      "            #   $77.8 = inplace_binop(rhs=$77.7, lhs=result.1, fn=+)  :: float64\n",
      "            #   del $77.7\n",
      "            #   result.1 = $77.8  :: float64\n",
      "            #   del $77.8\n",
      "            #   jump 74\n",
      "\n",
      "            result += arr[i,j]\n",
      "\n",
      "    # --- LINE 9 --- \n",
      "\n",
      "    return result\n",
      "\n",
      "\n",
      "\n",
      "================================================================================\n"
     ]
    }
   ],
   "source": [
    "sum2d_jit_big.inspect_types()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Al contrario de lo que teníamos antes, se han creado numerosas variables de tipo `pyobject`. Esto sucede en varios casos:\n",
    "\n",
    "* El algoritmo de inferencia de tipos falla\n",
    "* Utilizamos características del lenguaje no soportadas por numba\n",
    "* Utilizamos objetos de Python complejos"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Entendiendo numba: el modo *nopython*\n",
    "\n",
    "Como podemos leer en la documentación, [numba tiene dos modos de funcionamiento básicos](http://numba.pydata.org/numba-doc/0.17.0/user/jit.html#nopython): el modo *object* y el modo *nopython*.\n",
    "\n",
    "* El modo *object* genera código que gestiona todas las variables como objetos de Python y utiliza la API C de Python para operar con ellas. En general en este modo **no habrá ganancias de rendimiento** (e incluso puede ir más lento), con lo cual mi recomendación personal es directamente no utilizarlo. Hay casos en los que numba puede detectar los bucles y optimizarlos individualmente (*loop-jitting*), pero no le voy a prestar mucha atención a esto.\n",
    "* El modo *nopython* **genera código independiente de la API C de Python**. Esto tiene la desventaja de que no podemos usar todas las características del lenguaje, **pero tiene un efecto significativo en el rendimiento**. Otra de las restricciones es que **no se puede reservar memoria para objetos nuevos**.\n",
    "\n",
    "Por defecto numba usa el modo *nopython* siempre que puede, y si no pasa a modo *object*. Nosotros vamos a **forzar el modo nopython** (o «modo estricto» como me gusta llamarlo) porque es la única forma de aprovechar el potencial de numba."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Ámbito de aplicación\n",
    "\n",
    "El problema del modo *nopython* es que los mensajes de error son totalmente inservibles en la mayoría de los casos, así que antes de lanzarnos a compilar funciones con numba conviene hacer un repaso de qué no podemos hacer para anticipar la mejor forma de programar nuestro código. Podéis consultar en la documentación [el subconjunto de Python soportado por numba](http://numba.pydata.org/numba-doc/0.17.0/reference/pysupported.html) en modo *nopython*, y ya os aviso que, al menos de momento, no tenemos [*list comprehensions*](https://github.com/numba/numba/issues/504), [generadores](https://github.com/numba/numba/issues/984) ni algunas cosas más. Permitidme que resalte una frase sacada de la página principal de numba:\n",
    "\n",
    "> \"*With a few annotations, **array-oriented and math-heavy Python code** can be just-in-time compiled to native machine instructions, similar in performance to C, C++ and Fortran*\". [Énfasis mío]\n",
    "\n",
    "Siento decepcionar a la audiencia pero *numba no acelerará todo el código Python* que le echemos: está enfocado a operaciones matemáticas con arrays. Aclarado este punto, vamos a ponernos manos a la obra con un ejemplo aplicado :)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Solución de Navier de una placa plana"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from numpy import sin, pi\n",
    "\n",
    "def plate_displacement_py(xx, yy, ww, a, b, P, xi, eta, D, max_i, max_j, max_m, max_n):\n",
    "    for mm in range(1, max_m):\n",
    "        for nn in range(1, max_n):\n",
    "            for ii in range(max_i):\n",
    "                for jj in range(max_j):\n",
    "                    a_mn = 4 * P * sin(mm * pi * xi / a) * sin(nn * pi * eta / b) / (a * b)\n",
    "                    ww[ii, jj] += (a_mn / (mm**2 / a**2 + nn**2 / b**2)**2\n",
    "                                   * sin(mm * pi * xx[ii, jj] / a)\n",
    "                                   * sin(nn * pi * yy[ii, jj] / b)\n",
    "                                   / (pi**4 * D))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "plate_displacement_nb = numba.jit(plate_displacement_py, nopython=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Plate geometry\n",
    "a = 1.0  # m\n",
    "b = 1.0  # m\n",
    "h = 50e-3  # m\n",
    "\n",
    "# Material properties\n",
    "E = 69e9  # Pa\n",
    "nu = 0.35\n",
    "\n",
    "# Series terms\n",
    "max_m = 16\n",
    "max_n = 16\n",
    "\n",
    "# Computation points\n",
    "# NOTE: With an odd number of points the center of the place is included in\n",
    "# the grid\n",
    "NUM_POINTS = 101\n",
    "max_i = NUM_POINTS\n",
    "max_j = NUM_POINTS\n",
    "\n",
    "# Load\n",
    "P = 10e3  # N\n",
    "xi = a / 2\n",
    "eta = a / 2\n",
    "\n",
    "# Flexural rigidity\n",
    "D = h**3 * E / (12 * (1 - nu**2))\n",
    "\n",
    "# ---\n",
    "\n",
    "# Set up domain\n",
    "x = np.linspace(0, a, num=NUM_POINTS)\n",
    "y = np.linspace(0, b, num=NUM_POINTS)\n",
    "xx, yy = np.meshgrid(x, y)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Maximum displacement = 0.141317840389 mm\n",
      "alpha = 0.01158\n",
      "alpha * P a^2 / D = 0.1416 mm\n"
     ]
    }
   ],
   "source": [
    "# Compute displacement field\n",
    "ww = np.zeros_like(xx)\n",
    "plate_displacement_nb(xx, yy, ww, a, b, P, xi, eta, D, max_i, max_j, max_m, max_n)\n",
    "\n",
    "# Print maximum displacement\n",
    "w_max = np.abs(ww).max()\n",
    "print(\"Maximum displacement = %14.12f mm\" % (w_max * 1e3))\n",
    "print(\"alpha = %7.5f\" % (w_max / (P * a**2 / D)))\n",
    "print(\"alpha * P a^2 / D = %6.4f mm\" % (0.01160 * P * a**2 / D * 1e3))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1 loops, best of 3: 1.13 s per loop\n"
     ]
    }
   ],
   "source": [
    "%timeit plate_displacement_nb(xx, yy, ww, a, b, P, xi, eta, D, max_i, max_j, max_m, max_n)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Conseguimos multiplicar por 100 el rendimiento**, aunque como ya hemos visto este no es nuestro mejor resultado.\n",
    "\n",
    "Si probamos la versión de NumPy las cosas no mejoran respecto al original, porque se lanza el modo `object`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def plate_displacement_np(xx, yy, ww, a, b, P, xi, eta, D, max_m, max_n):\n",
    "    for mm in range(1, max_m):\n",
    "        for nn in range(1, max_n):\n",
    "            a_mn = 4 * P * sin(mm * pi * xi / a) * sin(nn * pi * eta / b) / (a * b)\n",
    "            ww += (a_mn / (mm**2 / a**2 + nn**2 / b**2)**2\n",
    "                   * sin(mm * pi * xx / a)\n",
    "                   * sin(nn * pi * yy / b)\n",
    "                   / (pi**4 * D))\n",
    "\n",
    "plate_displacement_nb2 = numba.jit(plate_displacement_np)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Maximum displacement = 0.141317840389 mm\n",
      "alpha = 0.01158\n",
      "alpha * P a^2 / D = 0.1416 mm\n"
     ]
    }
   ],
   "source": [
    "# Compute displacement field\n",
    "ww = np.zeros_like(xx)\n",
    "plate_displacement_nb2(xx, yy, ww, a, b, P, xi, eta, D, max_m, max_n)\n",
    "\n",
    "# Print maximum displacement\n",
    "w_max = np.abs(ww).max()\n",
    "print(\"Maximum displacement = %14.12f mm\" % (w_max * 1e3))\n",
    "print(\"alpha = %7.5f\" % (w_max / (P * a**2 / D)))\n",
    "print(\"alpha * P a^2 / D = %6.4f mm\" % (0.01160 * P * a**2 / D * 1e3))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1 loops, best of 3: 205 ms per loop\n"
     ]
    }
   ],
   "source": [
    "%timeit plate_displacement_nb2(xx, yy, ww, a, b, P, xi, eta, D, max_m, max_n)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**El rendimiento es igual al de NumPy, sin más**."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Conclusiones"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "numba es una herramienta muy poderosa con la que podemos exprimir el rendimiento del código haciendo muy pocos cambios y todavía escribiendo en un único lenguaje. Hay otras características poderosas, como liberación del GIL o aceleración en GPUs, que no vamos a ver aquí. No obstante, aunque la API se ha estabilizado bastante en las últimas versiones aún quedan algunas regresiones de rendimiento por encontrar y características del lenguaje por implementar. El desarrollo es muy rápido así que seguiremos viendo novedades en este frente durante los próximos meses."
   ]
  }
 ],
 "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.4.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}