Instala pypy 5.0 y numpypy en un virtualenv y juega con Jupyter

Hoy voy a mostrar como usar la última versión de pypy y numpypy en vuestro linux. Para instalar pypy usaremos la versión portable creada por squeaky-pl. Numpypy lo instalaremos en un entorno virtual juntamente con Jupyter para poder hacer las pruebas en un entorno más amigable que la consola de pypy.

Requerimientos

Necesitaremos tener instalada una versión reciente de virtualenv y git.

Al lío

¡Si queréis la versión TL;DR pinchad aquí! Si sois un poco más pacientes y queréis entender un poco lo que vamos a hacer seguid leyento.

Todos los comandos que vienen a continuación los tenéis que meter en un terminal. Primero creamos un directorio que se llamará pypy50 en vuestro $HOME

mkdir $HOME/pypy50

Ahora nos vamos al directorio recién creado y nos descargamos el fichero comprimido que contiene el pypy portable de 64 bits

cd $HOME/pypy50
wget https://bitbucket.org/squeaky/portable-pypy/downloads/pypy-5.0-linux_x86_64-portable.tar.bz2

Lo desempaquetamos:

tar xvfj pypy-5.0-linux_x86_64-portable.tar.bz2

Ahora creamos un directorio bin en nuestro $HOME. Si ya existe te puedes saltar este paso:

mkdir $HOME/bin

Creamos un enlace simbólico al ejecutable del pypy portable que hemos descargado que se encontrará en la carpeta bin del directorio $HOME:

ln -s $HOME/pypy50/pypy-5.0-linux_x86_64-portable/bin/pypy $HOME/bin

Cambiamos los permisos al ejecutable para darle permisos de ejecución:

chmod +x $HOME/pypy50/pypy-5.0-linux_x86_64-portable/bin/pypy

Al final de nuestro .bashrc vamos a añadir unas pocas líneas para que se añada el directorio bin de nuestro $HOME al $PATH:

echo "" >> $HOME/.bashrc
echo "# Added path to include pypy by $USER" >> $HOME/.bashrc
echo "export PATH=$PATH:$HOME/bin" >> $HOME/.bashrc
source $HOME/.bashrc

Creamos el virtualenv con pypy (en este paso necesitaréis tener virtualenv instalado). El virtualenv se creará en la carpeta bin de nuestro $HOME y se llamará pypyvenv:

virtualenv -p pypy $HOME/bin/pypyvenv

Instalamos numpypy (numpy para pypy) en el nuevo virtualenv creado (aquí necesitarás tener git instalado). Para ello usamos el pip del entorno virtual.

$HOME/bin/pypyvenv/bin/pip install git+https://bitbucket.org/pypy/numpy.git

Instalamos Jupyter haciendo algo parecido a lo anterior (aunque esta vez lo instalamos desde pypi, no confundir con pypy):

$HOME/bin/pypyvenv/bin/pip install jupyter

Y, por último, hacemos un poco de limpieza eliminando el fichero comprimido del pypy portable que hemos descargado anteriormente:

rm $HOME/pypy50/pypy*.tar.bz2

¡¡¡Listo!!!

Usando pypy

Para usar pypy (sin numpy) puedes lanzar una consola con pypy 5.0 (compatible con CPython 2.7) escribiendo en el terminal:

pypy

Usando pypy con numpy en un notebook de jupyter

Activamos el entorno virtual recien creado. Desde el terminal escribimos:

. ~/bin/pypyvenv/bin/activate

Y arrancamos jupyter:

jupyter notebook

Y después venís aquí y me contáis vuestras experiencias con pypy y numpypy o, si habéis encontrado fallos o queréis añadir mejoras, os vais a github y abrís un issue o mandáis un Pull Request y salimos ganando todos.

Ideas para mejorar el script (con vuestros pull requests)

  • Que pregunte donde instalar el pypy portable.
  • Que pregunte si queremos una carpeta bin o no.
  • Que pregunte cómo queremos llamar al entorno virtual y dónde lo queremos instalar.
  • Que pregunte si queremos instalar Jupyter y/u otras librerías.
  • ...

Saludos.

 

El producto de matrices y el nuevo operador @

Introducción.

El 13 de septiembre de 2015 fue lanzada la versión 3.5 de Python. Entre las novedades podemos encontrar la inclusión del PEP 465 que trata sobre el nuevo operador @ para la multipliación matricial y del que hablaremos en este post. Como bien sabrán los lectores de este blog, los arrays son la piedra angular de numerosísimas áreas de la programación científica y sirven para realizar operaciones de forma masiva y mucho más eficiente. Esto, sumado a la posibilidad de utilizarlos como matrices, proporciona una herramienta muy potente para llevar a cabo operaciones algebraicas. NumPy es la librería que nos permite utilizar esta maravillosa estructura de datos y según figura en el ya citado PEP, podría ser la librería fuera de la librería estándar más importada del mundo Python.

Continue reading

C elemental, querido numba

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 desembocó en mejoras de 70x usando 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).

Preparamos el setup inicial.

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:

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
%timeit -n 100 busca_min_np_jit(data)
100 loops, best of 3: 33 ms per loop
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.

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
%timeit busca_min(data)
1 loops, best of 3: 3.4 s per loop
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.

# antes cythonmagic
%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.

%%cython --annotate
import numpy as np
from cython cimport boundscheck, wraparound

cpdef char[:,::1] busca_min_cython10(double[:, ::1] malla):
    cdef unsigned int i, j
    cdef unsigned int ii = malla.shape[1]-1
    cdef unsigned int jj = malla.shape[0]-1
    cdef char[:,::1] minimos = np.zeros_like(malla, dtype = np.int8)
    #minimos[...] = 0
    cdef unsigned int start = 1
    #cdef float [:, :] malla_view = malla
    with boundscheck(False), wraparound(False):
        for j in range(start, ii):
            for i in range(start, jj):
                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







Generated by Cython 0.22

+01: import numpy as np
  __pyx_t_1 = __Pyx_Import(__pyx_n_s_numpy, 0, -1); if (unlikely(!__pyx_t_1)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 1; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
  __Pyx_GOTREF(__pyx_t_1);
  if (PyDict_SetItem(__pyx_d, __pyx_n_s_np, __pyx_t_1) < 0) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 1; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
  __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
/* … */
  __pyx_t_1 = PyDict_New(); if (unlikely(!__pyx_t_1)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 1; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
  __Pyx_GOTREF(__pyx_t_1);
  if (PyDict_SetItem(__pyx_d, __pyx_n_s_test, __pyx_t_1) < 0) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 1; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
  __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
 02: from cython cimport boundscheck, wraparound
 03: 
+04: cpdef char[:,::1] busca_min_cython10(double[:, ::1] malla):
static PyObject *__pyx_pw_46_cython_magic_de7594aedda59602146d5e749862b110_1busca_min_cython10(PyObject *__pyx_self, PyObject *__pyx_arg_malla); /*proto*/
static __Pyx_memviewslice __pyx_f_46_cython_magic_de7594aedda59602146d5e749862b110_busca_min_cython10(__Pyx_memviewslice __pyx_v_malla, CYTHON_UNUSED int __pyx_skip_dispatch) {
  unsigned int __pyx_v_i;
  unsigned int __pyx_v_j;
  unsigned int __pyx_v_ii;
  unsigned int __pyx_v_jj;
  __Pyx_memviewslice __pyx_v_minimos = { 0, 0, { 0 }, { 0 }, { 0 } };
  unsigned int __pyx_v_start;
  __Pyx_memviewslice __pyx_r = { 0, 0, { 0 }, { 0 }, { 0 } };
  __Pyx_RefNannyDeclarations
  __Pyx_RefNannySetupContext("busca_min_cython10", 0);
/* … */
  /* function exit code */
  __pyx_L1_error:;
  __Pyx_XDECREF(__pyx_t_1);
  __Pyx_XDECREF(__pyx_t_2);
  __Pyx_XDECREF(__pyx_t_3);
  __Pyx_XDECREF(__pyx_t_4);
  __Pyx_XDECREF(__pyx_t_5);
  __PYX_XDEC_MEMVIEW(&__pyx_t_6, 1);
  __pyx_r.data = NULL;
  __pyx_r.memview = NULL;
  __Pyx_AddTraceback("_cython_magic_de7594aedda59602146d5e749862b110.busca_min_cython10", __pyx_clineno, __pyx_lineno, __pyx_filename);

  goto __pyx_L2;
  __pyx_L0:;
  if (unlikely(!__pyx_r.memview)) {
    PyErr_SetString(PyExc_TypeError,"Memoryview return value is not initialized");
  }
  __pyx_L2:;
  __PYX_XDEC_MEMVIEW(&__pyx_v_minimos, 1);
  __Pyx_RefNannyFinishContext();
  return __pyx_r;
}

/* Python wrapper */
static PyObject *__pyx_pw_46_cython_magic_de7594aedda59602146d5e749862b110_1busca_min_cython10(PyObject *__pyx_self, PyObject *__pyx_arg_malla); /*proto*/
static PyObject *__pyx_pw_46_cython_magic_de7594aedda59602146d5e749862b110_1busca_min_cython10(PyObject *__pyx_self, PyObject *__pyx_arg_malla) {
  __Pyx_memviewslice __pyx_v_malla = { 0, 0, { 0 }, { 0 }, { 0 } };
  PyObject *__pyx_r = 0;
  __Pyx_RefNannyDeclarations
  __Pyx_RefNannySetupContext("busca_min_cython10 (wrapper)", 0);
  assert(__pyx_arg_malla); {
    __pyx_v_malla = __Pyx_PyObject_to_MemoryviewSlice_d_dc_double(__pyx_arg_malla); if (unlikely(!__pyx_v_malla.memview)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 4; __pyx_clineno = __LINE__; goto __pyx_L3_error;}
  }
  goto __pyx_L4_argument_unpacking_done;
  __pyx_L3_error:;
  __Pyx_AddTraceback("_cython_magic_de7594aedda59602146d5e749862b110.busca_min_cython10", __pyx_clineno, __pyx_lineno, __pyx_filename);
  __Pyx_RefNannyFinishContext();
  return NULL;
  __pyx_L4_argument_unpacking_done:;
  __pyx_r = __pyx_pf_46_cython_magic_de7594aedda59602146d5e749862b110_busca_min_cython10(__pyx_self, __pyx_v_malla);
  int __pyx_lineno = 0;
  const char *__pyx_filename = NULL;
  int __pyx_clineno = 0;

  /* function exit code */
  __Pyx_RefNannyFinishContext();
  return __pyx_r;
}

static PyObject *__pyx_pf_46_cython_magic_de7594aedda59602146d5e749862b110_busca_min_cython10(CYTHON_UNUSED PyObject *__pyx_self, __Pyx_memviewslice __pyx_v_malla) {
  PyObject *__pyx_r = NULL;
  __Pyx_RefNannyDeclarations
  __Pyx_RefNannySetupContext("busca_min_cython10", 0);
  __Pyx_XDECREF(__pyx_r);
  if (unlikely(!__pyx_v_malla.memview)) { __Pyx_RaiseUnboundLocalError("malla"); {__pyx_filename = __pyx_f[0]; __pyx_lineno = 4; __pyx_clineno = __LINE__; goto __pyx_L1_error;} }
  __pyx_t_1 = __pyx_f_46_cython_magic_de7594aedda59602146d5e749862b110_busca_min_cython10(__pyx_v_malla, 0); if (unlikely(!__pyx_t_1.memview)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 4; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
  __pyx_t_2 = __pyx_memoryview_fromslice(__pyx_t_1, 2, (PyObject *(*)(char *)) __pyx_memview_get_char, (int (*)(char *, PyObject *)) __pyx_memview_set_char, 0);; if (unlikely(!__pyx_t_2)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 4; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
  __Pyx_GOTREF(__pyx_t_2);
  __PYX_XDEC_MEMVIEW(&__pyx_t_1, 1);
  __pyx_r = __pyx_t_2;
  __pyx_t_2 = 0;
  goto __pyx_L0;

  /* function exit code */
  __pyx_L1_error:;
  __PYX_XDEC_MEMVIEW(&__pyx_t_1, 1);
  __Pyx_XDECREF(__pyx_t_2);
  __Pyx_AddTraceback("_cython_magic_de7594aedda59602146d5e749862b110.busca_min_cython10", __pyx_clineno, __pyx_lineno, __pyx_filename);
  __pyx_r = NULL;
  __pyx_L0:;
  __PYX_XDEC_MEMVIEW(&__pyx_v_malla, 1);
  __Pyx_XGIVEREF(__pyx_r);
  __Pyx_RefNannyFinishContext();
  return __pyx_r;
}
 05:     cdef unsigned int i, j
+06:     cdef unsigned int ii = malla.shape[1]-1
  __pyx_v_ii = ((__pyx_v_malla.shape[1]) - 1);
+07:     cdef unsigned int jj = malla.shape[0]-1
  __pyx_v_jj = ((__pyx_v_malla.shape[0]) - 1);
+08:     cdef char[:,::1] minimos = np.zeros_like(malla, dtype = np.int8)
  __pyx_t_1 = __Pyx_GetModuleGlobalName(__pyx_n_s_np); if (unlikely(!__pyx_t_1)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 8; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
  __Pyx_GOTREF(__pyx_t_1);
  __pyx_t_2 = __Pyx_PyObject_GetAttrStr(__pyx_t_1, __pyx_n_s_zeros_like); if (unlikely(!__pyx_t_2)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 8; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
  __Pyx_GOTREF(__pyx_t_2);
  __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
  __pyx_t_1 = __pyx_memoryview_fromslice(__pyx_v_malla, 2, (PyObject *(*)(char *)) __pyx_memview_get_double, (int (*)(char *, PyObject *)) __pyx_memview_set_double, 0);; if (unlikely(!__pyx_t_1)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 8; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
  __Pyx_GOTREF(__pyx_t_1);
  __pyx_t_3 = PyTuple_New(1); if (unlikely(!__pyx_t_3)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 8; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
  __Pyx_GOTREF(__pyx_t_3);
  PyTuple_SET_ITEM(__pyx_t_3, 0, __pyx_t_1);
  __Pyx_GIVEREF(__pyx_t_1);
  __pyx_t_1 = 0;
  __pyx_t_1 = PyDict_New(); if (unlikely(!__pyx_t_1)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 8; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
  __Pyx_GOTREF(__pyx_t_1);
  __pyx_t_4 = __Pyx_GetModuleGlobalName(__pyx_n_s_np); if (unlikely(!__pyx_t_4)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 8; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
  __Pyx_GOTREF(__pyx_t_4);
  __pyx_t_5 = __Pyx_PyObject_GetAttrStr(__pyx_t_4, __pyx_n_s_int8); if (unlikely(!__pyx_t_5)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 8; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
  __Pyx_GOTREF(__pyx_t_5);
  __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;
  if (PyDict_SetItem(__pyx_t_1, __pyx_n_s_dtype, __pyx_t_5) < 0) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 8; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
  __Pyx_DECREF(__pyx_t_5); __pyx_t_5 = 0;
  __pyx_t_5 = __Pyx_PyObject_Call(__pyx_t_2, __pyx_t_3, __pyx_t_1); if (unlikely(!__pyx_t_5)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 8; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
  __Pyx_GOTREF(__pyx_t_5);
  __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
  __Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
  __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
  __pyx_t_6 = __Pyx_PyObject_to_MemoryviewSlice_d_dc_char(__pyx_t_5);
  if (unlikely(!__pyx_t_6.memview)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 8; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
  __Pyx_DECREF(__pyx_t_5); __pyx_t_5 = 0;
  __pyx_v_minimos = __pyx_t_6;
  __pyx_t_6.memview = NULL;
  __pyx_t_6.data = NULL;
 09:     #minimos[...] = 0
+10:     cdef unsigned int start = 1
  __pyx_v_start = 1;
 11:     #cdef float [:, :] malla_view = malla
 12:     with boundscheck(False), wraparound(False):
+13:         for j in range(start, ii):
  __pyx_t_7 = __pyx_v_ii;
  for (__pyx_t_8 = __pyx_v_start; __pyx_t_8 < __pyx_t_7; __pyx_t_8+=1) {
    __pyx_v_j = __pyx_t_8;
+14:             for i in range(start, jj):
    __pyx_t_9 = __pyx_v_jj;
    for (__pyx_t_10 = __pyx_v_start; __pyx_t_10 < __pyx_t_9; __pyx_t_10+=1) {
      __pyx_v_i = __pyx_t_10;
+15:                 if (malla[j, i] < malla[j-1, i-1] and
      __pyx_t_12 = __pyx_v_j;
      __pyx_t_13 = __pyx_v_i;
      __pyx_t_14 = (__pyx_v_j - 1);
      __pyx_t_15 = (__pyx_v_i - 1);
      __pyx_t_16 = (((*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_malla.data + __pyx_t_12 * __pyx_v_malla.strides[0]) )) + __pyx_t_13)) ))) < (*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_malla.data + __pyx_t_14 * __pyx_v_malla.strides[0]) )) + __pyx_t_15)) )))) != 0);
      if (__pyx_t_16) {
      } else {
        __pyx_t_11 = __pyx_t_16;
        goto __pyx_L8_bool_binop_done;
      }
+16:                     malla[j, i] < malla[j-1, i] and
      __pyx_t_17 = __pyx_v_j;
      __pyx_t_18 = __pyx_v_i;
      __pyx_t_19 = (__pyx_v_j - 1);
      __pyx_t_20 = __pyx_v_i;
      __pyx_t_16 = (((*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_malla.data + __pyx_t_17 * __pyx_v_malla.strides[0]) )) + __pyx_t_18)) ))) < (*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_malla.data + __pyx_t_19 * __pyx_v_malla.strides[0]) )) + __pyx_t_20)) )))) != 0);
      if (__pyx_t_16) {
      } else {
        __pyx_t_11 = __pyx_t_16;
        goto __pyx_L8_bool_binop_done;
      }
+17:                     malla[j, i] < malla[j-1, i+1] and
      __pyx_t_21 = __pyx_v_j;
      __pyx_t_22 = __pyx_v_i;
      __pyx_t_23 = (__pyx_v_j - 1);
      __pyx_t_24 = (__pyx_v_i + 1);
      __pyx_t_16 = (((*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_malla.data + __pyx_t_21 * __pyx_v_malla.strides[0]) )) + __pyx_t_22)) ))) < (*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_malla.data + __pyx_t_23 * __pyx_v_malla.strides[0]) )) + __pyx_t_24)) )))) != 0);
      if (__pyx_t_16) {
      } else {
        __pyx_t_11 = __pyx_t_16;
        goto __pyx_L8_bool_binop_done;
      }
+18:                     malla[j, i] < malla[j, i-1] and
      __pyx_t_25 = __pyx_v_j;
      __pyx_t_26 = __pyx_v_i;
      __pyx_t_27 = __pyx_v_j;
      __pyx_t_28 = (__pyx_v_i - 1);
      __pyx_t_16 = (((*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_malla.data + __pyx_t_25 * __pyx_v_malla.strides[0]) )) + __pyx_t_26)) ))) < (*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_malla.data + __pyx_t_27 * __pyx_v_malla.strides[0]) )) + __pyx_t_28)) )))) != 0);
      if (__pyx_t_16) {
      } else {
        __pyx_t_11 = __pyx_t_16;
        goto __pyx_L8_bool_binop_done;
      }
+19:                     malla[j, i] < malla[j, i+1] and
      __pyx_t_29 = __pyx_v_j;
      __pyx_t_30 = __pyx_v_i;
      __pyx_t_31 = __pyx_v_j;
      __pyx_t_32 = (__pyx_v_i + 1);
      __pyx_t_16 = (((*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_malla.data + __pyx_t_29 * __pyx_v_malla.strides[0]) )) + __pyx_t_30)) ))) < (*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_malla.data + __pyx_t_31 * __pyx_v_malla.strides[0]) )) + __pyx_t_32)) )))) != 0);
      if (__pyx_t_16) {
      } else {
        __pyx_t_11 = __pyx_t_16;
        goto __pyx_L8_bool_binop_done;
      }
+20:                     malla[j, i] < malla[j+1, i-1] and
      __pyx_t_33 = __pyx_v_j;
      __pyx_t_34 = __pyx_v_i;
      __pyx_t_35 = (__pyx_v_j + 1);
      __pyx_t_36 = (__pyx_v_i - 1);
      __pyx_t_16 = (((*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_malla.data + __pyx_t_33 * __pyx_v_malla.strides[0]) )) + __pyx_t_34)) ))) < (*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_malla.data + __pyx_t_35 * __pyx_v_malla.strides[0]) )) + __pyx_t_36)) )))) != 0);
      if (__pyx_t_16) {
      } else {
        __pyx_t_11 = __pyx_t_16;
        goto __pyx_L8_bool_binop_done;
      }
+21:                     malla[j, i] < malla[j+1, i] and
      __pyx_t_37 = __pyx_v_j;
      __pyx_t_38 = __pyx_v_i;
      __pyx_t_39 = (__pyx_v_j + 1);
      __pyx_t_40 = __pyx_v_i;
      __pyx_t_16 = (((*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_malla.data + __pyx_t_37 * __pyx_v_malla.strides[0]) )) + __pyx_t_38)) ))) < (*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_malla.data + __pyx_t_39 * __pyx_v_malla.strides[0]) )) + __pyx_t_40)) )))) != 0);
      if (__pyx_t_16) {
      } else {
        __pyx_t_11 = __pyx_t_16;
        goto __pyx_L8_bool_binop_done;
      }
+22:                     malla[j, i] < malla[j+1, i+1]):
      __pyx_t_41 = __pyx_v_j;
      __pyx_t_42 = __pyx_v_i;
      __pyx_t_43 = (__pyx_v_j + 1);
      __pyx_t_44 = (__pyx_v_i + 1);
      __pyx_t_16 = (((*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_malla.data + __pyx_t_41 * __pyx_v_malla.strides[0]) )) + __pyx_t_42)) ))) < (*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_malla.data + __pyx_t_43 * __pyx_v_malla.strides[0]) )) + __pyx_t_44)) )))) != 0);
      __pyx_t_11 = __pyx_t_16;
      __pyx_L8_bool_binop_done:;
      if (__pyx_t_11) {
+23:                     minimos[i,j] = 1
        __pyx_t_45 = __pyx_v_i;
        __pyx_t_46 = __pyx_v_j;
        *((char *) ( /* dim=1 */ ((char *) (((char *) ( /* dim=0 */ (__pyx_v_minimos.data + __pyx_t_45 * __pyx_v_minimos.strides[0]) )) + __pyx_t_46)) )) = 1;
        goto __pyx_L7;
      }
      __pyx_L7:;
    }
  }
 24: 
+25:     return minimos
  __PYX_INC_MEMVIEW(&__pyx_v_minimos, 0);
  __pyx_r = __pyx_v_minimos;
  goto __pyx_L0;

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' con numpy arrays.
Vamos a ejecutar la nueva versión 100 veces, de la misma forma que hemos hecho con Numba:

%timeit -n 100 busca_min_cython10(data)
100 loops, best of 3: 27.6 ms per loop
Wow, virtualmente obtenemos la misma velocidad entre Numba y Cython y dos órdenes de magnitud de mejora con respecto a la versión Python.

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.

data = np.random.randn(500, 500)
%timeit -n 3 busca_min_np_jit(data)
%timeit -n 3 busca_min_cython10(data)
%timeit busca_min(data)
3 loops, best of 3: 2.04 ms per loop
3 loops, best of 3: 1.75 ms per loop
1 loops, best of 3: 209 ms per loop
data = np.random.randn(5000, 5000)
%timeit -n 3 busca_min_np_jit(data)
%timeit -n 3 busca_min_cython10(data)
%timeit busca_min(data)
3 loops, best of 3: 216 ms per loop
3 loops, best of 3: 174 ms per loop
1 loops, best of 3: 21.6 s per loop
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 has 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 es una librería que usa Numba y que pretende hacer lo mismo que bottleneck (una librería muy especializada para determinadas operaciones de Numpy), que usa Cython consiguiendo resultados comparables aunque levemente peores.

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.

Comprobando presencia de subconjuntos dentro de conjunto de forma eficiente

Hoy estaba trabajando con fechas y tenía que encontrar subconjuntos de fechas que estaban presentes o no en otros conjuntos. Numpy nos ofrece varias funciones de ayuda para comprobar este tipo de cosas. Ahora mismo se me vienen a la cabeza numpy.in1d y numpy.setdiff1d. El tema es que necesitaba hacer muchas operaciones de este tipo y he estado buscando la forma de poder optimizar el código un poco y resulta que, dependiendo de la operación, numpy puede resultar hasta 10 veces más lento (según mis pruebas de hoy) que usar puro CPython3. Veamos un ejemplo sencillo:

1) Imaginemos que mis datos son arrays de numpy y quiero conocer los datos de un array a que no están presentes en b:

import numpy as np
a = np.arange(1000000)
b = np.arange(10, 1000010)
%timeit np.setdiff1d(a, b)

Usando IPython y en mi máquina el anterior código me da el siguiente resultado:

10 loops, best of 3: 148 ms per loop

2) Si ahora hago lo mismo usando conjuntos (sets) obtendría el siguiente resultado:

a = set(np.arange(1000000))
b = set(np.arange(10, 1000010))
%timeit a.difference(b)

De la misma forma, usando IPython y en mi máquina el anterior código me da el siguiente resultado:

10 loops, best of 3: 30.2 ms per loop

Para este caso concreto numpy es ¡¡¡hasta 5 veces más lento!!!

Esto es una micro-optimización que me ha servido a mi en un caso concreto y si aplica espero que os resulte útil pero ya sabéis que la optimización prematura es la fuente de todos los males :-)

Saludos

Actualización:

Como bien comentan tanto Jaime como Chema más abajo, la comparación es tramposa ya que, normalmente, lo que mido no es lo que se hace en un programa completo, se debería de medir el proceso completo. La pretensión de este post era baja y solo mostrar la diferencia de implementaciones entre lo que hace numpy y lo que hace CPython. Os recomiendo que leáis los comentarios para aprender más.

Gracias, Jaime y Chema.

Dibujando una rosa de frecuencias (reloaded)

Esta entrada es una actualización a la entrada Dibujando una rosa de frecuencias dónde se rehace el código para usar nuevas funcionalidades de matplotlib que simplifica el script.
Imaginaos que estáis de vacaciones en Agosto en la playa y la única preocupación que tenéis es observar las nubes. Como sois un poco frikis y no podéis desconectar de vuestra curiosidad científica decidís apuntar las ocurrencias de la procedencia de las nubes y al final de las vacaciones decidís representar esos datos. La forma más normal de hacerlo sería usando una rosa de frecuencias.
Primero de todo vamos a importar los módulos que nos harán falta:
In [2]:
import sys

import numpy as np
import matplotlib.pyplot as plt
import matplotlib
import math

%matplotlib inline

print('Versión de Python usada: ', sys.version)
print('Versión de Numpy usada: ', np.__version__)
print('Versión de Matplotlib usada: ', matplotlib.__version__)
Versión de Python usada:  3.4.1 (v3.4.1:c0e311e010fc, May 18 2014, 10:38:22) [MSC v.1600 32 bit (Intel)]
Versión de Numpy usada:  1.8.1
Versión de Matplotlib usada:  1.3.1

A continuación creamos nuestra muestra de datos totalmente inventada:
In [3]:
## Creamos un conjunto de datos
datos = np.arange(10,90,10)
## Los datos los queremos en tanto por ciento
datos = datos * 100. / datos.sum()
## Direcciones en radianes empezando por el N
## A las direcciones les restamos 22.5º para que las barras
## estén centradas exactamente en 0, 45, 90,...
direcciones = (np.arange(0, 360, 45) - 22.5) * math.pi / 180.
sectores = ['N','NE','E','SE','S','SW','W','NW']
En el bloque anterior de código, lo único que hemos hecho es crear un conjunto de datos sin sentido y los hemos separado en 8 intervalos que pretenden ser las 8 direcciones de donde provienen las nubes empezando por el Norte y en el sentido de las agujas del reloj. Finalmente los datos los expresamos como frecuencia en tanto por ciento en cada una de las 8 direcciones.
Matplotlib nos permite hacer gráficos polares pero estos gráficos están pensados para gráficos en sentido contrario a las agujas del reloj y empezando a las tres en punto (o al este). Por ello debemos modificar como se verán los datos en el gráfico polar. Para ello definimos el tipo de gráfico, colocamos el nombre de la dirección en cada sector definido (en este caso hemos usado 8 sectores), ponemos un título a nuestro gráfico y hemos acabado.
In [4]:
fig = plt.figure(figsize = (10,10))
ax = fig.add_subplot(111, polar=True)
## La siguiente línea de código hace que los datos vayan en el 
## sentido de las agujas del reloj
ax.set_theta_direction(-1)
## La siguiente línea de código coloca el 'origen' de la rotación
## donde le indiquemos, en este caso em Norte.
ax.set_theta_zero_location('N')
## Título
ax.set_title('Procedencia de las nubes en agosto (%)')
## Dibujamos los datos
ax.bar(direcciones, datos)
## Colocamos las etiquetas del eje x
ax.set_thetagrids(np.arange(0, 360, 45), sectores, frac = 1.1, fontsize = 10)
Out[4]:
(<a list of 16 Line2D ticklines objects>,
 <a list of 8 Text major ticklabel objects>)
Y listo.
Saludos.

MicroEntradas: numpy.vectorize

Hoy vamos a ver numpy.vectorize que, como la documentación oficial indica, sirve para 'vectorizar' funciones que solo aceptan escalares como entrada. La entrada que podemos meter es una lista de objetos o un 'numpyarray' y nos devolverá como resultado un 'numpyarray'.

No os engañéis, normalmente, cuando hablamos de vectorizar pensamos en varios órdenes de magnitud de mejora en el rendimiento pero eso no es lo que hace esta función :-). Podéis pensar en esta función como algo parecido a usar map.

Por ejemplo, la función abs solo acepta un entero o un float. Por ejemplo, lo siguiente:

print(abs(-3))

Nos devolvería el valor 3. Hasta ahí bien. Pero si queremos hacer lo siguiente:

print(abs([-3, -2]))

Nos devolverá un TypeError: bad operand type for abs(): 'list'.

Podríamos usar map para ello pero no obtendríamos un 'numpyarray' como salida y en Python 3 devuelve un iterador que para transformarlo a, por ejemplo, una lista o un 'numpyarray' debemos añadir un paso más, la conversión a lista (en Python 2 devuelve directamente una lista) o dos pasos, la conversión a lista y esta a 'numpyarray'.

Vamos a ver cómo podemos 'vectorizar' la función abs sin usar map y usando numpy.vectorize (antes deberéis importar numpy como np):

vectabs = np.vectorize(abs)

que podemos usar de la siguiente forma:

print(vectabs([-3, -2]))

Et voilà. ya tenemos lo que queríamos. Veamos como va el rendimiento de lo que acabamos de hacer:

In [1]: import numpy as np
In [2]: kk = np.random.randn(1000000)
In [3]: timeit np.abs(kk)
100 loops, best of 3: 4.13 ms per loop
In [4]: vectabs = np.vectorize(abs)
In [5]: timeit vectabs(kk)
10 loops, best of 3: 87.6 ms per loop
In [6]: timeit np.array([abs(i) for i in kk])
1 loops, best of 3: 241 ms per loop

La última opción es equivalente a hacer np.array(list(map(abs, kk))) en tiempo. La versión vectorizada, vectabs, es 21 veces más lenta que la que podemos obtener usando numpy.abs por lo que no tendríamos una gran ganancia pero, sin embargo, vemos que es casi tres veces más rápida que la versión usando map (o una 'list comprehension') por lo que algo ganariamos respecto a Python 3 puro :-).

Como apuntes finales, si sabéis de alguna función en CPython que no existe equivalente en numpy y la necesitáis usar quizá podéis obtener una ganancia usando numpy.vectorize. Si tenéis que escribir la función vosotros, escribidla pensando en operaciones vectorizadas usando 'numpyarrays' y no os hará falta usar numpy.vectorize.

Como punto final. recordad que np.vectorize no es más que un decorador por lo que lo siguiente sería perfectamente válido:

@np.vectorize
def mi_funcion(*args):
    ...

Saludos.

[Editado: corrección de algún bug, disculpas!!]

¿Cómo trabajar con arrays de binarios en Python?

Esta semana vamos con una pregunta de Gonzalo, que nos dice por email:

Saludos. Últimamente he estado trabajado en python y me gusta mucho. En estos momentos estoy haciendo un algoritmo genético con codificación en numeros binarios utilizando bitset de C++. ¿Tendran información acerca de como trabajar con cadenas de binarios en python? Gracias y Felicidades por el blog.

He investigado un poco acerca de bitset en la referencia de C++, y leo que se trata básicamente de un array de valores booleanos, optimizados para su almacenamiento en memoria.

Continue reading

Leer datos de Arduino desde Python

Introducción

En este artículo vamos a ver cómo leer datos procedentes de una plataforma Arduino con Python. Para quienes no lo conozcáis, Arduino es una plataforma de hardware libre concebida para crear prototipos de manera rápida y fácil usando componentes electrónicos. Gracias a Arduino vamos a alejarnos un poco de lo que solemos ver en este blog, que es solo software, y vamos a poder interactuar con el mundo real de una manera más directa.

Arduino Uno

Este artículo nace gracias a mi reciente incorporación a Aerobot, el club de robótica de mi escuela, donde iré explorando las posibilidades de Arduino y Python :)

En esta entrada se han usado python 3.3.3, numpy 1.8.0, pyserial 2.7 y matplotlib 1.3.1.

Prefacio: ¿cómo funciona?

En este artículo no vamos a ver en detalle qué es Arduino o el lenguaje con el que se programa. Para ello os remito a gente que sabe mucho más que nosotros: al final del artículo incluyo unos enlaces que os pueden interesar si queréis profundizar en este tema. Sin embargo, sí que vamos a explicar brevemente cómo es el proceso de escribir un programa para Arduino, por razones que en seguida veremos.

Tal y como se detalla en la documentación, el proceso de compilación en Arduino funciona a grandes rasgos de la siguiente manera:

  1. Se escribe un programa (sketch) en el IDE de Arduino en C o C++.
  2. El IDE comprueba que la sintaxis es correcta y añade #include "Arduino.h" y una función main() propia de la placa.
  3. El programa se compila con avr-gcc y se manda el binario a la placa.
  4. Una vez Arduino tiene el programa y mientras esté alimentado, el programa se ejecutará, normalmente de manera indefinida.

Esta explicación viene para aclarar que, al menos en este artículo, no vamos a programar nuestra placa en Python. Por el proceso que hemos visto, las únicas maneras de hacer esto serían:

  • Traducir un subconjunto de Python a C/C++, compilarlo y subirlo a la placa. Después de una búsqueda rápida en Google no tengo noticias de que nadie haya intentado esto.
  • Escribir un sketch que defina un protocolo de comunicación entre Arduino y Python. Este es el enfoque tomado por estos proyectos:
    • https://github.com/vascop/Python-Arduino-Proto-API-v2
    • https://github.com/thearn/Python-Arduino-Command-API
    • https://github.com/lekum/pyduino

Con el segundo método no tenemos disponibles todas las funciones de Arduino, de modo que solo sirve para prototipar programas. Hoy me voy a olvidar de esto y voy a usar Arduino para programar la placa y Python para obtener los datos.

Continue reading

Informando de bugs en NumPy y SciPy en castellano a través de Pybonacci

No sé si se llegará a usar esto mucho, pero he pensado que tal vez haya gente encontrando comportamientos extraños en NumPy o SciPy que pueden ser bugs pero que no puedan enviar un informe de errores en inglés a los repositorios oficiales. Por eso he hecho un fork de ambos proyectos, y si pensáis que habéis encontrado un bug en NumPy o SciPy podéis decírnoslo en castellano aquí y nosotros nos encargamos de traducirlos.

https://github.com/Pybonacci/numpy/issues
https://github.com/Pybonacci/scipy/issues

NOTA IMPORTANTE: Estos repositorios no son para consultas sobre código o dudas personales. Como para informar de cualquier otro bug se exigirá seriedad y compromiso por parte de la persona que informa. En particular, aplican estas tres sencillas reglas:

  1. El informe debe contener en qué sistema operativo se trabaja y qué versiones de NumPy y/o SciPy se manejan. Para ello se pueden ejecutar los comandos:
$ uname -a
Linux nebulae 3.10.7-1-ARCH #1 SMP PREEMPT Thu Aug 15 11:55:34 CEST 2013 x86_64 GNU/Linux
$ python -c "import numpy; print(numpy.version.version)"
1.7.1
$ python -c "import scipy; print(scipy.version.version)"
0.12.0

  1. Se debe proporcionar el mínimo fragmento de código que reproduce el problema. Si es muy breve, puede incluirse en el informe; si no, es preferible usar http://gist.github.com u otros servicios. Mandar un programa de 200 líneas donde el problema está en la 138 no ayuda a descubrir el bug.

  2. Se debe indicar cuál era el comportamiento esperado, y cuál es el obtenido. De esta forma podemos diagnosticar fácilmente dónde está el problema. Frases como «no va», «se cuelga» sin más explicaciones no son útiles tampoco.

Se ruega difusión entre las personas que pudieran estar interesadas :)

¡Saludos!

Ajuste e interpolación unidimensionales básicos en Python con SciPy

Introducción

En este artículo vamos a ver una introducción a cómo hacer ajustes e interpolaciones en Python utilizando NumPy y los módulos interpolate y optimize de SciPy.

Ajustes de curvas e interpolaciones son dos tareas básicas que realizaremos con mucha frecuencia. Por ejemplo, cuando recojamos los datos de un experimento: sabemos que se tienen que comportar como una parábola, pero obviamente por errores de medición u otro tipo no obtenemos una parábola exactamente. En este caso necesitaremos realizar un ajuste de los datos, conocido el modelo (una curva de segundo grado en este caso).

En otras ocasiones dispondremos de una serie de puntos y querremos construir una curva que pase por todos ellos. En este caso lo que queremos es realizar una interpolación: si tenemos pocos puntos podremos usar un polinomio, y en caso contrario habrá que usar trazadores (splines en inglés). Vamos a empezar por este último método.

Si deseas consultar el código completo (incluyendo el que genera las figuras) puedes ver el notebook que usé para redactar el artículo.

En esta entrada se han usado python 3.3.2, numpy 1.7.1, scipy 0.12.0 y matplotlib 1.3.0.

Interpolación

Polinomios no, ¡gracias!

Lo primero que vamos a hacer va a ser desterrar la idea de que, sea cual sea el número de puntos que tengamos, podemos construir un polinomio que pase por todos ellos «y que lo haga bien». Si tenemos \(N\) puntos nuestro polinomio tendrá que ser de grado menor o igual que \(N - 1\), pero cuando \(N\) empieza a ser grande (del orden de 10 o más) a menos que los puntos estén muy cuidadosamente elegidos el polinomio oscilará salvajemente. Esto se conoce como fenómeno de Runge.

Continue reading