| Software | Version |
|---|
| Python | 3.4.3 64bit [GCC 4.4.7 20120313 (Red Hat 4.4.7-1)] |
| IPython | 3.0.0 |
| OS | Linux 3.18.6 1 ARCH x86_64 with arch |
| numpy | 1.9.2 |
| numba | 0.17.0 |
| cython | 0.22 |
| Fri Mar 13 13:44:39 2015 CET |
"
],
"text/latex": [
"\\begin{tabular}{|l|l|}\\hline\n",
"{\\bf Software} & {\\bf Version} \\\\ \\hline\\hline\n",
"Python & 3.4.3 64bit [GCC 4.4.7 20120313 (Red Hat 4.4.7-1)] \\\\ \\hline\n",
"IPython & 3.0.0 \\\\ \\hline\n",
"OS & Linux 3.18.6 1 ARCH x86\\_64 with arch \\\\ \\hline\n",
"numpy & 1.9.2 \\\\ \\hline\n",
"numba & 0.17.0 \\\\ \\hline\n",
"cython & 0.22 \\\\ \\hline\n",
"\\hline \\multicolumn{2}{|l|}{Fri Mar 13 13:44:39 2015 CET} \\\\ \\hline\n",
"\\end{tabular}\n"
],
"text/plain": [
"Software versions\n",
"Python 3.4.3 64bit [GCC 4.4.7 20120313 (Red Hat 4.4.7-1)]\n",
"IPython 3.0.0\n",
"OS Linux 3.18.6 1 ARCH x86_64 with arch\n",
"numpy 1.9.2\n",
"numba 0.17.0\n",
"cython 0.22\n",
"Fri Mar 13 13:44:39 2015 CET"
]
},
"execution_count": 3,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"%version_information numpy, numba, cython"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Vamos a empezar por importar los paquetes necesarios e inicializar la semilla del generador de números aleatorios:"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"import numpy as np\n",
"import numba\n",
"\n",
"np.random.seed(0)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Creamos nuestro array de datos:"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"data = np.random.randn(2000, 2000)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Y voy a copiar descaradamente la función de Kiko:"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"def busca_min(malla):\n",
" minimosx = []\n",
" minimosy = []\n",
" for i in range(1, malla.shape[1]-1):\n",
" for j in range(1, malla.shape[0]-1):\n",
" if (malla[j, i] < malla[j-1, i-1] and\n",
" malla[j, i] < malla[j-1, i] and\n",
" malla[j, i] < malla[j-1, i+1] and\n",
" malla[j, i] < malla[j, i-1] and\n",
" malla[j, i] < malla[j, i+1] and\n",
" malla[j, i] < malla[j+1, i-1] and\n",
" malla[j, i] < malla[j+1, i] and\n",
" malla[j, i] < malla[j+1, i+1]):\n",
" minimosx.append(i)\n",
" minimosy.append(j)\n",
"\n",
" return np.array(minimosx), np.array(minimosy)"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"(array([ 1, 1, 1, ..., 1998, 1998, 1998]),\n",
" array([ 1, 3, 11, ..., 1968, 1977, 1985]))"
]
},
"execution_count": 7,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"busca_min(data)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Paso 1: analizar el código\n",
"\n",
"Lo primero que pensé cuando vi esta función es que no me gustaba nada hacer `append` a esas dos listas tantas veces. Pero a continuación me pregunté si realmente tendrían tantos elementos... averigüémoslo:"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"0.11091025"
]
},
"execution_count": 8,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"mx, my = busca_min(data)\n",
"mx.size / data.size"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Tenemos que más de un 10 % de los elementos de la matriz cumplen la condición de ser «mínimos locales», así que no es nada despreciable. Esto en nuestro ejemplo hace *un total de más de 400 000 elementos*:"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"443641"
]
},
"execution_count": 9,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"mx.size"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Ahora la idea de crear dos listas y añadir los elementos uno a uno me gusta todavía menos, así que voy a cambiar de enfoque. Lo que voy a hacer va a ser crear otro array, de la misma forma que nuestros datos, y almacenar un valor `True` en aquellos elementos que cumplan la condición de mínimo local. De esta forma cumplo también una de las reglas de oro de Software Carpentry: \"*Always initialize from data*\"."
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"def busca_min_np(malla):\n",
" minimos = np.zeros_like(malla, dtype=bool)\n",
" for i in range(1, malla.shape[1]-1):\n",
" for j in range(1, malla.shape[0]-1):\n",
" if (malla[j, i] < malla[j-1, i-1] and\n",
" malla[j, i] < malla[j-1, i] and\n",
" malla[j, i] < malla[j-1, i+1] and\n",
" malla[j, i] < malla[j, i-1] and\n",
" malla[j, i] < malla[j, i+1] and\n",
" malla[j, i] < malla[j+1, i-1] and\n",
" malla[j, i] < malla[j+1, i] and\n",
" malla[j, i] < malla[j+1, i+1]):\n",
" minimos[i, j] = True\n",
"\n",
" return np.nonzero(minimos)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Encima puedo aprovechar la estupenda función `nonzero` de NumPy. Compruebo que las salidas son iguales:"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"np.testing.assert_array_equal(busca_min(data)[0], busca_min_np(data)[0])\n",
"np.testing.assert_array_equal(busca_min(data)[1], busca_min_np(data)[1])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Y evalúo los rendimientos:"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"1 loops, best of 3: 4.75 s per loop\n"
]
}
],
"source": [
"%timeit busca_min(data)"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"1 loops, best of 3: 4.62 s per loop\n"
]
}
],
"source": [
"%timeit busca_min_np(data)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Parece que los tiempos son más o menos parecidos, pero al menos ya no tengo dos objetos en memoria que van a crecer de manera aleatoria. Vamos a ver ahora cómo nos puede ayudar numba a acelerar este código."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Paso 2: aplicando `numba.jit(nopython=True)`\n",
"\n",
"Como hemos dicho antes, vamos a forzar que numba funcione en modo *nopython* para garantizar que obtenemos una mejora en el rendimiento. Si intentamos compilar la función definida en primer lugar va a fallar, porque ya hemos dicho más arriba que una de las condiciones es que *no se puede asignar memoria a objetos nuevos*:"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {
"collapsed": true
},
"outputs": [
{
"ename": "NotImplementedError",
"evalue": "Failed at nopython (nopython frontend)\n(