Como mejorar tu script fácilmente

Esta entrada apareció originalmente en inglés en mi blog.

Nos ha pasado a todos. Ese momento en el que descubres que sabes suficiente sobre un lenguage de programacion que quieres ponerlo en práctica y construir "algo", lo que sea.
Una de las mejores cosas de la comunidad de Python es no sólo su habilidad para construir cosas increíbles, sino también para compartirlas con todo el mundo, mejorando la comunidad en el proceso.

Sin embargo, llevo un tiempo fijándome en un patrón que se repite en algunos de estos proyectos. Seguro que has visto alguno así. Hablo de esos proyectos con 2 ó 3 componentes, donde el README tiene una pequeña descripción del proyecto, quizás un par de lineas explicando como ejecutar el proyecto, y frases del tipo, "Seguramente añadiré X o Y si tengo tiempo".

El caso es que muchos de estos proyectos son realmente interesantes, y tienen algún tipo de componentes que me gustaría usar sin tener que implementarlos yo mismo.

Te voy a mostrar 3 formas distintas de implementar uno de estos proyectos, cada una de ellas mejor (desde mi punto de vista) que la anterior:

Supongamos que queremos construir un script genial, donde la funcionalidad principal será que, dado un número entero por el usuario, realizará un calculo simple en base a ese entero, y devolverá el resultado.

Implementación 1

 
#!/usr/bin/env python

"""
Super awesome script
Asks the user for a number:
 - If the number is less or equal to 100, it returns the 1st tetration of the number (power of itself)
 - else, it returns the number squared
"""

__version__ = '0.1'

if __name__ == '__main__':

    while 1:
        user_number = input('Choose a number:\n') #raw_input() in python2
        if user_number.isdigit():
            user_number = int(user_number)
            break
        else:
            print('{} is not a valid number'.format(user_number))

    if user_number > 100:
        print(user_number**2)
    else:
        print(user_number**user_number)

Ésta suele ser la implementación de alquien que lleva poco tiempo en python. Funciona, pregunta al usuario por el input, realiza la operación, e imprime en pantalla el resultado.

Veo dos problemas en esta implementación:

1. No hay ningún tipo de separación entre la lógica de la interacción del usuario y la lógica del cálculo. Todo esta incluido en el mismo macro bloque. Pese a ser funcional, esta implementación hace que sea díficil el modificar o expandir este script (para hacerlo tendrías que leerte todo el código).

2. Estamos gestionando toda la validación por nuestra cuenta. Python tiene formas de hacer esto para que tú no te tengas que molestar en hacerlo :).

Para la siguiente implementación, usaremos el módulo mas simple de la libreria standard para trabajar con inputs del usuario, .

Implementación 2

 
#!/usr/bin/env python

"""
Super awesome script
Asks the user for a number:
 - If the number is less or equal to 100, it returns it to the power of itself
 - else, it returns the number squared
"""

import argparse

__version__ = '0.2'


if __name__ == '__main__':

    parser = argparse.ArgumentParser()
    parser.add_argument('--number', required=True, type=int,
                        help='number to perform calculation')
    values = parser.parse_args()
    user_number = values.number
    if user_number > 100:
        print(user_number**2)
    else:
        print(user_number**user_number)

En esta implementación hemos eliminado el problema #2 de la implementación anterior. En esta ocasión usamos argparse, de esta forma dejamos que la libreria estándar se encargue de la validación del input. Esta implementación no funciona a menos que el input sea válido.

Todavía tenemos el problema #1, la separación entre la lógica del input y la lógica primaria (la función de calculo).

En la siguiente implementación vemos como podemos arreglar esto.

Implementación 3

 
#!/usr/bin/env python

"""
Super awesome script
Asks the user for a number:
 - If the number is less or equal to 100, it returns it to the power of itself
 - else, it returns the number squared
"""

import argparse

__version__ = '0.3'



def calculation(number):
    """Performs awesome calculation"""
    if number > 100:
        return number**2
    else:
        return number**number

if __name__ == '__main__':

    parser = argparse.ArgumentParser()
    parser.add_argument('--number', required=True, type=int,
                        help='number to perform calculation')
    values = parser.parse_args()
    user_number = values.number
    calculation_result = calculation(user_number)
    print(calculation_result)

En esta implementación, hemos hecho dos cosas:

1. Hemos puesto la carga de la validación en un módulo bien mantenida como es argparse.
2. Hemos separado la lógica del input del usuario de la lógica del input de cálculo.

Éste último cambio tiene tres ventajas sobre #1 y #2.

- Ventaja 1: En primer lugar, si nos damos cuenta que por algún motivo queremos modificar el 100 por un 200, ahora podemos fácilmente modificar eso, sin tener que modificar ni leer todo el código. Siempre y cuando la función calculation siga teniendo los mismos inputs y outputs, el resto de código seguirá funcionando sin problemas.

- Ventaja 2: Otro efecto, y para mi el más significativo, es que si ahora yo leo este script que otra persona ha escrito, y me gusta tanto que quiero añadirlo a un proyecto mio, ¡ahora puedo importarlo sin problemas!.

En las implementacines #1 y #2, la única manera de usar el script era haciendo:

python calculation_script.py --number INTEGER

Ahora, en la implementación #3, tenemos una manera mucho mas útil de usar la lógica mas importante (la del cálculo). Si yo tengo otro script en el que quiero usar la funcion de cálculo, puedo usarla de la forma:

 
from calculation_script import calculation

number = 10
calculation_result = calculation(number)

¿Increíble, no? Simplemente haciendo una pequeña modificación a la estructura del proyecto, ahora cualquier persona se puede beneficiar del mismo.

- Ventaja 3: Supongamos que este simple proyecto empieza a crecer, más desarrolladores se interesan y empiezan a colaborar. El código empieza a crecer y alguien comenta que tendría sentido empezar a trabajar en el suite de testing. (si no sabes lo que es el testing, te recomiendo este artículo.)

Con la implementación #3, testear la funcionalidad de calculation es super fácil (gracias a /u/choffee en reddit por el apunte):

 
import pytest
from calculation_script import calculation

class TestCalculation:
    """Calculation function does funky things to number
    More above 100 than below
    """
    def test_zero():
        x = 0
        assert calculation(x) == 0

    def test_border():
        x = 100
        assert calculation(x) == 10000

    def test_one():
        x = 1
        assert calculation(x) == 1

Piensa en ello la próxima vez, no cuesta nada y hace que tu script sea mejor 🙂

Microentradas: Evitar ciertas etiquetas en la leyenda en Matplotlib

A veces, me llegan ficheros de datos con datos cada hora o cada día y los quiero representar en un plot. Para ello, podría acumular los ficheros en uno solo y luego pintarlo pero como lo debo hacer en 'tiempo casi-real' se puede meter todo en un bucle while que espera los ficheros cada hora/día/lo que sea y va pintando cada variable por tramos. Por ejemplo, una aproximación podría ser la siguiente:

import numpy as np
import matplotlib.pyplot as plt
plt.style.use('bmh')
%matplotlib inline

plt.figure(figsize = (12, 6))
for i in range(10):
    x = np.arange(i * 10, i * 10 + 10)
    y_var1 = np.random.randint(1, 5, 10)
    y_var2 = np.random.randint(5, 8, 10)
    plt.plot(x, y_var1, color = 'k', label = 'variable1')
    plt.plot(x, y_var2, color = 'g', label = 'variable2')
    plt.legend()
    plt.ylim(0, 9)

Como véis, en la gráfica anterior hay varios problemas pero como esta es una MicroEntrada solo nos vamos a centrar en el problema de las etiquetas repetidas en la leyenda.

¿Cómo podríamos evitar el meter tantas veces una etiqueta repetida?

Mi problema es que el bucle es o podría ser 'infinito' y tengo que inicializar las etiquetas de alguna forma. Si miro en esta respuesta encontrada en Stackoverflow dice que en la documentación se indica que "If label attribute is empty string or starts with “_”, those artists will be ignored." pero si busco aquí o en el enlace que indican en la respuesta en Stackoverflow no veo esa funcionalidad indicada en ningún sitio. Eso es porque aparecía en la versión 1.3.1 pero luego desapareció... Sin embargo podemos seguir usando esa funcionalidad aunque actualmente no esté documentada:

plt.figure(figsize = (12, 6))
for i in range(10):
    x = np.arange(i * 10, i * 10 + 10)
    y_var1 = np.random.randint(1, 5, 10)
    y_var2 = np.random.randint(5, 8, 10)
    plt.plot(x, y_var1, color = 'k', label = 'variable1' if i == 0 else "_esto_no_se_pintará")
    plt.plot(x, y_var2, color = 'g', label = 'variable2' if i == 0 else "_esto_tampoco")
    plt.legend()
    plt.ylim(0, 9)
Espero que a alguien le resulte útil.

Cómo llamar código C/C++ desde CPython (y Pypy) usando Cython y CFFI

Hace unas semanas surgió esta pregunta en StackOverflow en español: ¿Cómo llamar a código C++ desde Python?

Y la respuesta aceptada explica como hacer un wrapper sencillo usando Cython y CFFI. Como da la casualidad que la respuesta es mía voy a extenderla un poco para añadir más cosas y poder explicarla un poco mejor.

Prolegómenos

Antes de empezar a leer esta entrada deberías pasar a leer la entrada que hizo Juanlu hace un tiempo sobre CFFI titulada 'como crear extensiones en C para Python usando CFFI y Numba' donde se dan más detalles de todo el proceso a realizar con CFFI.

Antes de probar el código de la presente entrada deberías instalar cffi y cython:

conda install cffi cython # Válido en CPython

o

pip install cffi cython # Válido en CPython y Pypy

Todo lo que viene a continuación lo he probado en Linux solo usando CPython 3.5 y Pypy 5.1.1, compatible con CPython 2.7 e instalado usando esto.

Preliminares

Antes de pasar a la parte Cython y CFFI vamos a empezar creando los programas C/C++ que vamos a llamar desde Python.

Vamos a crear una librería que lo único que haga será sumar dos números enteros. Haremos una en C/C++ para Cython y una en C/C++ para CFFI.

C/C++ para Cython

C y C++ no son el mismo lenguaje pero para este caso el código se puede considerar el mismo. Para el caso C++ tendremos un fichero *.hpp y un fichero *.cpp (en C sería igual cambiando las extensiones a *.h y *.c, respectivamente).

El fichero *.hpp se llamará milibrería.hpp y contendrá el siguiente código:

long suma_enteros(long n, long m);

Mientras que el fichero *.cpp se llamará milibrería.cpp y contendrá el siguiente código:

long suma_enteros(long n, long m){
    return n + m;
}

Lo que hace el código es bastante simple.

C/C++ para CFFI

En este caso solo vamos a usar un fichero *.cpp y se llamará milibrería_cffi.cpp y contendrá el siguiente código:

long suma_enteros(long n, long m){
    return n + m;
}

extern "C"
{
    extern long cffi_suma_enteros(long n, long m)
    {
        return suma_enteros(n, m);
    }
}

El código es el mismo de antes más una segunda parte que nos permite hacer el código accesible desde Python.

Pegamento entre C/C++ y Python

En esta parte vamos a ver cómo unir el lenguaje compilado con el lenguaje interpretado.

Mediante Cython

Antes de nada necesitamos definir un fichero milibreria.pxd. Este fichero es parecido a lo que hacen los ficheros header en C/C++ o Fortran. Nos ayudará a 'encontrar' lo que hemos definido en c++ (más info sobre los ficheros pxd aquí):

cdef extern from "milibreria.hpp":
    long suma_enteros(long n, long m)

Un fichero *.pxd se puede importar en un fichero *.pyx usando la palabra clave cimport

Una vez 'enlazado' C/C++ con Cython mediante el fichero *.pxd necesitamos hacer que la parte C/C++ sea accesible desde Python. Para ello creamos el fichero pylibfromcpp.pyx, que es una especie de código Python un poco 'cythonizado' (cython es un superconjunto de Python):

cimport milibreria

def suma_enteros(n, m):
    return milibreria.suma_enteros(n, m)

Mediante CFFI

En este caso resulta un poco más sencillo, para este caso concreto. Hemos de crear el fichero Python que, mediante CFFI, enlazará C/C++ con Python. Este ficheros se llamará pylibfromCFFI.py y contendrá el siguiente código.:

import cffi


ffi = cffi.FFI()
ffi.cdef("long cffi_suma_enteros(long n, long m);")
C = ffi.dlopen("./milibreria.so")


def suma_enteros(n, m):
    return C.cffi_suma_enteros(n, m)

Setup

Compilando con Cython

Para poder acceder a la librería C/C++ hemos de crear un fichero setup.py que se encargará de la compilación que permitirá crear la extensión a la que accederemos desde Python. El fichero setup.py contendrá:

from distutils.core import setup, Extension
from Cython.Build import cythonize

ext = Extension("pylibfromcpp",
              sources=["pylibfromcpp.pyx", "milibreria.cpp"],
              language="c++",)

setup(name = "cython_pylibfromcpp",
      ext_modules = cythonize(ext))

Para crear la extensión en sí, en la misma carpeta donde hemos dejado todos los ficheros anteriores y desde la línea de comandos, hacemos (como siempre, recomiendo hacer esto desde un entorno virtual):

python setup.py build_ext -i

Y debería aparecer un fichero pylibfromcpp.cpp y otro fichero pylibfromcpp.pypy-41.so en la misma carpeta donde habéis ejecutado el comando anterior.

Compilando con CFFI

Para poder hacer accesible la funcionalidad definida en C/C++ desde Python podemos compilar usando:

g++ -o ./milibreria.so ./milibreria_cffi.cpp -fPIC -shared

Y deberíamos obtener el fichero milibreria.so.

Llamando desde Python

Usando nuestro 'wrapper' Cython

Ahora, si todo ha salido bien, dentro de un intérprete de python (como he comentado más arriba, lo he probado con CPython 3.5 y Pypy 5.1.1 y me ha funcionado en ambos) podemos hacer:

import pylibfromcpp
print(pylibfromcpp.suma_enteros(2, 3))

Usando nuestro 'wrapper' CFFI

De igual forma, si todo ha salido bien, podemos hacer:

import pylibfromcpp
print(pylibfromcpp.suma_enteros(2, 3))

Output completo en la consola pypy

Para el caso Cython

Python 2.7.10 (b0a649e90b6642251fb4a765fe5b27a97b1319a9, May 05 2016, 17:21:19)
[PyPy 5.1.1 with GCC 4.9.2] on linux2
Type "help", "copyright", "credits" or "license" for more information.
>>>> import pylibfromcpp
>>>> print(pylibfromcpp.suma_enteros(2, 3))
5

Para el caso CFFI

Python 2.7.10 (b0a649e90b6642251fb4a765fe5b27a97b1319a9, May 05 2016, 17:21:19)
[PyPy 5.1.1 with GCC 4.9.2] on linux2
Type "help", "copyright", "credits" or "license" for more information.
>>>> import pylibfromCFFI
>>>> print(pylibfromCFFI.suma_enteros(2, 3))
5

Comentarios finales

Un esquema, grosso modo, de lo que hemos hecho:

Esquema Cython - CFFI
Esquema Cython - CFFI

Pros y contras de cada una de las aproximaciones:

  • Cython permite usar Numpy sin problemas en CPython. Sin embargo, la última vez que intenté usar código Python con numpy arrays (Cython + Numpypy) reventaba todo en Pypy.
  • Cython lo podemos usar con CPython 2.x y 3.x. Cython funciona sin problemas en Pypy 5.1.1 (compatible con CPython 2.7). Numpypy NO funciona en Pypy3k.
  • El wrapper Cython que hemos hecho en este ejercicio es claramente más complejo que el que hemos hecho con CFFI (en este caso concreto).
  • Con Cython podemos usar el código compilado sin tocarlo mientras que con CFFI hemos de crear algo de código (muy simple) en el lenguaje compilado para acceder a su funcionalidad.
  • CFFI permite usar numpy arrays de forma sencilla, aunque, como con Cython, hay que 'ayudar con algo de código no Python' para que todo se pueda comunicar correctamente.

Documentación

Cython.

CFFI.

Joyitas en la stdlib: pathlib

El otro día estuvimos hablando de la biblioteca collections, una joya dentro de la librería estándar. Hoy vamos a hablar de una nueva biblioteca que se incluyó en la versión 3.4 de CPython llamada pathlib.

Solo python 3, actualízate!!!

Esta biblioteca nos da la posibilidad de usar clases para trabajar con las rutas del sistema de ficheros con una serie de métodos muy interesantes.

Algunas utilidades para configurar el problema

Vamos a crear un par de funciones que nos permiten crear y borrar un directorio de pruebas para poder reproducir el ejemplo de forma sencilla:

import os
import glob
import shutil
from random import randint, choice, seed
from string import ascii_letters

# función que nos crea un directorio de prueba en
# el mismo directorio del notebook
def crea_directorio():
    seed(1)
    base = os.path.join(os.path.curdir,
                        'pybonacci_probando_pathlib')
    os.makedirs(base, exist_ok = True)

    for i in range(0, randint(3, 5)):
        folder = ''.join([choice(ascii_letters) for _ in range(4)])
        path = os.path.join(base, folder)
        os.makedirs(path, exist_ok = True)
        for j in range(0, randint(2, 5)):
            ext = choice(['.txt', '.py', '.html'])
            name = ''.join([choice(ascii_letters) for _ in range(randint(5, 10))])
            filename = name + ext
            path2 = os.path.join(path, filename)
            open(path2, 'w').close()

# Función que nos permite hacer limpieza            
def borra_directorio():
    base = os.path.join(os.path.curdir,
                        'pybonacci_probando_pathlib')
    shutil.rmtree(base + os.path.sep)

Si ahora ejecutamos la función crea_directorio:

crea_directorio()

Nos debería quedar una estructura parecida a lo siguiente:

pybonacci_probando_pathlib/
├── KZWe
│   ├── CrUZoLgubb.txt
│   ├── IayRnBUbHo.txt
│   ├── WCEPyYng.txt
│   └── yBMWX.py
├── WCFJ
│   ├── GBGQmtsLFG.html
│   ├── PglOUshVv.py
│   └── RoWDsb.py
└── zLcE
    ├── AQlxJSXR.html
    ├── fCQGgXk.html
    └── xFUbEctT.html


Ejemplo usando lo disponible hasta hace poco

Pensemos en un problema que consiste en identificar todos los ficheros .py disponibles en determinada ruta y dejarlos en una nueva carpeta, que llamaremos python, todos juntos eliminándolos de la carpeta original en la que se encuentren.

De la forma antigua esto podría ser así:

# Suponemos que ya has creado los directorios y ficheros
# de prueba usando crea_directorio()

# recolectamos todos los ficheros *.py con sus rutas
base = os.path.join(os.path.curdir,
                    'pybonacci_probando_pathlib')
ficheros_py = glob.glob(os.path.join(base, '**', '*.py'))

# creamos la carpeta 'python' 
# dentro de 'pybonacci_probando_pathlib'
os.makedirs(os.path.join(base, 'python'), exist_ok = True)

# y movemos los ficheros a la nueva carpeta 'python'
for f in ficheros_py:
    fich = f.split(os.path.sep)[-1]
    shutil.move(f, os.path.join(base, 'python'))

Nuestra nueva estructura de ficheros debería ser la siguiente:

pybonacci_probando_pathlib/
├── KZWe
│   ├── CrUZoLgubb.txt
│   ├── IayRnBUbHo.txt
│   └── WCEPyYng.txt
├── python
│   ├── PglOUshVv.py
│   ├── RoWDsb.py
│   └── yBMWX.py
├── WCFJ
│   └── GBGQmtsLFG.html
└── zLcE
    ├── AQlxJSXR.html
    ├── fCQGgXk.html
    └── xFUbEctT.html

En el anterior ejemplo hemos tenido que usar las bibliotecas glob, os y shutil para poder realizar una operación relativamente sencilla. Esto no es del todo deseable porque he de conocer tres librerías diferentes y mi cabeza no da para tanto.

Limpieza

Me cargo la carpeta pybonacci_probando_pathlib para hacer un poco de limpieza:

borra_directorio()

Y vuelvo a crear la estructura de ficheros inicial:

crea_directorio()

Después de la limpieza vamos a afrontar el problema usando pathlib.

El mismo ejemplo con pathlib

Primero importamos la librería y, como bonus, creamos una función que hace lo mismo que la función borra_directorio pero usando pathlib, que llamaremos borra_directorio_pathlib:

from pathlib import Path

def borra_directorio_pathlib(path = None):
    if path is None:
        p = Path('.', 'pybonacci_probando_pathlib')
    else:
        p = path
    for i in p.iterdir():
        if i.is_dir():
            borra_directorio_pathlib(i)
        else:
            i.unlink()
    p.rmdir()

La anterior función con shutil es un poco más sencilla que con pathlib. Esto es lo único que hecho de menos en pathlib, algunas utilidades de shutil que vendrían muy bien de serie. Algo negativo tenía que tener.

En la anterior función, borra_directorio_pathlib, podemos ver ya algunas cositas de pathlib.

p = Path('.', 'pybonacci_probando_pathlib') nos crea una ruta que ahora es un objeto en lugar de una cadena. Dentro del bucle usamos el método iterdir que nos permite iterar sobre los directorios de la ruta definida en el objeto p. el iterador nos devuelve nuevos objetos que disponen de métodos como is_dir, que nos permite saber si una ruta se refiere a un directorio, o unlink, que nos permite eliminar el fichero o enlace. Por último, una vez que no tenemos ficheros dentro del directorio definido en p podemos usar el método rmdir para eliminar la carpeta.

Ahora veamos cómo realizar lo mismo que antes usando pathlib, es decir, mover los ficheros .py a la carpeta python que hemos de crear.

# recolectamos todos los ficheros *.py con sus rutas
p = Path('.', 'pybonacci_probando_pathlib')
ficheros_py = p.glob('**/*.py')

# creamos la carpeta 'python' dentro de 'pybonacci_probando_pathlib'
(p / 'python').mkdir(mode = 0o777, exist_ok = True)

# y copiamos los ficheros a la nueva carpeta 'python'
for f in ficheros_py:
    target = p / 'python' / f.name
    f.rename(target)

Nuevamente, nuestra estructura de ficheros debería ser la misma que antes:

pybonacci_probando_pathlib/
├── KZWe
│   ├── CrUZoLgubb.txt
│   ├── IayRnBUbHo.txt
│   └── WCEPyYng.txt
├── python
│   ├── PglOUshVv.py
│   ├── RoWDsb.py
│   └── yBMWX.py
├── WCFJ
│   └── GBGQmtsLFG.html
└── zLcE
    ├── AQlxJSXR.html
    ├── fCQGgXk.html
    └── xFUbEctT.html

Repasemos el código anterior:
Hemos creado un objeto ruta p tal como habíamos visto antes en la función borra_directorio_pathlib. Este objeto ahora dispone de un método glob que nos devuelve un iterador con lo que le pidamos, en este caso, todos los ficheros con extensión .py. En la línea (p / 'python').mkdir(mode = 0o777, exist_ok = True) podemos ver el uso de / como operador para instancias de Path. El primer paréntesis nos devuelve una nueva instancia de Path que dispone del método mkdir que hace lo que todos esperáis. Como ficheros_py era un iterador podemos usarlo en el bucle obteniendo nuevas instancias de Path con las rutas de los ficheros python que queremos mover. en la línea donde se define target hacemos uso del atributo name,que nos devuelve la última parte de la ruta. Por último, el fichero con extensión .py definido en el Path f lo renombramos a una nueva ruta, definida en target.

Y todo esto usando una única librería!!!

Echadle un ojo a la documentación oficial para descubrir otras cositas interesantes.

Si además de usar una única librería usamos parte de la funcionalidad de shutil tenemos una pareja muy potente, pathlib + shutil.

Limpieza II

Y para terminar, limpiamos nuestra estructura de ficheros pero usando ahora la función borra_directorio_pathlib que habíamos creado pero no usado aún:

borra_directorio_pathlib()

Notas

Ya hay un nuevo PEP relacionado y aceptado.

Enjoy!!

Joyitas en la stdlib: collections

Dentro de la biblioteca estándar de Python dispones de auténticas joyas, muchas veces ignoradas u olvidadas. Es por ello que voy a empezar un breve pero intenso recorrido por algunas piezas de arte disponibles de serie.

Módulo collections

Con la ayuda de este módulo puedes aumentar las estructuras de datos típicas disponibles en Python (listas, tuplas, diccionarios,...). Veamos algunas utilidades disponibles:

ChainMap

Solo Python 3. Actualízate!!

Dicho en bruto, es un conglomerado de diccionarios (también conocidos como mappings o hash tables).

Para que puede ser útil:

Ejemplo, imaginemos que tenemos un diccionario de configuración dict_a, que posee las claves a y b, y queremos actualizar sus valores con otros pares clave:valor que están en el diccionario dict_b, que posee las claves b y c. Podemos hacer:

from collections import ChainMap

dict_a = {'a': 1, 'b': 10}
dict_b = {'b': 100, 'c': 1000}

cm = ChainMap(dict_a, dict_b)
for key, value in cm.items():
    print(key, value)
a 1
c 1000
b 10

Hemos añadido el valor de la clave c de dict_b sin necesidad de modificar nuestro diccionario original de configuración dict_a, es decir, hemos hecho un 'cambio' reversible. También podemos 'sobreescribir' las claves que están en nuestro diccionario original de configuración, dict_b variando los parámetros del constructor:

cm = ChainMap(dict_b, dict_a)
for key, value in cm.items():
    print(key, value)
b 100
a 1
c 1000

Vemos que, además de añadir la clave c, hemos sobreescrito la clave b.

Los diccionarios originales están disponibles haciendo uso del atributo maps:

cm.maps
[{'b': 100, 'c': 1000}, {'a': 1, 'b': 10}]

Ejercicio: haced un dir de cm y un dir de dict_a y veréis que los atributos y métodos disponibles son parecidos.

Más información en este hilo de stackoverflow en el que me he basado para el ejemplo anterior (¿basar y copiar no son sinónimos?).

Counter

Permite contar ocurrencias de forma simple. En realidad, su funcionalidad se podría conseguir sin problemas con algunas líneas extra de código pero ya que lo tenemos, está testeado e implementado por gente experta vamos a aprovecharnos de ello.

En la documentación oficial hay algunos ejemplos interesantes y en github podéis encontrar unos cuantos más. Veamos un ejemplo simple pero potente, yo trabajo mucho con datos meteorológicos y uno de los problemas recurrentes es tener fechas repetidas que no deberían existir (pero pasa demasiado a menudo). Una forma rápida de buscar problemas de estos en ficheros y lanzar una alarma cuando ocurra lo que buscamos, sería:

from io import StringIO
from collections import Counter

virtual_file = StringIO("""2010/01/01 2.7
2010/01/02 2.2
2010/01/03 2.1
2010/01/04 2.3
2010/01/05 2.4
2010/01/06 2.2
2010/01/02 2.2
2010/01/03 2.1
2010/01/04 2.3
""")

if Counter(virtual_file.readlines()).most_common(1)[0][1] > 1:
    print('fichero con fecha repetida')
fichero con fecha repetida

namedtuple

A veces me toca crear algún tipo de estructura que guarda datos y algunos metadatos. Una forma simple sin crear una clase ad-hoc sería usar un diccionario. Un ejemplo simple sería:

import numpy as np
import datetime as dt
from pprint import pprint

datos = {
    'valores': np.random.randn(100),
    'frecuencia': dt.timedelta(minutes = 10),
    'fecha_inicial': dt.datetime(2016, 1, 1, 0, 0),
    'parametro': 'wind_speed',
    'unidades': 'm/s'
}

pprint(datos)
{'fecha_inicial': datetime.datetime(2016, 1, 1, 0, 0),
 'frecuencia': datetime.timedelta(0, 600),
 'parametro': 'wind_speed',
 'unidades': 'm/s',
 'valores': array([-3.02664796, -0.59492715, -1.36233816, -0.27333458,  0.34971592,
        1.43105631,  1.12980511,  0.49542105,  0.37546829,  1.37230197,
       -1.00757915,  1.39334713,  0.73904326,  0.01129817,  0.12431242,
        0.4388826 , -0.49561972, -0.9777947 ,  0.6009799 ,  0.89101799,
        0.48529884,  1.80287157,  1.56321415, -0.62089358, -2.22113341,
       -0.04751354,  0.89715794, -0.23252567,  0.2259216 ,  0.35214745,
       -1.50915239, -1.46547279, -0.4260315 ,  0.20851012,  1.60555432,
        0.4221521 , -1.03399518,  1.68276277,  0.5010984 ,  0.01294853,
       -0.80004557,  1.72141514, -1.38314354,  0.41374512,  0.32861028,
       -2.22385654,  0.80125671, -0.84757451,  0.66896035, -0.26901047,
       -0.06195842, -0.60743183, -0.15538184,  1.16314508, -0.42198419,
        0.61174838,  0.97211057, -1.19791368, -0.68773007,  2.96956504,
       -1.13000346, -0.24523032,  1.6312053 ,  0.77060561, -1.69925633,
       -0.31417013,  0.44196826, -0.59763569,  0.91595894,  1.47587324,
        0.5520219 , -0.62321715,  0.32543574, -1.26181508,  0.94623275,
       -0.25690824,  1.36108942,  0.15445091, -1.25607974,  0.50635589,
        0.65698443, -0.82418166, -0.34054522,  0.23511397, -1.5096761 ,
       -1.12291338, -1.82440698, -0.47433931, -1.86537903,  1.29256869,
        1.78898905,  0.72081117, -0.15169929, -1.24106944,  0.68920997,
        0.36932816, -1.15901835, -0.93990956,  0.37258685, -0.41316085])}

Lo anterior es simple y rápido pero usando una namedtuple dispongo de algo parecido con algunas cosas extra. Veamos un ejemplo similar usando namedtuple:

from collections import namedtuple

Datos = namedtuple('Datos', 'valores frecuencia fecha_inicial parametro unidades')

datos = Datos(np.random.randn(100), 
              dt.timedelta(minutes = 10),
              dt.datetime(2016, 1, 1, 0, 0),
              'wind_speed',
              'm/s')
print(datos)
Datos(valores=array([ 1.50377059, -1.48083897, -0.76143985,  0.15346996, -0.01094251,
        0.42117233,  1.07136364, -0.24586714,  1.2001748 ,  0.56880926,
        0.56959121,  0.63811853,  0.4621489 ,  1.06636058,  0.32129287,
        2.42264145, -1.25830559, -0.27102862,  2.04853711,  2.07166845,
       -0.27138347, -0.07075163, -0.43547714,  1.69140984,  2.57150371,
        0.80336641, -0.78767876, -2.22281324,  0.23112338, -0.0605485 ,
        0.58304378,  3.33116997, -1.1285789 , -0.2047658 , -0.39240644,
       -1.69724959, -0.0313781 , -0.22892613, -0.06029154, -0.32368036,
       -0.12969429,  1.06231438,  0.05429922, -1.12206555,  1.33383161,
        0.92582424,  0.51615352,  0.93188459,  0.65273332,  0.39108396,
        1.56345696, -0.33158622, -0.27455745,  0.69101563,  1.61244861,
        0.7961402 ,  0.38661924, -0.99864208, -0.10720116,  0.40919342,
       -0.43784138, -3.06455306,  1.69280852,  1.82180641,  0.03604298,
        0.17515747,  1.4370723 , -0.47437528,  1.14510249,  1.36360776,
        0.34575948, -0.14623582,  1.1048332 , -0.2266261 ,  1.34319382,
        0.75608216, -0.62416011, -0.27821722,  0.45365802, -0.98537653,
        0.20172051,  1.70476797,  0.55529542, -0.07833625, -0.62619796,
       -0.02892921, -0.07349236,  0.94659497,  0.20823509,  0.91628769,
       -1.14603843, -0.20748714,  1.13008222, -0.93365802, -0.48125316,
        0.45564591, -0.03136778, -0.86333962,  1.04590165, -0.51757806]), frecuencia=datetime.timedelta(0, 600), fecha_inicial=datetime.datetime(2016, 1, 1, 0, 0), parametro='wind_speed', unidades='m/s')

Ventajas que le veo con respecto a lo anterior:

  • Puedo acceder a los 'campos' o claves del diccionario usando dot notation
print(datos.valores)
[ 1.50377059 -1.48083897 -0.76143985  0.15346996 -0.01094251  0.42117233
  1.07136364 -0.24586714  1.2001748   0.56880926  0.56959121  0.63811853
  0.4621489   1.06636058  0.32129287  2.42264145 -1.25830559 -0.27102862
  2.04853711  2.07166845 -0.27138347 -0.07075163 -0.43547714  1.69140984
  2.57150371  0.80336641 -0.78767876 -2.22281324  0.23112338 -0.0605485
  0.58304378  3.33116997 -1.1285789  -0.2047658  -0.39240644 -1.69724959
 -0.0313781  -0.22892613 -0.06029154 -0.32368036 -0.12969429  1.06231438
  0.05429922 -1.12206555  1.33383161  0.92582424  0.51615352  0.93188459
  0.65273332  0.39108396  1.56345696 -0.33158622 -0.27455745  0.69101563
  1.61244861  0.7961402   0.38661924 -0.99864208 -0.10720116  0.40919342
 -0.43784138 -3.06455306  1.69280852  1.82180641  0.03604298  0.17515747
  1.4370723  -0.47437528  1.14510249  1.36360776  0.34575948 -0.14623582
  1.1048332  -0.2266261   1.34319382  0.75608216 -0.62416011 -0.27821722
  0.45365802 -0.98537653  0.20172051  1.70476797  0.55529542 -0.07833625
 -0.62619796 -0.02892921 -0.07349236  0.94659497  0.20823509  0.91628769
 -1.14603843 -0.20748714  1.13008222 -0.93365802 -0.48125316  0.45564591
 -0.03136778 -0.86333962  1.04590165 -0.51757806]
  • Puedo ver el código usado para crear la estructura de datos usando verbose = True. Usa exec entre bambalinas (o_O). Puedo ver que todas las claves se transforman en property's. Puedo ver que se crea documentación... MAGIA en estado puro!!!

(Si no quieres usar la keyword verbose = True puedes seguir teniendo acceso en un objeto usando obj._source)

Datos = namedtuple('Datos', 'valores frecuencia fecha_inicial parametro unidades', verbose = True)
from builtins import property as _property, tuple as _tuple
from operator import itemgetter as _itemgetter
from collections import OrderedDict

class Datos(tuple):
    'Datos(valores, frecuencia, fecha_inicial, parametro, unidades)'

    __slots__ = ()

    _fields = ('valores', 'frecuencia', 'fecha_inicial', 'parametro', 'unidades')

    def __new__(_cls, valores, frecuencia, fecha_inicial, parametro, unidades):
        'Create new instance of Datos(valores, frecuencia, fecha_inicial, parametro, unidades)'
        return _tuple.__new__(_cls, (valores, frecuencia, fecha_inicial, parametro, unidades))

    @classmethod
    def _make(cls, iterable, new=tuple.__new__, len=len):
        'Make a new Datos object from a sequence or iterable'
        result = new(cls, iterable)
        if len(result) != 5:
            raise TypeError('Expected 5 arguments, got %d' % len(result))
        return result

    def _replace(_self, **kwds):
        'Return a new Datos object replacing specified fields with new values'
        result = _self._make(map(kwds.pop, ('valores', 'frecuencia', 'fecha_inicial', 'parametro', 'unidades'), _self))
        if kwds:
            raise ValueError('Got unexpected field names: %r' % list(kwds))
        return result

    def __repr__(self):
        'Return a nicely formatted representation string'
        return self.__class__.__name__ + '(valores=%r, frecuencia=%r, fecha_inicial=%r, parametro=%r, unidades=%r)' % self

    def _asdict(self):
        'Return a new OrderedDict which maps field names to their values.'
        return OrderedDict(zip(self._fields, self))

    def __getnewargs__(self):
        'Return self as a plain tuple.  Used by copy and pickle.'
        return tuple(self)

    valores = _property(_itemgetter(0), doc='Alias for field number 0')

    frecuencia = _property(_itemgetter(1), doc='Alias for field number 1')

    fecha_inicial = _property(_itemgetter(2), doc='Alias for field number 2')

    parametro = _property(_itemgetter(3), doc='Alias for field number 3')

    unidades = _property(_itemgetter(4), doc='Alias for field number 4')


# Lo mismo de antes
print(datos._source)
from builtins import property as _property, tuple as _tuple
from operator import itemgetter as _itemgetter
from collections import OrderedDict

class Datos(tuple):
    'Datos(valores, frecuencia, fecha_inicial, parametro, unidades)'

    __slots__ = ()

    _fields = ('valores', 'frecuencia', 'fecha_inicial', 'parametro', 'unidades')

    def __new__(_cls, valores, frecuencia, fecha_inicial, parametro, unidades):
        'Create new instance of Datos(valores, frecuencia, fecha_inicial, parametro, unidades)'
        return _tuple.__new__(_cls, (valores, frecuencia, fecha_inicial, parametro, unidades))

    @classmethod
    def _make(cls, iterable, new=tuple.__new__, len=len):
        'Make a new Datos object from a sequence or iterable'
        result = new(cls, iterable)
        if len(result) != 5:
            raise TypeError('Expected 5 arguments, got %d' % len(result))
        return result

    def _replace(_self, **kwds):
        'Return a new Datos object replacing specified fields with new values'
        result = _self._make(map(kwds.pop, ('valores', 'frecuencia', 'fecha_inicial', 'parametro', 'unidades'), _self))
        if kwds:
            raise ValueError('Got unexpected field names: %r' % list(kwds))
        return result

    def __repr__(self):
        'Return a nicely formatted representation string'
        return self.__class__.__name__ + '(valores=%r, frecuencia=%r, fecha_inicial=%r, parametro=%r, unidades=%r)' % self

    def _asdict(self):
        'Return a new OrderedDict which maps field names to their values.'
        return OrderedDict(zip(self._fields, self))

    def __getnewargs__(self):
        'Return self as a plain tuple.  Used by copy and pickle.'
        return tuple(self)

    valores = _property(_itemgetter(0), doc='Alias for field number 0')

    frecuencia = _property(_itemgetter(1), doc='Alias for field number 1')

    fecha_inicial = _property(_itemgetter(2), doc='Alias for field number 2')

    parametro = _property(_itemgetter(3), doc='Alias for field number 3')

    unidades = _property(_itemgetter(4), doc='Alias for field number 4')


  • Puedo seguir obteniendo un diccionario (un OrderedDict, también incluido en el módulo collections) si así lo deseo:
datos._asdict()['valores']
array([ 1.50377059, -1.48083897, -0.76143985,  0.15346996, -0.01094251,
        0.42117233,  1.07136364, -0.24586714,  1.2001748 ,  0.56880926,
        0.56959121,  0.63811853,  0.4621489 ,  1.06636058,  0.32129287,
        2.42264145, -1.25830559, -0.27102862,  2.04853711,  2.07166845,
       -0.27138347, -0.07075163, -0.43547714,  1.69140984,  2.57150371,
        0.80336641, -0.78767876, -2.22281324,  0.23112338, -0.0605485 ,
        0.58304378,  3.33116997, -1.1285789 , -0.2047658 , -0.39240644,
       -1.69724959, -0.0313781 , -0.22892613, -0.06029154, -0.32368036,
       -0.12969429,  1.06231438,  0.05429922, -1.12206555,  1.33383161,
        0.92582424,  0.51615352,  0.93188459,  0.65273332,  0.39108396,
        1.56345696, -0.33158622, -0.27455745,  0.69101563,  1.61244861,
        0.7961402 ,  0.38661924, -0.99864208, -0.10720116,  0.40919342,
       -0.43784138, -3.06455306,  1.69280852,  1.82180641,  0.03604298,
        0.17515747,  1.4370723 , -0.47437528,  1.14510249,  1.36360776,
        0.34575948, -0.14623582,  1.1048332 , -0.2266261 ,  1.34319382,
        0.75608216, -0.62416011, -0.27821722,  0.45365802, -0.98537653,
        0.20172051,  1.70476797,  0.55529542, -0.07833625, -0.62619796,
       -0.02892921, -0.07349236,  0.94659497,  0.20823509,  0.91628769,
       -1.14603843, -0.20748714,  1.13008222, -0.93365802, -0.48125316,
        0.45564591, -0.03136778, -0.86333962,  1.04590165, -0.51757806])
  • Puedo crear subclases de forma simple para añadir funcionalidad. Por ejemplo, creamos una nueva clase con un nuevo método que calcula la media de los valores:
class DatosExtendidos(Datos):
    def media(self):
        "Calcula la media de los valores."
        return self.valores.mean()

datos_ext = DatosExtendidos(**datos._asdict())

print(datos_ext.media())
0.27764229179

deque

Otra joyita que quizá debería usar más a menudo sería deque. Es una secuencia mutable (parecido a una lista), pero con una serie de ventajas. Es una cola/lista cuyo principio y fin es 'indistinguible', es thread-safe y está diseñada para poder insertar y eliminar de forma rápida en ambos extremos de la cola (ahora veremos qué significa todo esto). Un uso evidente es el de usar, por ejemplo, una secuencia como stream de datos con un número de elementos fijo y/o rápidamente actualizable:

  • Podemos limitar su tamaño y si añadimos elementos por un lado se eliminan los del otro extremo.
  • Podemos rotar los datos de forma eficiente.
  • ...

Veamos un ejemplo:

from collections import deque

dq = deque(range(10), maxlen = 10)
lst = list(range(10))
print(dq)
print(lst)
deque([0, 1, 2, 3, 4, 5, 6, 7, 8, 9], maxlen=10)
[0, 1, 2, 3, 4, 5, 6, 7, 8, 9]
# los tres últimos elementos los anexa nuevamente al principio de la secuencia.
dq.rotate(3)
print(dq)

lst = lst[-3:] + lst[:-3]
print(lst)
deque([7, 8, 9, 0, 1, 2, 3, 4, 5, 6], maxlen=10)
[7, 8, 9, 0, 1, 2, 3, 4, 5, 6]

Veamos la eficiencia de esta operación:

tmp = deque(range(100000), maxlen = 100000)
%timeit dq.rotate(30000)
tmp = list(range(100000))
%timeit tmp[-30000:] + tmp[:-30000]
The slowest run took 9.62 times longer than the fastest. This could mean that an intermediate result is being cached.
1000000 loops, best of 3: 519 ns per loop
100 loops, best of 3: 3.07 ms per loop

Con una queue podemos anexar de forma eficiente a ambos lados:

dq.append(100)
print(dq)
dq.appendleft(10000)
print(dq)
deque([8, 9, 0, 1, 2, 3, 4, 5, 6, 100], maxlen=10)
deque([10000, 8, 9, 0, 1, 2, 3, 4, 5, 6], maxlen=10)
dq.extend(range(10))
print(dq)
dq.extendleft([10, 100])
print(dq)
deque([0, 1, 2, 3, 4, 5, 6, 7, 8, 9], maxlen=10)
deque([100, 10, 0, 1, 2, 3, 4, 5, 6, 7], maxlen=10)

Etc.

Puedes hacer cosas similares a las hechas con listas pero de forma más eficiente y práctica en determinados casos!!

Recordad que, además, disponemos del módulo queue en la librería estándar.

Conclusión

Este módulo esconde cosas muy interesantes, algunas que no hemos visto. Por tanto, si no lo conocéis, deberíais explorar el módulo collections, si lo conocéis nos podéis indicar como lo usáis en los comentarios que puedes encontrar más abajo.

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

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

Requerimientos

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

Al lío

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

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

mkdir $HOME/pypy50

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

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

Lo desempaquetamos:

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

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

mkdir $HOME/bin

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

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

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

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

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

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

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

virtualenv -p pypy $HOME/bin/pypyvenv

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

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

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

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

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

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

¡¡¡Listo!!!

Usando pypy

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

pypy

Usando pypy con numpy en un notebook de jupyter

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

. ~/bin/pypyvenv/bin/activate

Y arrancamos jupyter:

jupyter notebook

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

Ideas para mejorar el script (con vuestros pull requests)

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

Saludos.

 

Fórmula para el amor

Esta entrada se proyectó hace unos doscientos cuarenta y pico días.

Vamos a representar la siguiente fórmula:

\({x}^2 + (y - \sqrt{x^2})^2 = 1\)

Si despejamos la \(y\) nos quedarán las siguientes soluciones:

\(y_{1} = \sqrt{x^2} + \sqrt{1 - x^2}\)
\(y_{2} = \sqrt{x^2} - \sqrt{1 - x^2}\)

En código Python usando Numpy y Matplotlib tendremos lo siguiente:

import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
x = np.linspace(-1,1,50)
y1 = np.sqrt(x * x) + np.sqrt(1 - x * x)
y2 = np.sqrt(x * x) - np.sqrt(1 - x * x)
plt.plot(x, y1, c='r', lw = 3)
plt.plot(x, y2, c='r', lw = 3)
plt.show()

Felicidades a quien corresponda.

Idea copiada literalmente de aquí.

Cómo crear extensiones en C para Python usando CFFI y numba

Introducción

En este artículo vamos a ver cómo crear extensiones en C para Python usando CFFI y aceleradas con numba. El proyecto CFFI ("C Foreign Function Interface") pretende ofrecer una manera de llamar a bibliotecas escritas en C desde Python de una manera simple, mientras que numba, como podéis leer en nuestro blog, es un compilador JIT para código Python numérico. Mientras que hay algo de literatura sobre cómo usar CFFI, muy poco se ha escrito sobre cómo usar funciones CFFI desde numba, una característica que estaba desde las primeras versiones pero que no se completó hasta hace cuatro meses. Puede parecer contradictorio mezclar estos dos proyectos pero en seguida veremos la justificación y por qué hacerlo puede abrir nuevos caminos para escribir código Python extremadamente eficiente.

Este trabajo ha surgido a raíz de mis intentos de utilizar funciones hipergeométricas escritas en C desde funciones aceleradas con numba para el artículo que estoy escribiendo sobre poliastro. El resultado, si bien no es 100 % satisfactorio aún, es bastante bueno y ha sido relativamente fácil de conseguir, teniendo en cuenta que partía sin saber nada de C ni CFFI hace tres días.

¿Por qué CFFI + numba?

Como decíamos CFFI y numba, aunque tienen que ver con hacer nuestros programas más rápidos, tienen objetivos bastante diferentes:

  • CFFI nos permite usar C desde Python. De este modo, si encontramos algún algoritmo que merece la pena ser optimizado, lo podríamos escribir en C y llamarlo gracias a CFFI.
  • numba nos permite acelerar código Python numérico. Si encontramos algún algoritmo que merece la pena ser optimizado, adecentamos un poco la función correspondiente y un decorador la compilará a LLVM al vuelo.

Continue reading

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.