Reglas para refactorizar funciones lambda

Un gran ejercicio que podéis hacer de vez en cuando es revisar la documentación oficial de Python. La misma me parece increiblemente completa aunque también un poco anárquica o sin un guión mínimamente claro para seguir diferentes tópicos.

Hoy, revisando el HOWTO de programación funcional, casi al final del documento y sin llamar la atención, he encontrado la siguiente documentación para refactorizar funciones lambda sugerida por Fredrik Lundh. Las reglas que propone para la refactorización de las funciones lambda dicen lo siguiente:

  1. Escribe una función Lambda.
  2. Escribe un comentario explicando qué se supone que hace la función lambda.
  3. Estudia el comentario durante un rato y piensa un nombre que capture la esencia del comentario.
  4. Convierte la función lambda a una declaración def usando el nombre pensado en el anterior paso.
  5. Elimina el comentario.

😛

Feliz año 2016.

Joyas Python del 2015

Este es un resumen de algunas joyas que he descubierto este 2015 dentro del mundo Python. Que las haya descubierto en el 2015 no significa que necesariamente sean cosas novedosas pero la mayoría siguen de actualidad. Tampoco es un resumen ordenado. de hecho, es un pequeño cajón de sastre. Tampoco es temático sobre ciencia, aunque la mayoría están relacionadas con ciencia ya que es a lo que me dedico. En las siguientes líneas nombro muchas cosas pero solo incluyo enlaces sobre las joyas de las que quiero hablar.

WEB:

  • En el pasado he trasteado algo con Django para hacer cosas que se puedan compartir con mucha otra gente. Django es un framework web muy completo o, como se suele decir, con baterías incluidas y el de más amplio uso dentro del mundo Python. El hecho de incluir tantas cosas de uso habitual en un desarrollo web es su fuerte para la mayoría de desarrolladores pero también su talón de Aquiles para gente que solo lo usa de vez en cuando para hacer cosas simples. Demasiado sobrecargado para acordarte de todo ello cuando lo usas muy eventualmente y demasiado condicionante para hacer cosas simples sin un guión claro. Es por ello que este año he empezado a trastear con Flask. Lo recomiendo para gente que quiere convertir una idea simple en algo usable en poco tiempo. He prestado algo de atención a wagtail y me gustaría encontrar un tutorial para gente especialmente lerda en desarrollo web (que no se lo vendan a un desarrollador Django, vamos) y que no tiene tiempo.
  • Relacionado con el trasteo anterior, he empezado a trastear también con SQLAlchemy. Como Flask no te aporta de serie su propia idea de ORM, como sí hace Django, puedes acoplar el ORM que elijas, usar SQL a capón, Mongo,... Facilita mucho el trabajo con BBDD. Y aquí puedes encontrar una serie de recursos relacionados con SQLAlchemy.
  • También relacionado con el uso de Flask, he estado trasteando con Babel para internacionalizar 'cosas' (poder hacer uso de distintos idiomas). Es increible la facilidad de uso pero he encontrado ciertos problemas que no he sabido resolver (aun no sé muy bien porqué, seguramente mi poca experiencia con la biblioteca).

GRÁFICAS:

  • ¿Quieres hacer un mapa interactivo con Python? Hasta ahora había usado mis propias soluciones. Mis soluciones son fáciles de usar y fácilmente portables a la web de forma independiente pero requieren aprender, por ejemplo, OpenLayers o Leaflet y para hacer algo simple puede resultar excesivo. Pero para otros casos de uso hay otras soluciones que pueden resultar más convenientes. Es por ello que en los últimos tiempos he estado usando Folium. Es muy simple de usar y para según que necesidad es muy apropiado. Por otra parte, quizá su diseño limite un poco las posibilidades reales. Es por ello que, después de investigar un poco, descubrí mplleaflet. Esta última librería sigue la misma filosofía que mpld3, usa matplotlib y exporta el código en algo que es capaz de interpretar la librería javascript de turno (d3js para el caso de mpld3 o leaflet para el caso de mplleaflet). Las posibilidades de uso que se me ocurren para mplleaflet son muchas.
  • Otra joyita para hacer análisis estadístico y visualización es Seaborn. Es una delicia ver como con tan poco se puede hacer tanto. Una capa sobre otra capa sobre otra capa,..., dan un gran poder con un mínimo esfuerzo. Se pierde poder de personalización pero se gana inmediatez y, en el mundo del 'loquieroparahoy', es una gran ventaja eso de la inmediatez.
  • Una pequeña tontería pero que te puede resultar útil en algún sistema donde es difícil usar un interfaz gráfico o quieres tener algo ligero para hacer gráficas de ¿baja calidad? sería bashplotlib (hasta el nombre mola).
  • He empezado a trastear algo con Plotly pero los resultados no han sido especialmente buenos (le tendré que dar una nueva oportunidad en 2016):

UTILIDADES:

  • Una pequeña tontada para el día a día sería tqdm, que te añade una barra de progreso a los bucles de tu código.

RENDIMIENTO:

  • La depuración y optimización de código siempre es algo árido y gris. La optimización prematura es la raiz de todo mal. Juntamos las churras con las merinas y nos sale que tienes que probar line_profiler sin ningún género de dudas. Date un paseo por tú código paso a paso y descubre qué es lo que está haciendo que tooooodo sea tan lento.
  • Para correr código más rápido en Python mis opciones de hoy en día serían, por orden de qué es lo que intentaría antes, numba (si es código científico) o pypy (si solo uso numpy mezclado con cosas más estándar que no dependan de bibliotecas que usan la C-API de CPython). Si Numba no funciona y pypy no resuleve la papeleta valoro si el código lo voy a necesitar ejecutar muchas veces y el tiempo que tarda es inasumible en la mayoría de ocasiones y, si es así, tiro de Cython que, con un poco de esfuerzo y afeando un poco el código Python original, permite obtener velocidades cercanas a o del mismo orden que C/C++/Fortran.

LIBRERÍA ESTÁNDAR Y ALGUNAS ALTERNATIVAS:

  • De la librería estándar he estado usando bastante argparse, collections e itertools. Los tres no tienen nada que ver, los tres son muy potentes y, sabiendo usarlos, los tres se hacen imprescindibles. Quizá para el año que viene me ponga como deberes mirar más a fondo click como mejora a argparse y functools, toolz y/o CyToolz en combinación con collections e itertools.

AÑO 2016 (DEBERES QUE ME PONGO):

  • dask
  • Más PyTables.
  • Creo que le voy a dar bastante más a d3js (por dictadura del navegador).
  • scikit-extremes, mi propia solución al análisis de extremos univariantes en Python (se aceptan ayudas).
  • PyMC y/o PyMC3.

¿Y vuestras joyitas? Animáos a compartirlas en los comentarios, independientemente que estén relacionadas con ciencia o no.

Saludos y feliz año!!

Aprende historia gracias a geocodificación inversa, mapping y wikipedia

El otro día, mientras esperaba en Juan Bravo a un amigo, tuve algo de tiempo para divagar y entre esas divagaciones junté Juan Bravo, Python, Internet, geolocalización, historia,... En fin, que estaba en Juan Bravo, no tenía ni idea de quien era ese señor (llamadme ignorante), tenía algo de tiempo y se me ocurrió poder obtener información de calles a partir de un mapa y de la wikipedia y, de aquellos polvos, estos lodos, y nació map2wiki.

¿Qué es map2wiki?

En pocas palabras, es una aplicación web que te permite buscar una calle/avenida/plaza/... en un mapa y obtener información sobre lo que le da nombre a esa dirección.

¿Por qué no buscarlo directamente en la wikipedia?

Porque eso no es tan divertido y no hubiera aprendido nada sobre Flask, Jinja2, geocodificación inversa, OpenStreetMap, Nominatim, OpenLayers, Javascript, Brython, la API de la wikipedia, mis nulas aptitudes para el diseño web (de forma expresa no quería usar cosas como bootstrap para mantenerlo lo más simple posible) aunque podría ser peor, el infierno CSS...

CSS hell

Pero vayamos por partes...

¿Qué es la geocodificación inversa?

La geocodificación inversa (reverse geocoding en inglés) permite, a partir de unas coordenadas geográficas, obtener información sobre una dirección, por ejemplo, u otras cosas que se encuentren en esas coordenadas o cerca.

OpenStreetMap ofrece una API, Nominatim, que permite hacer eso mismo, a partir de unas coordenadas geográficas se obtiene información de su base de datos.

¿Cómo funciona?

En este post voy a relatar un poco el Así se hizo sin ver código, que a veces es más interesante que la película en sí. En otro post comentaré un poco el código por si alguien quiere utilizar alguna parte para otros proyectos.

Existen varias partes que se conectan entre sí.

  • Accedes a una página servida por Flask que ofrece un mapa, gracias a openlayers + openstreetmap.
  • En el mapa, nos podemos mover hasta una dirección que debe estar en español, ya que solo está pensada para ser usada en español. En el frontend tenemos la dirección central del mapa almacenada en una variable para la latitud y otra para la longitud (gracias a Brython/JS).
  • Con la dirección definida, pulsamos sobre el botón de buscar en la wikipedia. Este botón conecta un formulario en HTML, en el cliente, con Flask, en el servidor. Tenemos un formulario con algunos campos ocultos al que le pasaremos las coordenadas obtenidas en el frontend para que sean manejadas en Python gracias a Flask.
  • Una vez que tenemos las coordenadas en Python hacemos varias cosas:
    • Primero, vamos a la API de Nominatim y le metemos las coordenadas para que nos devuelva un JSON con la información de la dirección relacionada con esas coordenadas.
    • De ese JSON extraemos la información relevante a la dirección y la 'limpiamos' (sanitizar no está en el diccionario de la RAE). En la dirección se eliminan los siguientes términos junto con los posibles preposiciones y/o artículos que pueden acompañar a esa dirección ('calle de ...', 'avenida de los ...',...)
 
["alameda", "avenida", "bulevar", "calle", "camino", 
 "carrera", "cuesta", "pasaje", "pasadizo", "paseo", "plaza", 
 "rambla", "ronda", "travesia", "via"]
  • Seguimos:
    • Con la dirección sanitizada solo nos debería quedar el nombre de la dirección ya limpio. Por ejemplo, 'Paseo de la Marquesa de Python' debería quedar como 'Marquesa de Python'. Esa dirección ya limpia se la pasamos a la API de la Wikipedia usando la librería wikipedia que podéis encontrar en pypi. Si es posible encontrar algo en la wikipedia usando el valor que le hemos pasado nos devolverá un objeto con la información relevante del artículo.
    • Con el objeto con la información de la wikipedia obtenido extraemos parte de la información y la formateamos para mostrarla en la página.
    • Una vez tenemos la información de Nominatim (el JSON con la información de la dirección) y la información devuelta por la Wikipedia tenemos todo listo para que, mediante Flask, pasar esa información a una plantilla Jinja2, que construirá el HTML final con la información del JSON obtenido y de la Wikipedia, en caso de que sea posible, o un mensaje de error, en el caso de que no sea posible.

Y este es, principalmente, todo el proceso.

En el próximo artículo nos meteremos un poco más en las tripas para poder entender mejor como se unen las piezas. Lo que veamos no pretenderá ser algo exhaustivo sobre Flask, Jinja2 u otras tecnologías.

Espero que a alguien le resulte útil:

  1. la aplicación en sí, para aprender algo de historia,
  2. la explicación del como se hizo la aplicación, para entender como se juntan las piezas del rompecabezas en una aplicación con una estructura extremadamente simple y sin base de datos detrás.

Hasta la próxima entrada.

¿Cómo funciona el método append de una lista en CPython?

Vamos a empezar con más preguntas que respuestas.

Como sabéis, las listas de Python son arrays dinámicos. Por otro lado, las tuplas son arrays estáticos.

¿Qué implica que las listas sean arrays dinámicos?

Al ser un array dinámico podemos modificar sus elementos así como extender el array (lista).

¿Cómo funciona lo de extender el array (lista)?

Cada vez que usamos el método append de las listas se crea una copia de la lista original y se añade un elemento a esa copia para luego borrar el array original.

¿Es esto último cierto?

Más o menos.

Todos estaréis conmigo que si cada vez que añadimos un nuevo elemento tenemos que crear una copia y luego eliminar el array original podríamos crear cierto coste/gasto de recursos (en memoria, principalmente, creando copias).

Veamos un poco de código:

import sys

lista = []
for i in range(100):
    lista.append(i)
    txt = 'número de elementos = {0:>3} , tamaño de la lista = {1:>4}'
    print(txt.format(i + 1, sys.getsizeof(lista)))
número de elementos =   1 , tamaño de la lista =   96
número de elementos =   2 , tamaño de la lista =   96
número de elementos =   3 , tamaño de la lista =   96
número de elementos =   4 , tamaño de la lista =   96
número de elementos =   5 , tamaño de la lista =  128
número de elementos =   6 , tamaño de la lista =  128
número de elementos =   7 , tamaño de la lista =  128
número de elementos =   8 , tamaño de la lista =  128
número de elementos =   9 , tamaño de la lista =  192
número de elementos =  10 , tamaño de la lista =  192
número de elementos =  11 , tamaño de la lista =  192
número de elementos =  12 , tamaño de la lista =  192
número de elementos =  13 , tamaño de la lista =  192
número de elementos =  14 , tamaño de la lista =  192
número de elementos =  15 , tamaño de la lista =  192
número de elementos =  16 , tamaño de la lista =  192
número de elementos =  17 , tamaño de la lista =  264
número de elementos =  18 , tamaño de la lista =  264
número de elementos =  19 , tamaño de la lista =  264
número de elementos =  20 , tamaño de la lista =  264
número de elementos =  21 , tamaño de la lista =  264
número de elementos =  22 , tamaño de la lista =  264
número de elementos =  23 , tamaño de la lista =  264
número de elementos =  24 , tamaño de la lista =  264
número de elementos =  25 , tamaño de la lista =  264
número de elementos =  26 , tamaño de la lista =  344
número de elementos =  27 , tamaño de la lista =  344
número de elementos =  28 , tamaño de la lista =  344
número de elementos =  29 , tamaño de la lista =  344
número de elementos =  30 , tamaño de la lista =  344
número de elementos =  31 , tamaño de la lista =  344
número de elementos =  32 , tamaño de la lista =  344
número de elementos =  33 , tamaño de la lista =  344
número de elementos =  34 , tamaño de la lista =  344
número de elementos =  35 , tamaño de la lista =  344
número de elementos =  36 , tamaño de la lista =  432
número de elementos =  37 , tamaño de la lista =  432
número de elementos =  38 , tamaño de la lista =  432
número de elementos =  39 , tamaño de la lista =  432
número de elementos =  40 , tamaño de la lista =  432
número de elementos =  41 , tamaño de la lista =  432
número de elementos =  42 , tamaño de la lista =  432
número de elementos =  43 , tamaño de la lista =  432
número de elementos =  44 , tamaño de la lista =  432
número de elementos =  45 , tamaño de la lista =  432
número de elementos =  46 , tamaño de la lista =  432
número de elementos =  47 , tamaño de la lista =  528
número de elementos =  48 , tamaño de la lista =  528
número de elementos =  49 , tamaño de la lista =  528
número de elementos =  50 , tamaño de la lista =  528
número de elementos =  51 , tamaño de la lista =  528
número de elementos =  52 , tamaño de la lista =  528
número de elementos =  53 , tamaño de la lista =  528
número de elementos =  54 , tamaño de la lista =  528
número de elementos =  55 , tamaño de la lista =  528
número de elementos =  56 , tamaño de la lista =  528
número de elementos =  57 , tamaño de la lista =  528
número de elementos =  58 , tamaño de la lista =  528
número de elementos =  59 , tamaño de la lista =  640
número de elementos =  60 , tamaño de la lista =  640
número de elementos =  61 , tamaño de la lista =  640
número de elementos =  62 , tamaño de la lista =  640
número de elementos =  63 , tamaño de la lista =  640
número de elementos =  64 , tamaño de la lista =  640
número de elementos =  65 , tamaño de la lista =  640
número de elementos =  66 , tamaño de la lista =  640
número de elementos =  67 , tamaño de la lista =  640
número de elementos =  68 , tamaño de la lista =  640
número de elementos =  69 , tamaño de la lista =  640
número de elementos =  70 , tamaño de la lista =  640
número de elementos =  71 , tamaño de la lista =  640
número de elementos =  72 , tamaño de la lista =  640
número de elementos =  73 , tamaño de la lista =  768
número de elementos =  74 , tamaño de la lista =  768
número de elementos =  75 , tamaño de la lista =  768
número de elementos =  76 , tamaño de la lista =  768
número de elementos =  77 , tamaño de la lista =  768
número de elementos =  78 , tamaño de la lista =  768
número de elementos =  79 , tamaño de la lista =  768
número de elementos =  80 , tamaño de la lista =  768
número de elementos =  81 , tamaño de la lista =  768
número de elementos =  82 , tamaño de la lista =  768
número de elementos =  83 , tamaño de la lista =  768
número de elementos =  84 , tamaño de la lista =  768
número de elementos =  85 , tamaño de la lista =  768
número de elementos =  86 , tamaño de la lista =  768
número de elementos =  87 , tamaño de la lista =  768
número de elementos =  88 , tamaño de la lista =  768
número de elementos =  89 , tamaño de la lista =  912
número de elementos =  90 , tamaño de la lista =  912
número de elementos =  91 , tamaño de la lista =  912
número de elementos =  92 , tamaño de la lista =  912
número de elementos =  93 , tamaño de la lista =  912
número de elementos =  94 , tamaño de la lista =  912
número de elementos =  95 , tamaño de la lista =  912
número de elementos =  96 , tamaño de la lista =  912
número de elementos =  97 , tamaño de la lista =  912
número de elementos =  98 , tamaño de la lista =  912
número de elementos =  99 , tamaño de la lista =  912
número de elementos = 100 , tamaño de la lista =  912

En el anterior código hemos creado una lista vacía y le hemos ido añadiendo elementos y hemos obtenido el tamaño de la lista usando la función getsizeof que nos indica el tamaño del objeto en bytes. Luego hemos mostrado en pantalla el número de elementos que tiene la lista y el tamaño que ocupa.

Pero, ¿qué ocurre?, ¿por qué aumentando el número de elementos, a veces, no aumenta el tamaño del objeto?, ¿por qué luego cambia?, ¿por qué a medida que hay más elementos en la lista tarda más en cambiar el tamaño de la misma?

Veamos qué dice el código original de las listas en el repo de Python localizado en Objects/listobject.c.

A partir de la línea 42 del código C podemos leer:

/* This over-allocates proportional to the list size, making room
 * for additional growth.  The over-allocation is mild, but is
 * enough to give linear-time amortized behavior over a long
 * sequence of appends() in the presence of a poorly-performing
 * system realloc().
 * The growth pattern is:  0, 4, 8, 16, 25, 35, 46, 58, 72, 88, ...
 */
new_allocated = (newsize >> 3) + (newsize < 9 ? 3 : 6);

La última línea traducida a Python sería algo así:

new_allocated = (newsize >> 3) + (3 if newsize < 9 else 6)

En el primer paréntesis tenemos el operador bitwise right shift, similar a la versión en C (no hay que olvidar que CPython está escrito en C) mientras que en el segundo paréntesis tenemos el operador ternario (sin duda, un poco más legible que la versión en C).

¿Qué está pasando aquí?

Los buenos de los core developers de CPython han pensado que si usas un array dinámico será porque quieres hacer 'perrerías' con él, como ampliarlo. Y si lo amplías una vez es probable que lo amplíes varias veces. Es por ello que, normalmente, se usa un un tamaño un poco mayor, basado en el tamaño y siquiendo la regla mostrada más arriba, para el array (lista) y, de esta forma, podemos ampliarlo sin necesidad de crear tantas copias.

Veamos esto gráficamente:

Importamos matplotlib para poder crear los gráficos.

import matplotlib.pyplot as plt
plt.style.use('ggplot')
%matplotlib inline

Creamos nuestra lista y otra lista que almacenará los tamaños en bytes, sizes.

lista = list([1])
sizes = [sys.getsizeof(lista)]

for i in range(2, 100000):
    lista.append(i)
    sizes.append(sys.getsizeof(lista))

Y ahora dibujamos los tamaños en función del número de elementos dentro de la lista:

plt.figure(figsize = (10,10))
plt.plot(lista, sizes)
plt.xlabel('Número de elementos en la lista')
plt.ylabel('Tamaño en bytes para la lista de tamaño $N$')
<matplotlib.text.Text at 0x7f0655169c88>

Vemos como los "escalones" en el tamaño de la lista con N elementos va aumentando y el escalón cada vez es más largo a medida que aumenta el tamaño de la lista.

Veamos como es el valor del tamaño dividido por el número de elementos de la lista a medida que va aumentando el mismo (los ejes de la gráfica tienen escala logarítmica):

increment = [s/l for s, l in zip(sizes, lista)]

plt.figure(figsize = (10,10))
plt.yscale('log')
plt.xscale('log')
plt.ylim(7, max(increment))
plt.plot(lista, increment)
plt.xlabel('Número de elementos en la lista')
plt.ylabel('Bytes por número de elementos')
<matplotlib.text.Text at 0x7f06552ada20>

Curioso, ¿no?

Espero que hayáis aprendido tanto como he aprendido yo elaborando esta entrada.

Vamos a comernos el mundo y [Argentina] en Python


 

Edit (2015/11/14):

Quiero agradecer a todos los que se han prestado al juego de la subasta. Al final hemos conseguido donar más de 120 € al proyecto Argentina en Python.

¡¡¡Sois grandes!!!!

Quiero agradecer, especialmente, a Yami y a Ariel su generosidad


 

Hoy nos alejamos un poco de la temática general del blog pero nos metemos de lleno en el tema de cómo crear comunidad, con amor, con dedicación y con mucho esfuerzo, en una apuesta personal de unas personas admirables.

Además del apasionante relato que tenemos hoy, si sigues leyendo, puedes llegar a conseguir una entrada para la PyConES 2015 en Valencia.

En la presente entrada vamos a charlar con Johanna, Manuel y Errante, equipo fullstack de Argentina en Python.

Pero, -¿quiénes son estos?-, os preguntáreis.

Todo esto empezó como la apuesta vital de una persona que ve la vida de una forma diferente y hermosa a la que después se le fueron sumando esfuerzos.

Pero, ¿sigo sin saber qué es todo esto?

En su propia página web podemos leer, -todo esto es un proyecto que promueve el uso del lenguaje de programación Python para resolver problemas cotidianos en usuarios comunes, como así también desarrollar programas poderosos y complejos de una forma sencilla, estimulando el aprendizaje colaborativo y la filosofía del Software Libre-.

Y, -¿cómo llevan esto adelante?-

-Para poder promover Python (entre otras cosas de interés personal) van viajando por diferentes ciudades, siguiendo la energía de las comunidades de Software Libre regionales y contactando con gente que esté interesada en un curso o una charla sobre el tema. Así como también, simplemente, compartir momentos en común charlando sobre otras actividades y conociendo diferentes culturas y formas de vida-.

Aprovechando que en el reciente viaje de Pybonacci a Argentina les pude ver en acción y quedarme, literalmente, alucinado con su buena energía ante la vida he decidido intentar aportar mi granito de arena a su proyecto promocionándolo en este humilde blog. Para ello, hemos mantenido una conversación con nuestros protagonistas a unos 11.000 km de distancia para que nos cuenten sus motivaciones, anhelos, planes, ideas, sueños,...

Empecemos:

Kiko: Primero de todo, Manuel, ¿por qué?, ¿qué te movió a dar el paso para empezar Argentina en Python?

Humitos: "La descentralización del conocimiento". ¿Qué quiero decir con esto? Me imagino que en muchos países del mundo sucede lo mismo. A mí me tocó nacer en Argentina y «en el interior». Si bien nací en la capital (Paraná) de una provincia (Entre Ríos), esta ciudad tiene muy poco movimiento en el ámbito informático, ya sea de programación como en cualquier otro.

La verdad de la Argentina es que si no vivís en una de las 3 grandes/principales ciudades del país, no existís. Estas ciudades son Buenos Aires, Córdoba y Rosario -en algún órden en particular.

Entonces, durante mis épocas de Universidad y aprendizaje feroz, me dediqué a viajar por Argentina a cuanto evento de programación/software libre había. Y claro, como te decía, todos esos eventos eran algunas de estas ciudades "principales". Así conocí mucha gente grosa de Python Argentina, creé lazos fuertes y me fui metiendo en el mundo de la informática y las comunidades.

Muchas veces intenté hacer "algo" en Paraná, donde nací, y nunca lo logré. Al menos, hablando desde la informática. Nadie utilizaba Software Libre, nadie programaba en Python, a nadie le gustaba el trabajo de sysadmin... Entonces, me sentía, y estaba, muy solo. Probé hacer algunas cosas y siempre terminé decepcionado. Incluso, ni siquiera llegaba a juntar un equipo de organizadores.

Sentía que en muchos lugares de Argentina (y el mundo) había mucha gente con mucho potencial que podía sentirse muy sola. Que tenía ganas de compartir lo que hacían con el mundo. Que querían crear una comunidad igual de buena a las de las grandes ciudades, y colaborar para que sigan creciendo con el tiempo e incluso tengan las mismas competencias que estos monstruos.

¡Salgamos a divulgar Python y a crear comunidades en todos lados!

Kiko: ¿Desde el principio pensaste en un proyecto en solitario?, ¿nadie te comentó que quisiera unirse a esta loca aventura desde la concepción de la misma?

Humitos: Bueno, en primer lugar, no es fácil sumarse (en nada) a un loco como yo. Sin embargo, tampoco pensaba en un proyecto solitario. Siempre pensé en la comunidad de Python Argentina que, aunque no estén aquí físicamente conmigo, siento que me están acompañando y apoyando en todo esto.

Pero como decía, es difícil sumarse a un loco como yo. ¿Por qué digo esto? Porque todo el mundo me decía que estaba loco y que nada de todo esto iba a funcionar. Muy pocos amigos me apoyaron desde el primer día y me dijeron que era una idea grandiosa.

Sin embargo, cuando contaba mi proyecto muchos soñaban con esto de "viajar y trabajar desde la playa" pero, por otro lado, ninguno tenía/tuvo/tiene la valentía de «dejarlo todo» y salir a comerse el mundo empezando una nueva vida. Y hoy en día, te puedo decir que dejé todo lo que tenía en Paraná, pero no fue en vano. Sino que ahora tengo mucho más; sobre todo, mucha más «vida». Sí, mi familia, mis amigos, mi ex novia y un sin fin de cosas quedaron en Paraná, pero el crecimiento personal, social y comunitario ha sido tan pero tan grande que sigue valiendo el esfuerzo.

A unas pocas semanas de emprender viaje, mientras pensaba si los que me decían que estaba loco tenían razón, me escribe un amigo de Rosario (Emiliano Dalla Verde Marcozzi) y me dice que le parecía una idea grandiosa y que sí, quería acompañarme en esta aventura. Yo, que soy un tipo bastante estructura y que tenía muchas cosas ya armadas en la cabeza, le dije: "Genial, juntémosno un fin de semana para conversar cuáles son tus expectativas y cuales las mías. Coordinemos todo y empecemos".

Luego, esta reunión se pospuso, fue pasando el tiempo y, por este motivo de "tener que dejarlo todo", se canceló. En cambio, yo como estaba dispuesto a dejarlo todo para encontrar algo mejor, seguí con mis ideas firmes y emprendí viaje...

Kiko: ¿Cuáles fueron tus expectativas iniciales cuando iniciaste el proyecto? A día de hoy, ¿estás contento con la aventura y el resultado?

Humitos: Sinceramente, nunca tuve muchas expectativas y el objetivo fue cambiando de a poquito desde el inicio hasta lo que es hoy el proyecto. En principio, cuando salí por primera vez, lo que quería era poder encontrar esa armonía de trabajo en cualquier lugar: tenía que trabajar en un bar, en una estación de gasolina, en restaurantes y hostels; no tenía una oficina fija. El trabajo me daba el dinero para poder seguir adelante (todo salía de mi bolsillo), y por eso, era muy importante.

Por eso al principio quería encontrar esa armonía laboral para luego empezar a hacer experimentos, mezclarla con los eventos, ver de hacer talleres entre semana, empezar a disfrutar los lugares en los que me encontraba y encontrar la «organización ideal». Todo esto fue cambiando con el tiempo, a tal punto que renuncié a mi trabajo por tener más tiempo para dedicarle al proyecto.

Al día de hoy, estoy muy contento con la aventura. Como decía antes, crecí mucho como persona, me relacioné con gente que, de no ser por Argentina en Python, nunca hubiese conocido, conocí muchas culturas y formas de vida diferentes, estuve en lugares maravillosos en los que muchas personas harían cualquier cosa por estar ahí, etc. Así que, por el lado de la aventura ha sido todo un éxito y no puedo pedir más 😀

Ahora, con respecto al resultado del proyecto Argentina en Python, creo que he logrado varias veces el objetivo, pero como también lo fui adaptando/modificando/incrementando todavía queda mucho por hacer. Sin embargo, mostré Python en lugares donde quizás hubiesen tenido que pasar 5~10 años para que llegue, hicimos cursos para gente que quizás le tuvo miedo a la computadora toda su vida, estuvimos programando en Django con personas mayores de edad (más de 60) y logramos que muchas mujeres se acerquen a la computación sin miedo o sin que se sientan discriminadas.

Por otro lado, ¡siempre quiero más! y siento/creo que a Argentina en Python le falta mucho. Principalmente, le falta ser más conocido. Esto nos permitiría comenzar a organizar más eventos a la distancia, sin tener que llegar al lugar para comenzar con la gestión. De esta forma, podríamos estar organizando varios eventos en paralelo y realizar más de 1 o 2 por mes en ciudades cercanas. ¡Y claro! Quiero hacer más cosas. Quiero terminar de diseñar el curso de pilas-engine que lo tengo en mente (y parte en código) para realizar un taller de videojuegos en 2d muy sencillo pero muy divertido y atrapante.

También quiero poder contagiar Argentina en Python y que se creen otros Bolivia en Python, Venezuela en Python, México en Python, Australia en Python y, ¿por qué no?, un España en Python?. ¡Quiero haya más locos de Python (como yo) por todo el mundo!

Kiko: Argentina en Python tuvo una intersección con Química & Nómada hace un tiempo. En España hay un dicho que dice 'la unión hace la fuerza' (no sé si por allí se usa esa expresión). La pregunta, ¿la unión hace la fuerza?

Johanna: Citar ese refrán me da pauta a responder con otro “¡Sin duda!, dos cabezas piensan más que una” y esto se puede ver reflejado en todo lo que hemos logrado hasta este momento y todo lo que queremos seguir logrando, pero sin dejar de resaltar que Manuel venía labrando un arduo camino desde hace un año atrás, antes de conocerlo y sumarme al proyecto.

Humitos: Primero que nada, Argentina en Python tuvo una intersección con Johanna (antes que con su proyecto) ya que ella fue de gran ayuda cuando yo estuve en Santiago del Estero, proveyéndome un alojamiento cómodo, una conexión a Internet para poder trabajar y mucha buena onda.

Luego, una vez que decidimos comenzar a viajar juntos y luego de un mes en Corrientes, Argentina, donde nos estábamos alojando en la casa de Pablo (por CouchSurfing) que vivía en frente de una plaza en un barrio muy carenciado donde había más de 15 niños (entre 5 y 14 años aprox) todos los días jugando en los árboles y demás ella dijo: "Tenemos que hacer algo con esos chicos que están ahí, sedientos de saber, de conocer y aprender"

Ahí, con esos chicos y en esa plaza, se realizó la primera Jornada de Ciencias: Jugando al científico loco que fue un éxito rotundo. Jugamos y nos divertimos por horas con esos niños. Luego de eso, ella creó su proyecto llamado Química & Nómada con la misma idea de compartir el conocimiento pero contagiando al mundo de las ciencias con esto de «compartir» que tiene el Software Libre. Así nació su proyecto un tiempo después que comenzamos a viajar juntos.

En fin, eso a modo de introducción y anécdota.

"La unión hace la fuerza", eso está más que claro en los resultados: en un año que hace que estoy viajando con ella, hicimos sólo con el proyecto de Argentina en Python más de 18 eventos en 4 países y todo el año que viajé solo, hice nada más que 4 eventos y muchísimo más chicos 😀 . Sin embargo, en todo ese año que viajé solo, aprendí a manejar el 90% de las situaciones que hemos atravesado en el año que viajé con ella. Por eso digo, si bien no pude hacer muchos eventos durante el primer tiempo, aprendí muchísimas cosas sobre el proyecto y sobre todo como hacerlo crecer.

Kiko: ¿Qué habéis ganado cada uno con la unión?

Johanna: Creo que hemos ganado demasiado ambos en esta experiencia de recorrer Latinoamérica con el objetivo principal de compartir conocimiento, el aprendizaje que trae cada día es inmenso en todos los aspectos en ocasiones inimaginables, sin embargo a nivel personal puedo contarte que mi forma de pensar, de visualizar ideas, proyecto, metas y objetivos personales han dado un giro total; si antes pensaba lograr cosas sola o con la ayuda de un director de investigación, ahora sé que puedo generar infinitos proyectos que tendrán continuidad y un acercamiento próximo a la sociedad. Me cuesta mucho salirme de los parámetros invisibles de la ciencia en solitario y ser más relajada al comunicar ideas y propuestas, pero este se convierte en un objetivo constante de cada día: aprender a comunicar con contenido y de forma sencilla.

Humitos: Como decía en la respuesta anterior, Johanna ha colaborado en todos los eventos del último año: desde hacer el Flyer, ayudar con la difusión, darme miles de opiniones sobre la organización, haciendo la registración el día del evento y poniéndose la camiseta de los eventos todos los días.

Además, participó en todos los Django Girls como Coach ayudando a otras chicas a resolver sus problemas con sus computadoras. Eso está buenísimo porque siempre tenemos, al menos, una mujer como Coach lo cual no es menor porque si hacemos eventos destinados a mujeres y no tenemos mujeres Coach "se ve raro", ¿no?

Un detalle no menor, incluso probablemente podría ser el mayor, es su indispensable ayuda en cuanto a toda la organización del viaje, matar la soledad y aumentar la socialización. ¿Qué quiero decir con esto? Cuando estaba solo, absolutamente todo pasaba por mí, decisiones, búsqueda de hospedaje, conversaciones y demás. Era muy difícil y agotador -quita mucho tiempo. Hoy en día, tenemos las tareas muy bien diferenciadas sobre quién se ocupa de cada cosa y ella está muy ligada a lo social, ya que tiene una tonada Colombiana que enamora a cualquiera que le habla y siempre conseguimos, de alguna u otra forma, lo que necesitamos 😀

Kiko: Johanna, ahora que no nos escucha Manuel y puedes hablar mal 😛 , ¿cuéntanos por qué te uniste a la aventura?

Johanna: Varias razones me llevaron a unirme en este viajar por Latinoamérica y el mundo con el proyecto Argentina en Python, pero lo que realmente afirma y concreta esta decisión de unirme finalmente al proyecto, fue la experiencia vivida en el cierre del primer día de la PyCon Ar 2014 en Rafaela, Argentina. Fue allí donde vi con mis propios ojos como trabajan en comunidad. Alguien presentaba una propuesta de algo con lo que tenia inconvenientes o quería mejorar e invitaba al resto a participar. Luego los participantes se unían al proyecto que más les gustaba y entre todos trabajaban resolviendo, mejorando y creando cosas fantásticas. Dije "¿Qué es esto?" Ten en cuenta que vengo del mundo académico de la ciencia en donde los congresos científicos tienen otra modalidad y otra forma de compartir conocimiento. Así que pensé mientras los observaba colaborar "¡Yo quiero ser parte de esta comunidad!" y aquí estoy aprendiendo a colaborar.

Kiko: Johanna, ¿nos puedes contar un poco sobre Química & Nómada?

Johanna: Química & Nómada nace de una búsqueda personal de compartir “ciencia” de una manera informal, con el objetivo principal de acercar a la comunidad al entorno científico, desde la esencia del descubrir y aprender a través de vivir la experiencia de realizar un experimento en donde se aprende a través del juego. Para llevar a cabo este compartir de ciencia, realizamos: jornadas de ciencia para chicos y cafés científicos.

El proyecto tomó forma en el mes de febrero del 2015 en el barrio obrero “Laguna Seca” de la provincia de Corrientes, Argentina. Inspirado a entrar en acción inmediata por 6 niños que jugaban en la plaza ubicada en frente de la casa que nos hospedábamos en esa estación del viaje, ahora nos encontramos en Perú, pero hemos recorrido 3 países y queremos seguir recorriendo Latinoamérica y el mundo compartiendo ciencia.

Cada experiencia que tenemos al realizar los eventos es maravillosa, en especial en las jornadas de ciencia cuando los chicos proponen modificaciones a un experimento que permite realizar ese dialogo abierto de hipótesis y análisis de resultados para comprobar el objetivo de investigación.

En conclusión, más ciencia para la comunidad, empezando por los más chicos que serán en un futuro la población pensante y para los científicos actuales una invitación a ser parte del compartir y acercar ciencia a la comunidad.

Kiko: Tengo entendido que organizásteis un teen track en la ScipyLa 2015 en Posadas (Argentina). ¿Nos podéis comentar la experiencia?

Johanna: El teen track en la ScipyLA, fue una experiencia importante y maravillosa, en el que dos aspectos esenciales se unían para que deseara participar como co-organizadora de esta jornada.

El primer aspecto es a nivel personal, era primordial el participar en la SciPy Latinoamérica, teniendo en cuenta que este evento reúne dos ámbitos de mi interés: ciencia y programación, que lo hacen ideal para mí. Teniendo en cuenta que soy una química interesada en la programación.

El segundo aspecto fue tener la posibilidad de generar un aporte al evento del cual participaría y esto era factible siendo parte de la organización de la jornada del teen track. Esta jornada además cumplía con el objetivo del proyecto Química & Nómada de acercar ciencia a la comunidad. El cual fue posible junto a Manuel Kaufmann cuando asumimos y llevamos a cabo la jornada.

En una apreciación personal sobre el resultado de la jornada Track Teen se puede definir en felicidad y motivación para  continuar gestionando estos espacios en el marco de otras conferencias sobre programación o ciencia.

Previo a eso debo agradecer a los disertantes que compartieron toda su pasión y amor por lo que hacen en su día a día en diferentes profesiones, “¡capos!”. Es la expresión que me viene a la mente en el momento de pensar en como puedo definir su maravillosa gestión dentro de la jornada. Motivamos a más de 50 adolescentes de ultimo grado de secundario de diferentes instituciones académicas de la ciudad de Posadas, Argentina, a ser los científicos o programadores del futuro.

Puedes leer más sobre el evento aquí.

Humitos: "Increíble" creo que queda chico. Cuando desde la organización de SciPyLA nos propusieron hacer esto teníamos mucho miedo y lo dudamos al principio. Lo pensamos unos días, hablamos entre nosotros y tratamos de pre-organizar lo que queríamos organizar. Pusimos mil ideas sobre la mesa, proyectos de desarrollo y un sin fin de opciones sobre qué hacer ese día.

Decidimos hacerlo y emprendernos en esa locura de organización de un evento grande, de forma remota y en el que mientras tanto íbamos a estar organizando otros eventos en otros países también. Esta parte fue muy tensa y estresante pero...

¡Los resultados fueron increíbles y valió muchísimo el esfuerzo! Tuvimos un día entero, dentro de una conferencia internacional de ciencia y tecnología, a ~30 chicos de secundaria a los que le pudimos hablar de programación, astronomía, química, electrónica, circo, comunidad y Software Libre. Además, tuvieron la oportunidad de hablar con profesionales en estos campos durante todo ese día. Logramos mezclar adolescentes con programadores profesionales del mundo real.

Personalmente, aprendí un millón de cosas y me puso en contacto con muchísima gente maravillosa: desde los docentes que fueron acompañando a los chicos hasta con Juan Luis Cano y Kiko (de @pybonacci), quienes me volaron la cabeza en diferentes universos. Incluso, Juan nos dió una mano enorme durante el Track Teen hablándole a los chicos sobre las cosas locas que pasan fuera de la tierra y captó la atención de todos los chicos, como así también la mía y se me pasaron algunas tareas que tenía asignada para hacer, Jajaja!

Hubo varios chicos que luego de la jornada se acercaron a preguntarme "¿Cómo seguir?", "¿Porqué en la escuela nos enseñan PASCAL si Python está mucho mejor?", "¿Cómo podemos hacer para participar más de este mundo?" y cosas similar. Traté de motivar a todos esos diamantes en bruto a que se sumen a todo lo que puedan. Espero haberlo logrado con algunos.

Los docentes que fueron a acompañar a los chicos, maravilloso. Todavía no puedo creer que un preceptor que fue con un curso sabía tanto de tantas cosas y participaba muy bien. Habló de programación hasta química con un conocimiento del tema que me dejó helado. A él también le dí toda la información y recursos que pude, incluso me contactó por email unos días más tarde.

Entonces, "Increíble" queda muy chico. Creo que con esa jornada logramos abrirle la cabeza a una pequeña parte de la próxima generación y eso me pone muy contento.

Kiko: Para ir terminando, Manuel, ¿por qué te llaman Humitos?

Humitos: ¡Qué buena pregunta! Nunca escribí la historia oficial en ningún lado. La he contado verbalmente, pero nunca lo registré. Así que el blog de Pybonacci va a tener la primicia 😛

De adolescente, digamos unos 13 años, cuando usábamos MS-DOS y esas computadoras (que no cualquiera podía acceder) corrían juegos como el Monkey Island, Prince of Persia, Asterix y Obelix; nosotros nos juntábamos con mi primo Gustavo en su casa en Rosario, Santa Fe y nos pasábamos muchas horas jugando y debatiendo sobre cómo avanzar al siguiente nivel -sobre todo en los de aventura gráfica.

Una de las vacaciones que viajo a Rosario para visitarlos, me encuentro que tenía un juego nuevo, también de aventura gráfica, llamado "Flight of the Amazon Queen" que en algún momento mostraba una animación donde el avión donde estaba viajando el personaje principal se estrellaba en una isla (creo que así comenzaba el juego). Entonces, uno comenzaba a caminar por toda la pantalla y en momento, cuando estaba muy lejos, se mostraba el mapa completo desde arriba en el cuál uno podía seleccionar a dónde ir para seguir con la aventura.

Ahí, en ese mapa, en el lugar dónde se había estrellado el avión, salía mucho humo y a mí eso me llamaba la atención. Pasamos horas jugando con mi primo y un amigo de él, Adrián, durante todo el día. Yo ya venía diciendo que en el "humito" estaba la solución sobre cómo seguir avanzando en el videojuego, pero como era más chico que ellos, no me prestaban mucha atención.

Luego, en un momento cuando estaban trabados y no podían seguir adelante en el juego, volví a insistir: "Andá al humito", y nada. Pasó el tiempo y de nuevo: "Andá al humito" y me dijeron: "No se puede ir al humito, nene", pero nunca habían hecho click, o sea, no sabían realmente. Así, cada vez que se preguntaban cómo resolver el acertijo yo decía: "Andá al humito, andá al humito" repetidas veces de una forma insoportable, hasta que se cansaron de mí y me echaron de la habitación por molesto.

Días más tarde, se ve que la historia se popularizó dentro del grupo de amigos de mi primo (todos más grandes que yo) ya que todos me empezaron a decir "humito" luego de esa historia y además, se reían y la contaban una y otra vez.

Finalmente, cuando me empecé a involucrar en el mundo de la informática y la tecnología vi que todos tenía un nickname y yo no tenía -si bien en Rosario me decían "humito", ese grupo era el único que me llamaba así- asi que decidí aprovecharlo para empezar a navegar la internet.

"Pero, ellos te decían «humito» y no «humitos» con 's' al final", dirás, y bueno, la historia de la 's' la voy a guardar para cuando escriba un libro 😛

Kiko: Y la última pregunta sería para Errante, ¿cómo logras soportar a estos dos durante tantas horas?

Errante: Ante todo, gracias a vos, Kiko, ya que eres el primero que me entrevista, esto me hace entregarte la primicia de toda la aventura.

Te diría que no es una tarea nada fácil, pero soy consciente del compromiso que tengo y el cambio que podemos generar con todo esto.

Por otro lado, estoy muy agradecido de estos dos locos que me han llevado por lugares extraordinarios y me hicieron conocer paisajes que de otra forma no hubiese sido posible.

Sin embargo, agradezco cuando llegamos a las metrópolis -que no me gustan para nada- ya que me dejan muchos días estacionado y no tengo que escuchar las boludeces que hablan entre ellos. Me puedo relajar, dormir tranquilo y recargar baterías para la próxima aventura...

"Mi apariencia es de un auto familiar, pero mi alma y potencia tienen todo el carácter de un 4x4."

¿Y la entrada para la PyConES 2015?

Este año adquirí mi entrada para la PyConES pero, por diversas circunstancias, no voy a poder ir. Mi entrada queda libre a disposición del que quiera participar en una ¡¡subasta!!

¿De qué se trata eso de la subasta? Como he visto que hay varias personas que se habían quedado sin entrada la idea es que la gente que quiera la entrada puje por ella empezando por 35 € (el precio original de la entrada).

¿Dónde irá el dinero de la entrada? Los 35 € que yo pagué no irán a mi bolsillo, los dono a Argentina en Python. Todo lo que pase de los 35 € lo donará el ganador de la subasta también a Argentina en Python. Por tanto, el ganador de la puja, si la misma llega, por ejemplo, a los a 50€, 35 € los pagará para la entrada y otros 15 € los donará a un proyecto tan bonito como Argentina en Python.

Para participar deja un comentario con el valor de tu puja (introduciendo un correo electrónico válido) e indicando #pujaPyConES en el comentario, por ejemplo:

#pujaPyConES

5000€

La subasta termina el día 13 a las 00.00 (es decir, podéis pujar hasta las 23.59 del día 12). Hora local de la Península Ibérica (España).

Muchas gracias a Manuel y a Johanna.

Análisis Cluster (III): Clasificación no supervisada mediante K-medias

(Este es el tercer capítulo de la mini-serie de artículos sobre análisis cluster que estamos haciendo en pybonacci, si todavía no has leído los anteriores artículos les puedes echar un ojo ahora. Escribir esta tercera parte solo ha tardado ¡¡¡tres años en salir!!!).
El algoritmo de k-medias es uno de los algoritmos más sencillos de agrupamiento dentro del campo del aprendizaje automático (machine learning para los hipsters). Sirve para agrupar $N$-elementos en $K$-grupos distintos.

Agrupamiento

(Breve resumen de lo visto en anteriores artículos).
Las técnicas de agrupamiento o análisis Cluster son un tipo de aprendizaje no supervisado ya que no requieren un aprendizaje previo a partir de los datos.
Son métodos para agrupar datos sin etiqueta en subgrupos. Estos subgrupos pretenden reflejar algún tipo de estructura interna.

K-medias

El algoritmo de K-medias es un algoritmo de agrupamiento particional que nos crea $K$ particiones de los $N$ datos (con $N \geqslant K$). El algoritmo se puede resumir en los siguientes pasos:

  • Seleccionar $K$ puntos de forma aleatoria en un dominio que comprende todos los puntos del conjunto de datos. estos $K$ puntos representan los centros iniciales de los grupos (lo veremos en el método _init_centroids de la clase KMeans que luego explicaremos más en detalle). No vamos a implementar ningún método de posicionamiento inicial complejo, los centros iniciales serán valores aleatorios colocados entre los umbrales mínimos y máximos del conjunto de puntos (entre las $x$, $y$ mínimas y máximas). El problema que resolveremos será en dos dimensiones, como ya podéis intuir.
  • Asignar cada punto del conjunto de datos al grupo con el centro más próximo (lo veremos en el método _assign_centroids de la clase KMeans, llegaremos en breve).
  • Recalculamos los centros de los grupos a partir de todos los datos asignados a cada grupo, es decir, calculamos el centroide de cada grupo método _recalculate_centroids de la clase (KMeans).
  • Repetir los dos pasos anteriores hasta que se cumpla un condición que puede ser:
    • Se ha alcanzado un número máximo de iteraciones (no deseable).
    • Los centroides ya no cambian más allá de un rango.

    Nosotros solo comprobaremos el segundo ya que solo vamos a manejar conjuntos pequeños de datos. Le pondremos un límite de variación a nuestra clase KMeans que deberán cumplir todos los centroides. La comprobación sobre si parar las iteraciones se realizarán en el método _check_stop de la clase KMeans.

Todo esto pretende ser más educativo que riguroso

La implementación que vamos a hacer será en Python puro ya que pretende ser educativa de cara a ayudar a conocer el funcionamiento del algoritmo paso a paso. Tenéis implementaciones mucho más desarrolladas, precisas, con mejor rendimiento,..., en paquetes como scipy, scikit-learn u otros.

La implementación

[AVISO PARA NAVEGANTES: Todo el código está pensado para funcionar en Python 3. Si usas Python 2 deberías empezar a pensar en actualizarte].
Primero escribo la clase a capón y luego la voy desgranando poco a poco. La he escrito como un iterador de forma que nos permita iterar fácilmente paso a paso (la iteración paso a paso la veremos de forma visual para intentar aportar aun mayor claridad). Para ver más sobre iteradores podéis DuckDuckGoear o ver este enlace.

class KMeans:
    def __init__(self, x, y, n_clusters = 1, limit = 10):
        self.x = x
        self.x_min = min(x)
        self.x_max = max(x)
        self.y = y
        self.y_min = min(y)
        self.y_max = max(y)
        self.n_clusters = n_clusters
        self.limit = limit
        self._init_centroids()        
        self.iterations = 0

    def _init_centroids(self):
        self.x_centroids = []
        self.y_centroids = []
        self.colors = []
        for i in range(self.n_clusters):
            self.x_centroids.append(randint(self.x_min, self.x_max))
            self.y_centroids.append(randint(self.y_min, self.y_max))
            r = randint(0,255)
            g = randint(0,255)
            b = randint(0,255)
            color = 'rgb({0},{1},{2})'.format(r, g, b)
            self.colors.append(color)
        self._assign_centroids()

    def _assign_centroids(self):
        self.c = []
        # Maximum possible distance to a centroid
        dmax  = sqrt((self.x_min - self.x_max)**2 + 
                     (self.y_min - self.y_max)**2)
        for xi, yi in zip(self.x, self.y):
            cluster = 0
            d0 = dmax
            for i in range(self.n_clusters):
                d = sqrt((xi - self.x_centroids[i])**2 + 
                         (yi - self.y_centroids[i])**2)
                if d < d0:
                    cluster = i
                    d0 = d
            self.c.append(cluster)

    def _recalculate_centroids(self):
        self._x_centroids = self.x_centroids[:]
        self._y_centroids = self.y_centroids[:]
        for n in range(self.n_clusters):
            x0 = 0
            y0 = 0
            cont = 0
            for i, c in enumerate(self.c):
                if c == n:
                    cont += 1
                    x0 += self.x[i]
                    y0 += self.y[i]
            self.x_centroids[n] = x0 / cont
            self.y_centroids[n] = y0 / cont
        self._assign_centroids()

    def _check_stop(self):
        for i in range(self.n_clusters):
            d = sqrt(
                (self._x_centroids[i] - self.x_centroids[i])**2 +
                (self._y_centroids[i] - self.y_centroids[i])**2
                )
            if d > self.limit:
                return False
        return True

    def __iter__(self):
        return self

    def __next__(self):
        self.iterations += 1
        self._recalculate_centroids()
        stop = self._check_stop()
        if stop == True:
            raise StopIteration
        return self

La clase anterior paso a paso:

método __init__

def __init__(self, x, y, n_clusters = 1, limit = 10):
        self.x = x
        self.x_min = min(x)
        self.x_max = max(x)
        self.y = y
        self.y_min = min(y)
        self.y_max = max(y)
        self.n_clusters = n_clusters
        self.limit = limit
        self._init_centroids()        
        self.iterations = 0

Inicializamos la clase con los valores x de los puntos, los valores y de los puntos, el número de grupos a usar, n_clusters y el límite del desplazamiento de los grupos a partir del cual consideraremos que podemos parar de iterar, limit.
Además, a partir del los valores introducidos, extraemos las ventana umbral donde se colocan los puntos (atributos x_min, x_max, y_min e y_max) que después usaremos para inicializar los centroides (mediante _init_centroids).

método _init_centroids

    def _init_centroids(self):
        self.x_centroids = []
        self.y_centroids = []
        self.colors = []
        for i in range(self.n_clusters):
            self.x_centroids.append(randint(self.x_min, self.x_max))
            self.y_centroids.append(randint(self.y_min, self.y_max))
            r = randint(0,255)
            g = randint(0,255)
            b = randint(0,255)
            color = 'rgb({0},{1},{2})'.format(r, g, b)
            self.colors.append(color)
        self._assign_centroids()

Mediante este método, que se usa al inicializar la clase, creamos los centroides iniciales a partir de valores aleatorios situados entre los valores de x_min, x_max, y_min e y_max. Además, asignamos unos colores aleatorios a cada centroide (que luego usaremos en la visualización). Pensándolo fríamente, ahora que estoy escribiendo esto, la verdad que colors podría estar fuera de esta clase pero lo vamos a dejar así.
Una vez que se han creado los centroides hacemos la asignación de cada punto del conjunto de puntos $x$ e $y$ a cada centroide mediante el método _assign_centroids.

método _assign_centroids

    def _assign_centroids(self):
        self.c = []
        # Maximum possible distance to a centroid
        dmax  = sqrt((self.x_min - self.x_max)**2 + 
                     (self.y_min - self.y_max)**2)
        for xi, yi in zip(self.x, self.y):
            cluster = 0
            d0 = dmax
            for i in range(self.n_clusters):
                d = sqrt((xi - self.x_centroids[i])**2 + 
                         (yi - self.y_centroids[i])**2)
                if d < d0:
                    cluster = i
                    d0 = d
            self.c.append(cluster)

en el atributo c (que es una lista) almacenamos el valor del grupo al que pertenece cada punto del conjunto de datos $x$ e $y$. Para ello, primero tenemos que calcular la distancia de cada punto a cada centroide y el que centroide que tenga menos distancia al punto será el asignado. Por tanto, en este paso, hemos de calcular $N \cdot K$ distancias.

método _recalculate_centroids

    def _recalculate_centroids(self):
        self._x_centroids = self.x_centroids[:]
        self._y_centroids = self.y_centroids[:]
        for n in range(self.n_clusters):
            x0 = 0
            y0 = 0
            cont = 0
            for i, c in enumerate(self.c):
                if c == n:
                    cont += 1
                    x0 += self.x[i]
                    y0 += self.y[i]
            self.x_centroids[n] = x0 / cont
            self.y_centroids[n] = y0 / cont
        self._assign_centroids()

En este paso, recalculamos los centroides. Cada nuevo centroide será el centroide de los puntos asignados a ese centroide. Los antiguos centroides los conservamos para poder compararlos con los nuevos y ver si han variado poco o mucho las nuevas posiciones de los centroides. Una vez que hemos calculado los nuevos centroides y que mantenemos los antiguos asignamos los puntos a los nuevos centroides mediante el método _assign_centroids explicado anteriormente.

método _check_stop

    def _check_stop(self):
        for i in range(self.n_clusters):
            d = sqrt(
                (self._x_centroids[i] - self.x_centroids[i])**2 +
                (self._y_centroids[i] - self.y_centroids[i])**2
                )
            if d > self.limit:
                return False
        return True

En este método calculamos si la diferencia entre la posición de los centroides antiguos y de los nuevos es superior o inferior al límite o umbral que hemos definido al instanciar la clase. Si cualquiera de los centroides se ha movido más del umbral definido seguiremos iterando (_check_stop devolverá False), si ninguno supera el umbral le diremos que pare de iterar (_check_stop devolverá True).

métodos __iter__ y __next__

    def __iter__(self):
        return self

    def __next__(self):
        self.iterations += 1
        self._recalculate_centroids()
        stop = self._check_stop()
        if stop == True:
            raise StopIteration
        return self

Si os habéis leído el enlace sobre iteradores que he dejado más arriba espero que esto sea sencillo de entender.

  • El método __iter__ no necesita mayor explicación.
  • El método __next__ incrementa el atributo iterations, que no vamos a usar pero que podríamos usar para limitar, por ejemplo, el número de iteraciones, os lo dejo como ejercicio :-D, cada vez que damos un nuevo paso recalculamos los centroides y chequeamos si hay que parar la iteración porque hemos cumplido las condiciones fijadas (sí, el stop == True es redundante pero he pretendido que sea lo más claro posible).

Visualización

El hecho de crear la clase como un iterador es porque he pensado que podríamos iterar en cada paso y hacer una visualización interactiva. Como la visualizazión quiero que funciones cuando el notebook está estático he recurrido a usar Brythonmagic.
Si queréis saber más sobre Brythonmagic podéis leer el README del enlace anterior o ver esta entrada con vídeo y todo que muestra el funcionamiento.
Como resumen, básicamente es una magic function para el notebook que nos permite usar brython (una implementación de Python 3 hecha en javascript que funciona en el navegador sin necesidad de ningún kernel detrás).
Si no lo tienes instalado lo puedes instalar mediante pip.

!pip install brythonmagic
Downloading/unpacking brythonmagic
  Downloading brythonmagic-0.1.1.tar.gz
  Running setup.py (path:/home/kiko/pyprojs/venv-scipy/build/brythonmagic/setup.py) egg_info for package brythonmagic
    
Requirement already satisfied (use --upgrade to upgrade): ipython>=1.0 in /home/kiko/pyprojs/venv-scipy/lib/python3.4/site-packages (from brythonmagic)
Installing collected packages: brythonmagic
  Running setup.py install for brythonmagic
    
Successfully installed brythonmagic
Cleaning up...
Para poder hacer uso de la extensión la hemos de cargar mediante:

%load_ext brythonmagic
El paso siguiente nos permite usar toda la maquinaria de brython en el navegador dentro del notebook y es el último paso para que todo funcione.

%%HTML
<script type="text/javascript" src="http://brython.info/src/brython_dist.js"></script>
Podemos comprobar si funciona con el siguiente ejemplo. Después de ejecutar la celda debrías de ver un mensaje en pantalla (esto seré un poco enojante cuando lo veamos en estático puesto que siempre saltará esta ventana...):

%%brython
from browser import alert
alert('it works!')
A continuación metemos el código de la clase KMeans que ya hemos detallado más arriba y la llamaremos desde otras celdas brython a partir del nombre que le hemos dado en esta celda (kmeans_class, ver primera línea de la celda)

%%brython -s kmeans_class
from math import sqrt
from random import randint

class KMeans:
    def __init__(self, x, y, n_clusters = 1, limit = 10):
        self.x = x
        self.x_min = min(x)
        self.x_max = max(x)
        self.y = y
        self.y_min = min(y)
        self.y_max = max(y)
        self.n_clusters = n_clusters
        self.limit = limit
        self._init_centroids()        
        self.iterations = 0
        
    def _init_centroids(self):
        self.x_centroids = []
        self.y_centroids = []
        self.colors = []
        for i in range(self.n_clusters):
            self.x_centroids.append(randint(self.x_min, self.x_max))
            self.y_centroids.append(randint(self.y_min, self.y_max))
            r = randint(0,255)
            g = randint(0,255)
            b = randint(0,255)
            color = 'rgb({0},{1},{2})'.format(r, g, b)
            self.colors.append(color)
        self._assign_centroids()
    
    def _assign_centroids(self):
        self.c = []
        # Maximum possible distance to a centroid
        dmax  = sqrt((self.x_min - self.x_max)**2 + 
                     (self.y_min - self.y_max)**2)
        for xi, yi in zip(self.x, self.y):
            cluster = 0
            d0 = dmax
            for i in range(self.n_clusters):
                d = sqrt((xi - self.x_centroids[i])**2 + 
                         (yi - self.y_centroids[i])**2)
                if d < d0:
                    cluster = i
                    d0 = d
            self.c.append(cluster)
    
    def _recalculate_centroids(self):
        self._x_centroids = self.x_centroids[:]
        self._y_centroids = self.y_centroids[:]
        for n in range(self.n_clusters):
            x0 = 0
            y0 = 0
            cont = 0
            for i, c in enumerate(self.c):
                if c == n:
                    cont += 1
                    x0 += self.x[i]
                    y0 += self.y[i]
            self.x_centroids[n] = x0 / cont
            self.y_centroids[n] = y0 / cont
        self._assign_centroids()
    
    def _check_stop(self):
        for i in range(self.n_clusters):
            d = sqrt(
                (self._x_centroids[i] - self.x_centroids[i])**2 +
                (self._y_centroids[i] - self.y_centroids[i])**2
                )
            if d > self.limit:
                return False
        return True
    
    def __iter__(self):
        return self
    
    def __next__(self):
        self.iterations += 1
        self._recalculate_centroids()
        stop = self._check_stop()
        if stop == True:
            raise StopIteration
        return self
El código que figura a continuación es parte de una librería que empecé hace un tiempo para hacer gráficos en el canvas del navegador de forma simple a partir de Brython. La funcionalidad básica está en el módulo base.py y lo que figura a continuación es parte de ese módulo modificado en algún sitio.
NO VOY A EXPLICAR ESTE CÓDIGO EN DETALLE PUESTO QUE NO ES PARTE DEL PROPÓSITO DE LA ENTRADA Y NO QUIERO DESVIAR LA ATENCIÓN. SI TENÉIS ALGUNA DUDA O INTERÉS PODÉIS USAR LOS COMENTARIOS DEL BLOG PARA ELLO.
El código a continuación permite dibujar sobre un canvas HTML5 diferentes formas. Solo vamos a definir formas para dibujar círculos, para los puntos y los centroides de los grupos, y líneas, para mostrar qué puntos se asignan a qué centroide en cada paso.

%%brython -s canvas_utils
from browser import document as doc
import math

## Base classes for higher level objects
class Figure:
    """
    Base class to create other elements.
    """
    def __init__(self, canvasid, 
                       facecolor = "white", 
                       edgecolor = "black", 
                       borderwidth = None):
        """        
        Parameters
        ----------
        *canvasid*: String
            String indicating the canvas id where the image should be 
            rendered.
        *facecolor*: String
            String value containing a valid HTML color
        *edgecolor*: String
            String value containing a valid HTML color
        *borderwidth*: Integer
            Value indicating the width of the border in pixels.
            If not provided it will 0 and the edgecolor will not be
            visible
        """

        if isinstance(canvasid, str):
            self.id = canvasid
        else:
            raise Exception("The canvasid parameter should be a string")
             
        try:
            self.canvas = doc[self.id]
        except:
            raise Exception("No HTML element with id=%s" %
                            self.id)
        
        try:
            self._W = self.canvas.width
            self._H = self.canvas.height
            self._ctx = self.canvas.getContext("2d")
        except:
            raise Exception("You must provide the ID of a <canvas> element")
        
        self.facecolor = facecolor
        self.borderwidth = borderwidth
        self.edgecolor = edgecolor
        self.clf()
    
    def clf(self):
        "clear the figure"
        self._ctx.save()
        
        # The following line should clear the canvas but I found a
        # problem when I use beginPath ¿¿¿???
        #self._ctx.clearRect(0, 0, self._W, self._H)
        # So I use the following line that is less performant but
        # this operation shouldn't be done very often...
        self.canvas.width = self.canvas.width
        
        self._ctx.fillStyle = self.facecolor
        self._ctx.fillRect(0, 0, self._W, self._H)
        self._ctx.fill()
        if self.borderwidth:
            self._ctx.lineWidth = self.borderwidth
            self._ctx.strokeStyle = self.edgecolor
            self._ctx.strokeRect(0, 0, self._W, self._H)
            self._ctx.stroke()
        self._ctx.restore()
        

class Shape:
    """
    Base class to create other elements.
    """
    def __init__(self, context, x, y,
                       facecolor = "black", 
                       edgecolor = "black",
                       #alpha = 1,
                       borderwidth = None):
        """        
        Parameters
        ----------
        *context*: a canvas context
            a valid canvas context where the text will be rendered
        *x*: int or float
            x value for location in pixels
        *y*: int or float
            y value for location in pixels
        *facecolor*: String
            String value containing a valid HTML color
        *edgecolor*: String
            String value containing a valid HTML color
        *alpha*: int or float
            Value between 0 (transparent) and 1 (opaque) to set the
            transparency of the text
        *borderwidth*: Integer
            Value indicating the width of the border in pixels.
            If not provided it will 0 and the edgecolor will not be
            visible
        """
        self._ctx = context
        self.x = x
        self.y = y
        self.facecolor = facecolor
        self.borderwidth = borderwidth
        self.edgecolor = edgecolor
        #self.alpha = alpha

class Circle(Shape):
    def __init__(self, *args, radius = 10, **kwargs):
        """
        Parameters
        ----------
        *radius*: int or float
            radius of the circle in pixels.
        """
        Shape.__init__(self, *args, **kwargs)
        self.r = radius
        self.draw()
    
    def draw(self):
        self._ctx.save()
        #self._ctx.globalAlpha = self.alpha
        self._ctx.beginPath()
        self._ctx.fillStyle = self.facecolor
        self._ctx.arc(self.x, self.y, self.r, 0, 2 * math.pi)
        self._ctx.fill()
        if self.borderwidth:
            self._ctx.lineWidth = self.borderwidth
            self._ctx.strokeStyle = self.edgecolor
            self._ctx.arc(self.x, self.y, self.r, 0, 2 * math.pi)
            self._ctx.stroke()
        self._ctx.closePath()
        self._ctx.restore()

class Line(Shape):
    def __init__(self, *args, polygon = False, borderwidth = 2, **kwargs):
        Shape.__init__(self, *args, **kwargs)
        self.borderwidth = borderwidth
        self.polygon = polygon
        self.draw()
    
    def draw(self):
        self._ctx.save()
        #self._ctx.globalAlpha = self.alpha
        self._ctx.beginPath()
        self._ctx.moveTo(self.x[0], self.y[0])
        for i in range(len(self.x)):
            self._ctx.lineTo(self.x[i], self.y[i])
        if self.polygon:
            self._ctx.closePath()
            if self.facecolor:
                self._ctx.fillStyle = self.facecolor
                self._ctx.fill()
        if self.borderwidth:
            self._ctx.lineWidth = self.borderwidth
            self._ctx.strokeStyle = self.edgecolor
            self._ctx.stroke()
        self._ctx.restore()
El siguiente código solo nos va a valer para establecer una mínima disposición de los elementos contenidos en la visualización.
Se incluye un canvas, donde se dibujarán las cosas, y un botón, para poder iterar el algoritmo.

html = """
<div id="main">
 <p>
  <canvas id="cnvs01" width=500 height=500></canvas>
 </p>
 <p>
  <button id="button" class="btn">Next step</button>
 </p>
</div>"""
Por último, este es el código que realiza todo a partir de los demás. Voy a intentar explicarlo un poco más en detalle:

fig = Figure('cnvs01', borderwidth = 2)

n_points = 50
x = [randint(10, fig._W - 10) for value in range(n_points)]
y = [randint(10, fig._H - 10) for value in range(n_points)]

kmeans = KMeans(x, y, n_clusters = 4, limit = 1)

def plot(obj):
    fig._ctx.save()
    fig._ctx.fillStyle= "#ffffff"
    fig._ctx.globalAlpha = 0.3
    fig._ctx.fillRect(2,2,fig._W-4,fig._H-4)
    fig._ctx.restore()
    x = obj.x
    y = obj.y
    npoints = len(x)
    colors = obj.colors
    xc = obj.x_centroids
    yc = obj.y_centroids
    c = obj.c
    for i in range(npoints):
        color = colors[c[i]]
        Line(fig._ctx, [x[i], xc[c[i]]], [y[i], yc[c[i]]],
             facecolor = color, edgecolor = color)
        Circle(fig._ctx, x[i], y[i], 
               facecolor = color, edgecolor = 'black',
               borderwidth = 1, radius = 4)
    for xci, yci, color in zip(xc, yc, colors):
        Circle(fig._ctx, xci, yci, 
               facecolor = color, edgecolor = 'black',
               borderwidth = 1, radius = 8)

def update(ev):
    plot(kmeans)
    try:
        next(kmeans)
    except:
        #doc['button'].disabled = True
        del doc['button']

doc['button'].bind('click', update)

  • fig = Figure('cnvs01', borderwidth = 2) inicializamos el canvas con un borde negro usando la clase Figure creada más arriba.
  • Definimos el número de puntos a usar y las posiciones de los puntos.
  • Inicializamos el objeto mediante kmeans = KMeans(x, y, n_clusters = 4, limit = 1). En este caso queremos que tenga 4 clusters (n_clusters) y que el límite (limit) para dejar de iterar sea cuando los centroides se muevan un pixel o menos.
  • La función plot hace varias cosas.
    1. Primero suaviza los colores de la imagen previa antes de la actual iteración (esto está más relacionado con el canvas HTML5 y con javascript y no hace falta entenderlo mucho más allá de lo comentado en esta línea)
      fig._ctx.save()
      fig._ctx.fillStyle= "#ffffff"
      fig._ctx.globalAlpha = 0.3
      fig._ctx.fillRect(2,2,fig._W-4,fig._H-4)
      fig._ctx.restore()
    2. Después extraemos todos los datos del objeto (kmeans) que vamos a usar dentro de la función.
      x = obj.x
      y = obj.y
      npoints = len(x)
      colors = obj.colors
      xc = obj.x_centroids
      yc = obj.y_centroids
      c = obj.c
    3. Dibujamos las líneas entre los puntos y los centroides y los puntos con colores asociados a cada grupo
      for i in range(npoints):
         color = colors[c[i]]
         Line(fig._ctx, [x[i], xc[c[i]]], [y[i], yc[c[i]]],
              facecolor = color, edgecolor = color)
         Circle(fig._ctx, x[i], y[i], 
                facecolor = color, edgecolor = 'black',
                borderwidth = 1, radius = 4)
    4. Y finalmente se dibujan los centroides de cada grupo con un círculo un poco más grande que los propios puntos
      for xci, yci, color in zip(xc, yc, colors):
         Circle(fig._ctx, xci, yci, 
                facecolor = color, edgecolor = 'black',
                borderwidth = 1, radius = 8)
  • Creamos, además, una función update que será la que llama a la función plot y la que llama a una nueva iteración. El código que figura en el bloque except, del doc['button'], es más de la parte de brython y sirve para que el botón desaparezca una vez que el iterador ha quedado exhausto (se han agotado las iteraciones).
  • La última línea de código, doc['button'].bind('click', update), asocia el evento click del botón a la función update que he comentado anteriormente.

Si ahora ejecutamos la siguiente celda de código deberíamos ver un canvas de 500px x 500px con un botón debajo. Si vamos pulsando al botón veremos como nueva iteración en acción, además de ver las anteriores iteraciones de forma difuminada. Una vez que hemos alcanzado la condición de que los centroides no se mueven más desaparecerá el botón (para no teneros ahí pulsando y que me digáis que eso no sigue funcionando...).

%%brython -S kmeans_class canvas_utils -h html

fig = Figure('cnvs01', borderwidth = 2)

n_points = 50

x = [randint(10, fig._W - 10) for value in range(n_points)]
y = [randint(10, fig._H - 10) for value in range(n_points)]

kmeans = KMeans(x, y, n_clusters = 4, limit = 1)

def plot(obj):
    fig._ctx.save()
    fig._ctx.fillStyle= "#ffffff"
    fig._ctx.globalAlpha = 0.3
    fig._ctx.fillRect(2,2,fig._W-4,fig._H-4)
    fig._ctx.restore()
    x = obj.x
    y = obj.y
    npoints = len(x)
    colors = obj.colors
    xc = obj.x_centroids
    yc = obj.y_centroids
    c = obj.c
    for i in range(npoints):
        color = colors[c[i]]
        Line(fig._ctx, [x[i], xc[c[i]]], [y[i], yc[c[i]]],
             facecolor = color, edgecolor = color)
        Circle(fig._ctx, x[i], y[i], 
               facecolor = color, edgecolor = 'black',
               borderwidth = 1, radius = 4)
    for xci, yci, color in zip(xc, yc, colors):
        Circle(fig._ctx, xci, yci, 
               facecolor = color, edgecolor = 'black',
               borderwidth = 1, radius = 8)

def update(ev):
    plot(kmeans)
    try:
        next(kmeans)
    except:
        #doc['button'].disabled = True
        del doc['button']
        
doc['button'].bind('click', update)
Podéis ver todo el código contenido en una página web aquí.

Conclusiones

Espero que haya quedado claro el funcionamiento básico del algoritmo K-Medias. También espero que si veis alguna errata en el código me aviséis o hagáis un Pull Request a nuestro repositorio de notebooks.
Podéis ver este notebook funcionando incluso en estático en el siguiente enlace.

Trabajando con Python y R

Trabajando de forma conjunta con Python y con R.

Hoy vamos a ver como podemos juntar lo bueno de R, algunas de sus librerías, con Python usando rpy2.

Pero, lo primero de todo, ¿qué es rpy2? rpy2 es una interfaz que permite que podamos comunicar información entre R y Python y que podamos acceder a funcionalidad de R desde Python. Por tanto, podemos estar usando Python para todo nuestro análisis y en el caso de que necesitemos alguna librería estadística especializada de R podremos acceder a la misma usando rpy2.

Para poder usar rpy2 necesitarás tener instalado tanto Python (CPython versión >= 2.7.x) como R (versión >=3), además de las librerías R a las que quieras acceder. Conda permite realizar todo el proceso de instalación de los intérpretes de Python y R, además de librerías, pero no he trabajado con Conda y R por lo que no puedo aportar mucho más en este aspecto. Supongo que será parecido a lo que hacemos con Conda y Python.

Para este microtutorial voy a hacer uso de la librería extRemes de R que permite hacer análisis de valores extremos usando varias de las metodologías más comúnmente aceptadas.

Como siempre, primero de todo, importaremos la funcionalidad que necesitamos para la ocasión.

# Importamos pandas y numpy para manejar los datos que pasaremos a R
import pandas as pd
import numpy as np

# Usamos rpy2 para interactuar con R
import rpy2.robjects as ro

# Activamos la conversión automática de tipos de rpy2
import rpy2.robjects.numpy2ri
rpy2.robjects.numpy2ri.activate()

import matplotlib.pyplot as plt
%matplotlib inline
En el anterior código podemos ver una serie de cosas nuevas que voy a explicar brevemente:

  • import rpy2.robjects as ro, esto lo explicaremos un poquito más abajo.
  • import rpy2.robjects.numpy2ri, importamos el módulo numpy2ri. Este módulo permite que hagamos conversión automática de objetos numpy a objetos rpy2.
  • rpy2.robjects.numpy2ri.activate(), hacemos uso de la función activate que activa la conversión automática de objetos que hemos comentado en la línea anterior.

Brevísima introducción a algunas de las cosas más importantes de rpy2.

Para evaluar directamente código R podemos hacerlo usando rpy2.robjects.r con el código R expresado como una cadena (rpy2.robjects lo he importado como ro en este caso, como podéis ver más arriba):

codigo_r = """
saluda <- function(cadena) {
  return(paste("Hola, ", cadena))
}
"""
ro.r(codigo_r)
<SignatureTranslatedFunction - Python:0x03096490 / R:0x03723E98>
En la anterior celda hemos creado una función R llamada saluda y que ahora está disponible en el espacio de nombres global de R. Podemos acceder a la misma desde Python de la siguiente forma:

saluda_py = ro.globalenv['saluda']
Y podemos usarla de la siguiente forma:

res = saluda_py('pepe')
print(res[0])
Hola,  pepe
En la anterior celda véis que para acceder al resultado he tenido que usar res[0]. En realidad, lo que nos devuleve rpy2 es:

print(type(res))
print(res.shape)
<class 'numpy.ndarray'>
(1,)
En este caso un numpy array con diversa información del objeto rpy2. Como el objeto solo devuelve un string pues el numpy array solo tiene un elemento.

Podemos acceder al código R de la función de la siguiente forma:

print(saluda_py.r_repr())
function (cadena) 
{
    return(paste("Hola, ", cadena))
}
Hemos visto como acceder desde Python a nombres disponibles en el entorno global de R. ¿Cómo podemos hacer para que algo que creemos en Python este accesible en R?

variable_r_creada_desde_python = ro.FloatVector(np.arange(1,5,0.1))
Veamos como es esta variable_r_creada_desde_python dentro de Python

variable_r_creada_desde_python
<FloatVector - Python:0x09D5A7B0 / R:0x07FE8900>
[1.000000, 1.100000, 1.200000, ..., 4.700000, 4.800000, 4.900000]
¿Y lo que se tendría que ver en R?

print(variable_r_creada_desde_python.r_repr())
c(1, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2, 2.1, 2.2, 
2.3, 2.4, 2.5, 2.6, 2.7, 2.8, 2.9, 3, 3.1, 3.2, 3.3, 3.4, 3.5, 
3.6, 3.7, 3.8, 3.9, 4, 4.1, 4.2, 4.3, 4.4, 4.5, 4.6, 4.7, 4.8, 
4.9)
Pero ahora mismo esa variable no está disponible desde R y no la podríamos usar dentro de código R que permanece en el espacio R (vaya lío, ¿no?)

ro.r('variable_r_creada_desde_python')
---------------------------------------------------------------------------
RRuntimeError                             Traceback (most recent call last)
<ipython-input-10-524753a78365> in <module>()
----> 1 ro.r('variable_r_creada_desde_python')

d:\users\X003621\AppData\Local\Continuum\Miniconda3\lib\site-packages\rpy2\robjects\__init__.py in __call__(self, string)
    251     def __call__(self, string):
    252         p = rinterface.parse(string)
--> 253         res = self.eval(p)
    254         return res
    255 

d:\users\X003621\AppData\Local\Continuum\Miniconda3\lib\site-packages\rpy2\robjects\functions.py in __call__(self, *args, **kwargs)
    168                 v = kwargs.pop(k)
    169                 kwargs[r_k] = v
--> 170         return super(SignatureTranslatedFunction, self).__call__(*args, **kwargs)
    171 
    172 pattern_link = re.compile(r'\\link\{(.+?)\}')

d:\users\X003621\AppData\Local\Continuum\Miniconda3\lib\site-packages\rpy2\robjects\functions.py in __call__(self, *args, **kwargs)
     98         for k, v in kwargs.items():
     99             new_kwargs[k] = conversion.py2ri(v)
--> 100         res = super(Function, self).__call__(*new_args, **new_kwargs)
    101         res = conversion.ri2ro(res)
    102         return res

RRuntimeError: Error in eval(expr, envir, enclos) : 
  object 'variable_r_creada_desde_python' not found
Vale, tendremos que hacer que sea accesible desde R de la siguiente forma:

ro.globalenv["variable_ahora_en_r"] = variable_r_creada_desde_python
print(ro.r("variable_ahora_en_r"))
[ 1.   1.1  1.2  1.3  1.4  1.5  1.6  1.7  1.8  1.9  2.   2.1  2.2  2.3  2.4
  2.5  2.6  2.7  2.8  2.9  3.   3.1  3.2  3.3  3.4  3.5  3.6  3.7  3.8  3.9
  4.   4.1  4.2  4.3  4.4  4.5  4.6  4.7  4.8  4.9]
Ahora que ya la tenemos accesible la podemos usar desde R. Por ejemplo, vamos a usar la función sum en R que suma los elementos pero directamente desde R:

print(ro.r('sum(variable_ahora_en_r)'))
print(np.sum(variable_r_creada_desde_python))
[ 118.]
118.0
Perfecto, ya sabemos, de forma muy sencilla y básica, como podemos usar R desde Python, como podemos pasar información desde R hacia Python y desde Python hacia R. ¡¡¡Esto es muy poderoso!!!, estamos juntando lo mejor de dos mundos, la solidez de las herramientas científicas de Python con la funcionalidad especializada que nos pueden aportar algunas librerías de R no disponibles en otros ámbitos.

Trabajando de forma híbrida entre Python y R

Vamos a empezar importando la librería extRemes de R:

# Importamos la librería extRemes de R
from rpy2.robjects.packages import importr
extremes = importr('extRemes')
En la anterior celda hemos hecho lo siguiente:

  • from rpy2.robjects.packages import importr, La función importr nos servirá para importar las librerías R
  • extremes = importr('extRemes'), de esta forma importamos la librería extRemes de R, sería equivalente a hacer en R library(extRemes).
Leemos datos con pandas. En el mismo repo donde está este notebook está también un fichero de texto con datos que creé a priori. Supuestamente son datos horarios de velocidad del viento por lo que vamos a hacer análisis de valores extremos de velocidad del viento horaria.

data = pd.read_csv('datasets/Synthetic_data.txt',
                   sep = '\s*', skiprows = 1, parse_dates = [[0, 1]],
                   names = ['date','time','wspd'], index_col = 0)
data.head(3)
wspd
date_time
1983-01-01 00:00:00 7.9
1983-01-01 01:00:00 8.2
1983-01-01 02:00:00 8.5
Extraemos los máximos anuales los cuales usaremos posteriormente dentro de R para hacer cálculo de valores extremos usando la distribución generalizada de valores extremos (GEV):

max_y = data.wspd.groupby(pd.TimeGrouper(freq = 'A')).max()
Dibujamos los valores máximos anuales usando Pandas:

max_y.plot(kind = 'bar', figsize = (12, 4))
<matplotlib.axes._subplots.AxesSubplot at 0x10a923d0>
Referenciamos la funcionalidad fevd (fit extreme value distribution) dentro del paquete extremes de R para poder usarla directamente con los valores máximos que hemos obtenido usando Pandas y desde Python.

fevd = extremes.fevd
Como hemos comentado anteriormente, vamos a calcular los parámetros de la GEV usando el método de ajuste GMLE (Generalised Maximum Lihelihood Estimation) y los vamos a guardar directamente en una variable Python.

Veamos la ayuda antes:

print(fevd.__doc__)
Python representation of an R function.
description
-----------


 Fit a univariate extreme value distribution functions (e.g., GEV, GP, PP, Gumbel, or Exponential) to data; possibly with covariates in the parameters.
 


fevd(
    x,
    data,
    threshold = rinterface.NULL,
    threshold_fun = ~,
    location_fun = ~,
    scale_fun = ~,
    shape_fun = ~,
    use_phi = False,
    type = c,
    method = c,
    initial = rinterface.NULL,
    span,
    units = rinterface.NULL,
    time_units = days,
    period_basis = year,
    na_action = <rpy2.rinterface.SexpVector - Python:0x116D6C50 / R:0x0C4BB100>,
    optim_args = rinterface.NULL,
    priorFun = rinterface.NULL,
    priorParams = rinterface.NULL,
    proposalFun = rinterface.NULL,
    proposalParams = rinterface.NULL,
    iter = 9999.0,
    weights = 1.0,
    blocks = rinterface.NULL,
    verbose = False,
)

x :  `fevd`: `x` can be a numeric vector, the name of a column of `data` or a formula giving the data to which the EVD is to be fit.  In the case of the latter two, the `data` argument must be specified, and must have appropriately named columns.`plot` and `print` method functions: any list object returned by `fevd`. ,

object :  A list object of class \dQuote{fevd} as returned by `fevd`. ,

data :  A data frame object with named columns giving the data to be fit, as well as any data necessary for modeling non-stationarity through the threshold and/or any of the parameters. ,

threshold :  numeric (single or vector).  If fitting a peak over threshold (POT) model (i.e., `type` = \dQuote{PP}, \dQuote{GP}, \dQuote{Exponential}) this is the threshold over which (non-inclusive) data (or excesses) are used to estimate the parameters of the distribution function.  If the length is greater than 1, then the length must be equal to either the length of `x` (or number of rows of `data`) or to the number of unique arguments in `threshold.fun`. ,

threshold.fun :  formula describing a model for the thresholds using columns from `data`.  Any valid formula will work.  `data` must be supplied if this argument is anything other than ~ 1.  Not for use with `method` \dQuote{Lmoments}. ,

location.fun :  formula describing a model for each parameter using columns from `data`.  `data` must be supplied if any of these arguments are anything other than ~ 1. ,

scale.fun :  formula describing a model for each parameter using columns from `data`.  `data` must be supplied if any of these arguments are anything other than ~ 1. ,

shape.fun :  formula describing a model for each parameter using columns from `data`.  `data` must be supplied if any of these arguments are anything other than ~ 1. ,

use.phi :  logical; should the log of the scale parameter be used in the numerical optimization (for `method` \dQuote{MLE}, \dQuote{GMLE} and \dQuote{Bayesian} only)?  For the ML and GML estimation, this may make things more stable for some data. ,

type :  `fevd`: character stating which EVD to fit.  Default is to fit the generalized extreme value (GEV) distribution function (df).`plot` method function: character describing which plot(s) is (are) desired.  Default is \dQuote{primary}, which makes a 2 by 2 panel of plots including the QQ plot of the data quantiles against the fitted model quantiles (`type` \dQuote{qq}), a QQ plot (\dQuote{qq2}) of quantiles from model-simulated data against the data, a density plot of the data along with the model fitted density (`type` \dQuote{density}) and a return level plot (`type` \dQuote{rl}). In the case of a stationary (fixed) model, the return level plot will show return levels calculated for return periods given by `return.period`, along with associated CIs (calculated using default `method` arguments depending on the estimation method used in the fit.  For non-stationary models, the data are plotted as a line along with associated effective return levels for return periods of 2, 20 and 100 years (unless `return.period` is specified by the user to other values.  Other possible values for `type` include \dQuote{hist}, which is similar to \dQuote{density}, but shows the histogram for the data and \dQuote{trace}, which is not used for L-moment fits.  In the case of MLE/GMLE, the trace yields a panel of plots that show the negative log-likelihood and gradient negative log-likelihood (note that the MLE gradient is currently used even for GMLE) for each of the estimated parameter(s); allowing one parameter to vary according to `prange`, while the others remain fixed at their estimated values.  In the case of Bayesian estimation, the \dQuote{trace} option creates a panel of plots showing the posterior df and MCMC trace for each parameter. ,

method :  `fevd`: character naming which type of estimation method to use.  Default is to use maximum likelihood estimation (MLE). ,

initial :  A list object with any named parameter component giving the initial value estimates for starting the numerical optimization (MLE/GMLE) or the MCMC iterations (Bayesian).  In the case of MLE/GMLE, it is best to obtain a good intial guess, and in the Bayesian case, it is perhaps better to choose poor initial estimates.  If NULL (default), then L-moments estimates and estimates based on Gumbel moments will be calculated, and whichever yields the lowest negative log-likelihood is used.  In the case of `type` \dQuote{PP}, an additional MLE/GMLE estimate is made for the generalized Pareto (GP) df, and parameters are converted to those of the Poisson Process (PP) model.  Again, the initial estimates yielding the lowest negative log-likelihoo value are used for the initial guess. ,

span :  single numeric giving the number of years (or other desired temporal unit) in the data set.  Only used for POT models, and only important in the estimation for the PP model, but important for subsequent estimates of return levels for any POT model.  If missing, it will be calculated using information from `time.units`. ,

units :  (optional) character giving the units of the data, which if given may be used subsequently (e.g., on plot axis labels, etc.). ,

time.units :  character string that must be one of \dQuote{hours}, \dQuote{minutes}, \dQuote{seconds}, \dQuote{days}, \dQuote{months}, \dQuote{years}, \dQuote{m/hour}, \dQuote{m/minute}, \dQuote{m/second}, \dQuote{m/day}, \dQuote{m/month}, or \dQuote{m/year}; where m is a number.  If `span` is missing, then this argument is used in determining the value of `span`.  It is also returned with the output and used subsequently for plot labelling, etc. ,

period.basis :  character string giving the units for the period.  Used only for plot labelling and naming output vectors from some of the method functions (e.g., for establishing what the period represents for the return period). ,

rperiods :  numeric vector giving the return period(s) for which it is desired to calculate the corresponding return levels. ,

period :  character string naming the units for the return period. ,

burn.in :  The first `burn.in` values are thrown out before calculating anything from the MCMC sample. ,

a :  when plotting empirical probabilies and such, the function `ppoints` is called, which has this argument `a`. ,

d :  numeric determining how to scale the rate parameter for the point process.  If NULL, the function will attempt to scale based on the values of `period.basis` and `time.units`, the first of which must be \dQuote{year} and the second of which must be one of \dQuote{days}, \dQuote{months}, \dQuote{years}, \dQuote{hours}, \dQuote{minutes} or \dQuote{seconds}.  If none of these are the case, then `d` should be specified, otherwise, it is not necessary. ,

density.args :  named list object containing arguments to the `density` and `hist` functions, respectively. ,

hist.args :  named list object containing arguments to the `density` and `hist` functions, respectively. ,

na.action :  function to be called to handle missing values.  Generally, this should remain at the default (na.fail), and the user should take care to impute missing values in an appropriate manner as it may have serious consequences on the results. ,

optim.args :  A list with named components matching exactly any arguments that the user wishes to specify to `optim`, which is used only for MLE and GMLE methods.  By default, the \dQuote{BFGS} method is used along with `grlevd` for the gradient argument.  Generally, the `grlevd` function is used for the `gr` option unless the user specifies otherwise, or the optimization method does not take gradient information. ,

priorFun :  character naming a prior df to use for methods GMLE and Bayesian.  The default for GMLE (not including Gumbel or Exponential types) is to use the one suggested by Martins and Stedinger (2000, 2001) on the shape parameter; a beta df on -0.5 to 0.5 with parameters `p` and `q`.  Must take `x` as its first argument for `method` \dQuote{GMLE}.  Optional arguments for the default function are `p` and `q` (see details section).The default for Bayesian estimation is to use normal distribution functions.  For Bayesian estimation, this function must take `theta` as its first argument.Note: if this argument is not NULL and `method` is set to \dQuote{MLE}, it will be changed to \dQuote{GMLE}. ,

priorParams :  named list containing any prior df parameters (where the list names are the same as the function argument names).  Default for GMLE (assuming the default function is used) is to use `q` = 6 and `p` = 9.  Note that in the Martins and Stedinger (2000, 2001) papers, they use a different EVD parametrization than is used here such that a positive shape parameter gives the upper bounded distribution instead of the heavy-tail one (as emloyed here).  To be consistent with these papers, `p` and `q` are reversed inside the code so that they have the same interpretation as in the papers.Default for Bayesian estimation is to use ML estimates for the means of each parameter (may be changed using `m`, which must be a vector of same length as the number of parameters to be estimated (i.e., if using the default prior df)) and a standard deviation of 10 for all other parameters (again, if using the default prior df, may be changed using `v`, which must be a vector of length equal to the number of parameters). ,

proposalFun :  For Bayesian estimation only, this is a character naming a function used to generate proposal parameters at each iteration of the MCMC.  If NULL (default), a random walk chain is used whereby if theta.i is the current value of the parameter, the proposed new parameter theta.star is given by theta.i + z, where z is drawn at random from a normal df. ,

proposalParams :  A named list object describing any optional arguments to the `proposalFun` function.  All functions must take argument `p`, which must be a vector of the parameters, and `ind`, which is used to identify which parameter is to be proposed.  The default `proposalFun` function takes additional arguments `mean` and `sd`, which must be vectors of length equal to the number of parameters in the model (default is to use zero for the mean of z for every parameter and 0.1 for its standard deviation). ,

iter :  Used only for Bayesian estimation, this is the number of MCMC iterations to do. ,

weights :  numeric of length 1 or n giving weights to be applied     in the likelihood calculations (e.g., if there are data points to     be weighted more/less heavily than others). ,

blocks :  An optional list containing information required to fit point process models in a computationally-efficient manner by using only the exceedances and not the observations below the threshold(s). See details for further information.       ,

FUN :  character string naming a function to use to estimate the parameters from the MCMC sample.  The function is applied to each column of the `results` component of the returned `fevd` object. ,

verbose :  logical; should progress information be printed to the screen?  If TRUE, for MLE/GMLE, the argument `trace` will be set to 6 in the call to `optim`. ,

prange :  matrix whose columns are numeric vectors of length two for each parameter in the model giving the parameter range over which trace plots should be made.  Default is to use either +/- 2 * std. err. of the parameter (first choice) or, if the standard error cannot be calculated, then +/- 2 * log2(abs(parameter)).  Typically, these values seem to work very well for these plots. ,

... :  Not used by most functions here.  Optional arguments to `plot` for the various `plot` method functions.In the case of the `summary` method functions, the logical argument `silent` may be passed to suppress (if TRUE) printing any information to the screen. ,

Y ahora vamos a hacer un cálculo sin meternos mucho en todas las opciones posibles.

res = fevd(max_y.values, type = "GEV", method = "GMLE")
¿Qué estructura tiene la variable res que acabamos de crear y que tiene los resultados del ajuste?

print(type(res))
<class 'rpy2.robjects.vectors.ListVector'>
print(res.r_repr)
<bound method ListVector.r_repr of <ListVector - Python:0x10AB8878 / R:0x0CA9B458>
[Vector, ndarray, ndarray, ..., ndarray, ListV..., ListV...]
  call: <class 'rpy2.robjects.vectors.Vector'>
  <Vector - Python:0x10AB8418 / R:0x0CB2FFB4>
[RNULLType, Vector, Vector, Vector]
  data.name: <class 'numpy.ndarray'>
  array(['structure(c(22.2, 25.5, 21.5, 22.5, 23.7, 22.5, 21.7, 29.7, 24.2, ',
       '23.8, 28.1, 23.4, 23.7, 25.6, 23.2, 24.9, 22.8, 24.6, 22.3, 25.5, ',
       '22.6, 24, 20.8, 23.5, 24.4, 24.1, 25.1, 19.4, 22.8, 24.2, 25, ',
       '25.3), .Dim = 32L)', ''], 
      dtype='<U66')
  weights: <class 'numpy.ndarray'>
  array([ 1.])
  ...
  call: <class 'numpy.ndarray'>
  array(['location', 'scale', 'shape'], 
      dtype='<U8')
<ListVector - Python:0x10AB8878 / R:0x0CA9B458>
[Vector, ndarray, ndarray, ..., ndarray, ListV..., ListV...]
<ListVector - Python:0x10AB8878 / R:0x0CA9B458>
[Vector, ndarray, ndarray, ..., ndarray, ListV..., ListV...]>
Según nos indica lo anterior, ahora res es un vector que está compuesto de diferentes elementos. Los vectores pueden tener un nombre para todos o algunos de los elementos. Para acceder a estor nombres podemos hacer:

res.names
array(['call', 'data.name', 'weights', 'in.data', 'x', 'priorFun',
       'priorParams', 'method', 'type', 'period.basis', 'par.models',
       'const.loc', 'const.scale', 'const.shape', 'n', 'na.action',
       'parnames', 'results', 'initial.results'], 
      dtype='<U15')
Según el output anterior, parece que hay un nombre results, ahí es donde se guardan los valores del ajuste, los estimadores. Para acceder al mismo podemos hacerlo de diferentes formas. Con Python tendriamos que saber el índice y acceder de forma normal (__getitem__()). Existe una forma alternativa usando el método rx que nos permite acceder directamente con el nombre:

results = res.rx('results')
print(results.r_repr)
<bound method ListVector.r_repr of <ListVector - Python:0x10ABBCB0 / R:0x0CBFBC40>
[ListVector]
<ListVector - Python:0x10ABBCB0 / R:0x0CBFBC40>
[ListVector]>
Parece que tenemos un único elemento:

results = results[0]
results.r_repr
<bound method ListVector.r_repr of <ListVector - Python:0x10ABF490 / R:0x0C851BA0>
[ndarray, ndarray, ndarray, ..., RNULL..., ndarray, ListV...]
  par: <class 'numpy.ndarray'>
  array([ 23.06394152,   1.75769129,  -0.16288164])
  value: <class 'numpy.ndarray'>
  array([  1.00000000e+16])
  counts: <class 'numpy.ndarray'>
  array([1, 1], dtype=int32)
  ...
  par: <class 'rpy2.rinterface.RNULLType'>
  rpy2.rinterface.NULL
  value: <class 'numpy.ndarray'>
  array([[ 0.,  0.,  0.],
       [ 0.,  0.,  0.],
       [ 0.,  0.,  0.]])
<ListVector - Python:0x10ABF490 / R:0x0C851BA0>
[ndarray, ndarray, ndarray, ..., RNULL..., ndarray, ListV...]>
Vemos ahora que results tiene un elemento con nombre par donde se guardan los valores de los estimadores del ajuste a la GEV que hemos obtenido usando GMLE. Vamos a obtener finalmente los valores de los estimadores:

location, scale, shape = results.rx('par')[0][:]
print(location, scale, shape)
23.0639415199 1.75769128743 -0.162881636772

Funcion mágica para R (antigua rmagic)

Usamos la antigua función mágica rmagic que ahora se activará en el notebook de la siguiente forma:

%load_ext rpy2.ipython
Veamos como funciona la functión mágica de R:

help(rpy2.ipython.rmagic.RMagics.R)
Help on function R in module rpy2.ipython.rmagic:

R(self, line, cell=None, local_ns=None)
    ::
    
      %R [-i INPUT] [-o OUTPUT] [-n] [-w WIDTH] [-h HEIGHT] [-p POINTSIZE]
             [-b BG] [--noisolation] [-u {px,in,cm,mm}] [-r RES]
             
] Execute code in R, optionally returning results to the Python runtime. In line mode, this will evaluate an expression and convert the returned value to a Python object. The return value is determined by rpy2's behaviour of returning the result of evaluating the final expression. Multiple R expressions can be executed by joining them with semicolons:: In [9]: %R X=c(1,4,5,7); sd(X); mean(X) Out[9]: array([ 4.25]) In cell mode, this will run a block of R code. The resulting value is printed if it would printed be when evaluating the same code within a standard R REPL. Nothing is returned to python by default in cell mode:: In [10]: %%R ....: Y = c(2,4,3,9) ....: summary(lm(Y~X)) Call: lm(formula = Y ~ X) Residuals: 1 2 3 4 0.88 -0.24 -2.28 1.64 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 0.0800 2.3000 0.035 0.975 X 1.0400 0.4822 2.157 0.164 Residual standard error: 2.088 on 2 degrees of freedom Multiple R-squared: 0.6993,Adjusted R-squared: 0.549 F-statistic: 4.651 on 1 and 2 DF, p-value: 0.1638 In the notebook, plots are published as the output of the cell:: %R plot(X, Y) will create a scatter plot of X bs Y. If cell is not None and line has some R code, it is prepended to the R code in cell. Objects can be passed back and forth between rpy2 and python via the -i -o flags in line:: In [14]: Z = np.array([1,4,5,10]) In [15]: %R -i Z mean(Z) Out[15]: array([ 5.]) In [16]: %R -o W W=Z*mean(Z) Out[16]: array([ 5., 20., 25., 50.]) In [17]: W Out[17]: array([ 5., 20., 25., 50.]) The return value is determined by these rules: * If the cell is not None (i.e., has contents), the magic returns None. * If the final line results in a NULL value when evaluated by rpy2, then None is returned. * No attempt is made to convert the final value to a structured array. Use %Rget to push a structured array. * If the -n flag is present, there is no return value. * A trailing ';' will also result in no return value as the last value in the line is an empty string. optional arguments: -i INPUT, --input INPUT Names of input variable from shell.user_ns to be assigned to R variables of the same names after calling self.pyconverter. Multiple names can be passed separated only by commas with no whitespace. -o OUTPUT, --output OUTPUT Names of variables to be pushed from rpy2 to shell.user_ns after executing cell body (rpy2's internal facilities will apply ri2ro as appropriate). Multiple names can be passed separated only by commas with no whitespace. -n, --noreturn Force the magic to not return anything. Plot: Arguments to plotting device -w WIDTH, --width WIDTH Width of plotting device in R. -h HEIGHT, --height HEIGHT Height of plotting device in R. -p POINTSIZE, --pointsize POINTSIZE Pointsize of plotting device in R. -b BG, --bg BG Background of plotting device in R. SVG: SVG specific arguments --noisolation Disable SVG isolation in the Notebook. By default, SVGs are isolated to avoid namespace collisions between figures.Disabling SVG isolation allows to reference previous figures or share CSS rules across a set of SVGs. PNG: PNG specific arguments -u <{px,in,cm,mm}>, --units <{px,in,cm,mm}> Units of png plotting device sent as an argument to *png* in R. One of ["px", "in", "cm", "mm"]. -r RES, --res RES Resolution of png plotting device sent as an argument to *png* in R. Defaults to 72 if *units* is one of ["in", "cm", "mm"]. code
A veces, será más simple usar la función mágica para interactuar con R. Veamos un ejemplo donde le pasamos a R el valor obtenido de la función fevd del paquete extRemes de R que he usado anteriormente y corremos cierto código directamente desde R sin tener que usar ro.r.

%R -i res plot.fevd(res)
En la anterior celda de código le he pasado como parámetro de entrada (- i res) la variable res que había obtenido anteriormente para que esté disponible desde R. y he ejecutado código R puro (plot.fevd(res)).

Si lo anterior lo quiero hacer con rpy2 puedo hacer lo siquiente:

CUIDADO, la siguiente celda de código puede provocar que se reinicialice el notebook y se rompa la sesión. Si has hecho cambios en el notebook guárdalos antes de ejecutar la celda, por lo que pueda pasar...

ro.globalenv['res'] = res
ro.r("plot.fevd(res)")
rpy2.rinterface.NULL
Lo anterior me bloquea el notebook y me 'rompe' la sesión (en windows, al menos) ya que la ventana de gráficos se abre de forma externa... Por tanto, una buena opción para trabajar de forma interactiva con Python y R de forma conjunta y que no se 'rompa' nada es usar tanto rpy2 como su extensión para el notebook de Jupyter (dejaremos de llamarlo IPython poco a poco).

Usando Python y R combinando rpy2 y la función mágica

Vamos a combinar las dos formas de trabajar con rpy2 en el siguiente ejemplo:

metodos = ["MLE", "GMLE"]
tipos = ["GEV", "Gumbel"]
Lo que vamos a hacer es calcular los parámetros del ajuste usando la distribución GEV y Gumbel, que es un caso especial de la GEV. El ajuste lo calculamos usando tanto MLE como GMLE. Además de mostrar los valores resultantes del ajuste para los estimadores vamos a mostrar el dibujo de cada uno de los ajustes y algunos test de bondad. Usamos Python para toda la maquinaria de los bucles, usamos rpy2 para obtener los estimadores y usamos la función mágica de rpy2 para mostrar los gráficos del resultado.

for t in tipos:
    for m in metodos:
        print('tipo de ajuste: ', t)
        print('método de ajuste: ', m)
        res = fevd(max_y.values, method = m, type = t)
        if m == "Bayesian":
            print(res.rx('results')[0][-1][0:-2])
        elif m == "Lmoments":
            print(res.rx('results')[0])
        else:
            print(res.rx('results')[0].rx('par')[0][:])
        %R -i res plot.fevd(res)
tipo de ajuste:  GEV
método de ajuste:  MLE
[ 23.05170779   1.80858528  -0.14979836]
tipo de ajuste:  GEV
método de ajuste:  GMLE
[ 23.06394152   1.75769129  -0.16288164]
tipo de ajuste:  Gumbel
método de ajuste:  MLE
[ 22.90587606   1.81445179]
tipo de ajuste:  Gumbel
método de ajuste:  GMLE
[ 22.90587606   1.81445179]

Comentarios finales

Espero que este microtutorial os valga, al menos, para conocer rpy2 y la potencia que os puede llegar a aportar a vuestros análisis 'pythónicos'. Como resumen:

  • Tenemos en nuestras manos una herramienta muy poderosa.
  • Rpy2 puede estar poco madura en algún aspecto aunque ha mejorado bastante con respecto a alguna versión de rpy2 que usé anteriormente.
  • También podéis usar directamente R como kernel, aunque perdéis la interacción con Python. También puede ocurrir que la instalación os haga perder mucho el tiempo para poder hacerlo funcionar si os veis obligados a usarlo desde windows.
  • En la elaboración de este microtutorial la consola de R donde iba haciendo algunas pruebas simples se me ha 'roto' muchísimas más veces de las que consideraría aceptables. No se puede quedar colgada, cerrar,..., seis o siete veces en media hora una consola haciendo cosas simples. Eso hace que si quieres usar R de forma interactiva debas usar alternativas como Jupyter, RStudio u otros que desconozco ya que la consola oficial no está 'ni pa pipas' (por lo menos en Windows, el sistema operativo con más usuarios potenciales, mal que me pese).
  • Sigo manteniendo muchas reservas respecto a R como Lenguaje de Programación (en mayúsculas) por lo que si puedo limitar su uso a alguna librería especializada que necesito y a la que pueda acceder con rpy2 es lo que seguiré haciendo (If you are using R and you think you're in hell, this is a map for you.)

Y el notebook...

En el caso de que queráis trastear con el notebook lo podéis descargar desde aquí.

También podéis descargar todos los notebooks desde nuestro repo oficial de notebooks.

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 😛
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.

tutormagic: Jupyter + pythontutor

Esta será una microentrada para presentar una extensión para el notebook que estoy usando en un curso interno que estoy dando en mi empresa.
Si a alguno más os puede valer para mostrar cosas básicas de Python (2 y 3, además de Java y Javascript) para muy principiantes me alegro.

Nombre en clave: tutormagic

Esta extensión lo único que hace es embeber dentro de un IFrame la página de pythontutor usando el código que hayamos definido en una celda de código precedida de la cell magic %%tutor.
Como he comentado anteriormente, se puede escribir código Python2, Python3, Java y Javascript, que son los lenguajes soportados por pythontutor.

Ejemplo

Primero deberemos instalar la extensión. Está disponible en pypi por lo que la podéis instalar usando pip install tutormagic. Una vez instalada, dentro de un notebook de IPython la deberías cargar usando:

%load_ext tutormagic
Una vez hecho esto ya deberiamos tener disponible la cell magic para ser usada. Podéis ver un ejemplo en este notebook.

Y eso es todo

Lo dicho, espero que sea útil para alguien.

Saludos.