#!/usr/bin/env python
# coding: utf-8

# ## Volvemos al torneo del rendimiento!!!

# Recapitulando. Un artículo sobre Cython donde conseguíamos [mejoras de velocidad de código Python con numpy arrays de 40x usando Cython](http://pybonacci.org/2015/03/09/c-elemental-querido-cython/) desembocó [en mejoras de 70x usando numba](http://pybonacci.org/2015/03/13/como-acelerar-tu-codigo-python-con-numba/). En esta tercera toma vamos a ver si con Cython conseguimos las velocidades de numba tomando algunas ideas de la implementación de JuanLu y definiendo una función un poco más inteligente que mi implementación con Cython ([busca_min_cython9](http://pybonacci.org/2015/03/09/c-elemental-querido-cython/#Cythonizando,-que-es-gerundio-%28toma-9%29.)).

# Preparamos el *setup inicial*.

# In[1]:


import numpy as np
import numba

np.random.seed(0)

data = np.random.randn(2000, 2000)


# JuanLu hizo alguna trampa usando un numpy array en lugar de dos listas y devolviendo el resultado usando `numpy.nonzero`. En realidad no es trampa, es pura envidia mía al ver que ha usado una forma más inteligente de conseguir lo mismo que hacía mi función original :-P
# 
# Usando esa implementación considero que es más inteligente tener un numpy array de salida por lo que el uso de `np.nonzero` sería innecesario y añadiría algo de pérdida de rendimiento si luego vamos a seguir trabajando con numpy arrays. Por tanto, la implementación de JuanLu eliminando el uso de `numpy.nonzero` sería:

# In[2]:


def busca_min_np_jit(malla):
    minimos = np.zeros_like(malla, dtype=bool)
    _busca_min(malla, minimos)
    return minimos  # en lugar de 'return np.nonzero(minimos)'

@numba.jit(nopython=True)
def _busca_min(malla, minimos):
    for i in range(1, malla.shape[1]-1):
        for j in range(1, malla.shape[0]-1):
            if (malla[j, i] < malla[j-1, i-1] and
                malla[j, i] < malla[j-1, i] and
                malla[j, i] < malla[j-1, i+1] and
                malla[j, i] < malla[j, i-1] and
                malla[j, i] < malla[j, i+1] and
                malla[j, i] < malla[j+1, i-1] and
                malla[j, i] < malla[j+1, i] and
                malla[j, i] < malla[j+1, i+1]):
                minimos[i, j] = True


# In[3]:


get_ipython().run_line_magic('timeit', '-n 100 busca_min_np_jit(data)')


# Ejecutándolo 100 veces obtenemos un valor más bajo de 33.6 ms devolviendo un numpy.array de 1's y 0's con los 1's indicando la posición de los máximos.
# 
# La implementación original la vamos a modificar un poco para que devuelva lo mismo.

# In[4]:


def busca_min(malla):
    minimos = np.zeros_like(malla)
    for i in range(1, malla.shape[1]-1):
        for j in range(1, malla.shape[0]-1):
            if (malla[j, i] < malla[j-1, i-1] and
                malla[j, i] < malla[j-1, i] and
                malla[j, i] < malla[j-1, i+1] and
                malla[j, i] < malla[j, i-1] and
                malla[j, i] < malla[j, i+1] and
                malla[j, i] < malla[j+1, i-1] and
                malla[j, i] < malla[j+1, i] and
                malla[j, i] < malla[j+1, i+1]):
                minimos[i, j] = 1

    return minimos


# In[5]:


get_ipython().run_line_magic('timeit', 'busca_min(data)')


# Los tiempos son similares a la función original y, aunque estamos usando más memoria, tenemos una mejora con numba que ya llega a los dos órdenes de magnitud (alrededor de 100x!!) y una función más usable para trabajar con numpy.

# Vamos a modificar la opción Cython más rápida que obtuvimos para que se comporte igual que las de Numba y Python.
# 
# Primero cargamos la extensión Cython.

# In[6]:


# antes cythonmagic
get_ipython().run_line_magic('load_ext', 'Cython')


# Vamos a usar la opción `annotate` para ver cuanto 'blanco' tenemos y la nueva versión Cython la vamos a llamar `busca_min_cython10`.

# In[7]:


get_ipython().run_cell_magic('cython', '--annotate', 'import numpy as np\nfrom cython cimport boundscheck, wraparound\n\ncpdef char[:,::1] busca_min_cython10(double[:, ::1] malla):\n    cdef unsigned int i, j\n    cdef unsigned int ii = malla.shape[1]-1\n    cdef unsigned int jj = malla.shape[0]-1\n    cdef char[:,::1] minimos = np.zeros_like(malla, dtype = np.int8)\n    #minimos[...] = 0\n    cdef unsigned int start = 1\n    #cdef float [:, :] malla_view = malla\n    with boundscheck(False), wraparound(False):\n        for j in range(start, ii):\n            for i in range(start, jj):\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] = 1\n\n    return minimos\n')


# Vemos que la mayor parte está en 'blanco'. Eso significa que estamos evitando usar la C-API de CPython y la mayor parte sucede en C. Estoy usando `typed memoryviews` que permite trabajar de forma 'transparente' numpy arrays.
# 
# Vamos a ejecutar la nueva versión 100 veces, de la misma forma que hemos hecho con Numba:

# In[9]:


get_ipython().run_line_magic('timeit', '-n 100 busca_min_cython10(data)')


# Wow, virtualmente obtenemos la misma velocidad entre Numba y Cython y dos órdenes de magnitud de mejora con respecto a la versión Python.

# In[12]:


res_numba = busca_min_np_jit(data)
res_cython = busca_min_cython10(data)
res_python = busca_min(data)

np.testing.assert_array_equal(res_numba, res_cython)
np.testing.assert_array_equal(res_numba, res_python)
np.testing.assert_array_equal(res_cython, res_python)


# Parece que el resultado es el mismo en todo momento

# Probemos con arrays de menos y más tamaño.

# In[15]:


data = np.random.randn(500, 500)
get_ipython().run_line_magic('timeit', '-n 3 busca_min_np_jit(data)')
get_ipython().run_line_magic('timeit', '-n 3 busca_min_cython10(data)')
get_ipython().run_line_magic('timeit', 'busca_min(data)')


# In[14]:


data = np.random.randn(5000, 5000)
get_ipython().run_line_magic('timeit', '-n 3 busca_min_np_jit(data)')
get_ipython().run_line_magic('timeit', '-n 3 busca_min_cython10(data)')
get_ipython().run_line_magic('timeit', 'busca_min(data)')


# Parece que las distintas versiones escalan de la misma forma y el rendimiento parece, más o menos, lineal.

# ## Conclusiones de este nuevo capítulo.

# Las conclusiones que saco yo de este mano a mano que hemos llevado a cabo JuanLu (featuring Numba) y yo (featuring Cython):
# 
# * Cython: Si te restringes a cosas sencllas, es relativamente sencillo de usar. Básicamente habría que optimizar bucles y, solo en caso de que sea necesario, añadir tipos a otras variables para evitar pasar por la C-API de CPython en ciertas operaciones puesto que puede tener un coste elevado en el rendimiento. Para cosas más complejas, a pesar de que sigue siendo más placentero que C se puede complicar un poco más (pero no mucho más, una vez que lo entendido cómo usarlo).
# * Numba: Es bastante sorprendente lo que se puede llegar a conseguir con poco esfuerzo. Parece que siempre introducirá un poco de *overhead* puesto que hace muchas cosas entre bambalinas y de la otra forma (Cython) hace lo que le digamos que haga. También es verdad que muchas cosas no están soportadas, que los errores que obtenemos puede ser un poco crípticos y se hace difícil depurar el código. Pero a pesar de todo lo anterior y conociendo el historial de la gente que está detrás del proyecto Numba creo que su futuro será brillante. Por ejemplo, [Numbagg](https://github.com/shoyer/numbagg) es una librería que usa Numba y que pretende hacer lo mismo que [bottleneck](https://github.com/kwgoodman/bottleneck) (una librería muy especializada para determinadas operaciones de Numpy), que usa Cython consiguiendo [resultados comparables aunque levemente peores](https://github.com/shoyer/numbagg#benchmarks).
# 
# No sé si habrá algún capítulo más de esta serie... Lo dejo en manos de JuanLu o de cualquiera que nos quiera enviar un nuevo post relacionado.