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.

SciPy Latin America 2015 en Argentina: ¡manda tu charla!

Ya lo hemos comentado por Twitter algunas veces, pero para quienes aún no lo sepáis: del 20 al 22 de mayo se celebrará la SciPy Latin America 2015 en Posadas, Misiones (Argentina), y gracias a la generosa invitación de los organizadores tanto Kiko como el que os escribe iremos en calidad de speakers invitados. ¡Muchísimas gracias! :D

SciPy Latin America 2015

Este evento surge después de dos años de celebrar conferencias sobre Python científico en Argentina, la primera en 2013 en Puerto Madryn, Chubut y la segunda en 2014 en Bahía Blanca, Buenos Aires. Este año los organizadores han sumado fuerzas con gente de Brasil para llevar adelante un evento mucho más ambicioso a nivel de todo Latinoamérica. Desde aquí no puedo más que quitarme el sombrero por semejante esfuerzo (que incluye personas entendiéndose en portugués y español :) ) y estoy seguro de que la SciPyLA 2015 será todo un éxito.

El plazo para presentar propuestas se cierra el día 6 de abril, ¡así que no pierdas un minuto para mandar tu charla, taller o póster a la conferencia! En el llamado a propuestas de la SciPyLA 2015 tienes toda la información sobre qué tipo de formatos se aceptan. La conferencia es gratuita, abierta a todo el mundo y habrá charlas en español, inglés y portugués. Nosotros ya hemos mandado nuestras propuestas, ¿a qué esperas tú? :)

Para mí es una ocasión muy emocionante para conocer a gente con la que normalmente solo puedo hablar a través de las redes sociales, y a la vez un reconocimiento excepcional a la labor que llevamos haciendo ya tres años en Pybonacci. Tendremos el honor de coincidir con personas que han tenido un papel crucial en el desarrollo de la comunidad Python argentina, hispana en general y brasileña, como Damián Ávila, Raniere Silva, Martín Gaitán, Filipe Saraiva, Mariano Reingart, Celia Cintas o Manuel Kaufmann, nos reencontraremos con Juan Bautista después de dos años de cruzarnos en la PyConES 2013 y desvirtualizaremos a fans que nos siguen desde el primer día como Cristian Hernán Schmidt. Puedes ver la lista pública de asistentes a la SciPyLA 2015 en la web de la conferencia.

Cartel SciPyLA 2015

No te puedes perder este evento: manda ya tu charla, regístrate (que es gratis) y ¡nos vemos en Argentina! ;)

MicroEntradas: scipy.constants

¿No te acuerdas de la constante de gravitación universal?, ¿no sabes cuanta área es un acre?, ¿por qué me pasan esos resultados en pulgadas?, ¿a qué altura estamos volando cuando me dicen que estamos a 10.000 pies?,... Puedes responder a todo eso usando el, creo que infrautilizado, módulo constants dentro del paquete scipy.

Primero de todo, vamos a importat el módulo en cuestión:

from scipy import constants as constantes

En este módulo disponemos de varias constantes físicas y matemáticas de uso muy habitual en determinados campos. El número pi, la constante de gravitación universal, la constante de Plank o la masa del electrón están en la punta de tus dedos.

print(constantes.pi,
      constantes.gravitational_constant,
      constantes.Plank,
      constantes.m_e)

Nos dará el siguiente resultado:

3.141592653589793 6.67384e-11 6.62606957e-34 9.10938291e-31

Podemos acceder a otras constantes (no tan constantes) usando un diccionario con el nombre de la constante según la base de datos del Committee on Data for Science and Technology (CODATA):

constantes.value('standard atmosphere')

Podemos obtener el valor de varias unidades en el sistema internacional simplemente poniendo su nombre

print(constantes.foot, constantes.inch)

Incluso tenemos funciones para hacer conversiones de unidades

print(constantes.C2K(10))

Nos daría el valor de 10ºC en grados Kelvin.

Saludos.

 

Teoría de control en Python con SciPy (II): Control PID

Introducción

En esta serie de artículos vamos a estudiar cómo podemos aplicar Python al estudio de la teoría de control, en este caso utilizando SciPy. La teoría de control se centra en los sistemas dinámicos con entradas: sistemas físicos cuyo estado evoluciona con el tiempo en función de la información que reciben del exterior. Como puedes ver, esta definición es enormemente amplia: el control toca aspectos de la ingeniería y de las matemáticas, y tiene aplicaciones también en las ciencias sociales: psicología, sociología, finanzas...

  1. Conceptos básicos
  2. Control PID

En esta segunda parte, una vez vistos algunos conceptos básicos en la primera, vamos a retomar el problema del control de crucero del coche exactamente donde lo dejamos:

La conclusión que podemos extraer de este gráfico es que para pasar de 0 a 100 km/h, nuestro coche necesita casi un minuto. ¡Fatal! ¿Cómo arreglamos esto? Aquí entra la belleza de la teoría de control, pero lo vamos a dejar para la segunda parte ;)

A las referencias que ya recomendé en la primera parte voy a añadir [Ogata, 2010], un excelente libro sobre ingeniería de control traducido al español por profesores de España y Argentina. De ahí he estudiado la parte matemática del control PID y me he inspirado para las figuras. Cualquier aspecto teórico se puede consultar en este libro.

En esta entrada se han usado python 3.3.2, numpy 1.8.0, scipy 0.13.0 y matplotlib 1.3.0.

Concepto de realimentación

Si recordamos el modelo de nuestro coche, teníamos un sistema con una entrada (la fuerza de tracción que genera el motor) y una salida o variable a controlar (la velocidad del coche). Este tipo de sistemas se denominan en lazo abierto y se pueden esquematizar con un diagrama de bloques de este estilo:

Lazo abierto

En este caso la planta sería el motor y el controlador un sistema que consiguiese esa tracción constante que consideramos en la primera parte. Este tipo de sistemas son poco útiles porque no podemos disponer de la información de la salida para controlar la entrada. No les prestaremos más atención.

Si queremos tener en cuenta las posibles perturbaciones de la salida deberíamos medirla continuamente para comprobar que cumple con nuestros requisitos. A esto se le denomina control en lazo cerrado y se puede esquematizar de la siguiente manera:

Lazo cerrado

Ya vemos que se van complicando un poco las cosas. En este caso, la entrada de la planta la proporciona el actuador, y la entrada del actuador es la diferencia entre la señal de referencia y la salida. Esta diferencia se denomina señal error por razones obvias. Con este cambio de esquema y de filosofía tenemos en cada instante información sobre la salida del sistema y podemos por tanto ajustar el control del mismo. Veamos un método para llevar a cabo este control.

Continue reading

Teoría de control en Python con SciPy (I): Conceptos básicos

Introducción

En esta serie de artículos vamos a estudiar cómo podemos aplicar Python al estudio de la teoría de control, en este caso utilizando SciPy. La teoría de control se centra en los sistemas dinámicos con entradas: sistemas físicos cuyo estado evoluciona con el tiempo en función de la información que reciben del exterior. Como puedes ver, esta definición es enormemente amplia: el control toca aspectos de la ingeniería y de las matemáticas, y tiene aplicaciones también en las ciencias sociales: psicología, sociología, finanzas...

  1. Conceptos básicos
  2. Control PID

En esta primera parte vamos a hacer una breve introducción matemática para centrar el tema y vamos a ver el manejo básico de sistemas LTI.

Cuando uno piensa en estudiar sistemas dinámicos con un ordenador, automáticamente se le viene a la cabeza MATLAB, y no sin motivo. Este programa tiene unas capacidades extraordinarias en este campo, y aunque nos duela decirlo Python no está al mismo nivel. Sin embargo, queremos mostrar en este artículo que Python tiene el potencial de ser una alternativa real a MATLAB, enseñando los fundamentos del análisis de sistemas dinámicos utilizando el paquete scipy.signal. Yo mismo he trabajado un poco en este paquete en los últimos meses, así que he tenido la oportunidad de ver cómo funciona y también de conocer sus carencias; algunas de mis contribuciones han visto la luz en la recién liberada versión 0.13 de SciPy, pero aún queda mucho por mejorar.

Equivalencia entre los dominios del tiempo y de la frecuencia a través de la transformada de Laplace
Equivalencia entre los dominios del tiempo y de la frecuencia a través de la transformada de Laplace

Los ejemplos para este artículo los he sacado de [Sedra y Smith, 2004], un excelente libro de electrónica, y de [Messner et al. 2011], unos tutoriales para MATLAB y Simulink. Para la teoría, recomiendo el excelente [Gil y Rubio 2009], un libro editado por la Universidad de Navarra y disponible para visualización, impresión y copia para uso personal sin fines de lucro (¡gracias @Alex__S12!).

En esta entrada se han usado python 3.3.2, numpy 1.8.0, scipy 0.13.0 y matplotlib 1.3.0.

Sistemas lineales invariantes en el tiempo o LTI

Vamos a centrarnos en los sistemas lineales e invariantes en el tiempo (en inglés, sistemas LTI), que como su propio nombre indica tienen estas propiedades:

  • Linealidad: el sistema cumple el principio de superposición.
  • Invarianza en el tiempo: el comportamiento del sistema no varía con el tiempo: una misma entrada en dos instantes de tiempo diferentes siempre producirá la misma salida.

Concretamente, muchos de estos sistemas pueden modelarse como un sistema de ecuaciones diferenciales ordinarias lineales de coeficientes constantes. Restringiendo aún más nuestro objeto de estudio, para sistemas de una sola entrada y una sola salida (SISO en inglés) tendremos:

\(\displaystyle a_n \frac{d^n y}{dt^n} + a_{n-1} \frac{d^{n-1} y}{dt^{n-1}} + \dots + y(t) = b_m \frac{d^m x}{dt^m} + b_{m-1} \frac{d^{m-1} x}{dt^{m-1}} + \dots + x(t)\)

donde \(y(t)\) es la respuesta del sistema y \(x(t)\) es la entrada, conocida. Aplicando la transformada de Laplace a la ecuación tendremos:

\(\displaystyle Y(s) = \frac{b_m s^m + b_{m-1} s^{m-1} + \dots + 1}{a_n s^n + a_{n-1} s^{n-1} + \dots + 1} X(s) = H(s) X(s) \)

siendo \(Y(s)\) y \(X(s)\) las transformadas de laplace de \(y(t)\) y \(x(t)\) y \(H(s)\) la función de transferencia del sistema. Vemos que la ecuación resultante es algebraica y por tanto muy fácil de resolver.

Continue reading

Hoja de ruta para SciPy 1.0

Los desarrolladores de SciPy han anunciado, a través de la lista de correo, que se ha preparado un borrador con la hoja de ruta hacia la consolidación de SciPy 1.0.

Hoja de ruta:

https://github.com/rgommers/scipy/blob/roadmap/doc/ROADMAP.rst.txt

Anuncio en la lista de correo:

http://mail.scipy.org/pipermail/scipy-dev/2013-September/019237.html

Ralf Gommers, Pauli Virtanen y David Cournapeau se reunieron en la EuroSciPy que se celebró en Bruselas en agosto y han redactado un documento donde se recogen los pasos necesarios en términos de nueva funcionalidad, solución de fallos, etc. para llegar a SciPy 1.0. Cuando esto se consiga, significará que SciPy contiene lo esencial y que tiene una API y un código de buena calidad.

Como dice Jake Vanderplas, es una oportunidad única para participar e influir en la dirección que tomará el desarrollo de SciPy, que es hoy por hoy la piedra angular del ecosistema Python científico junto con NumPy. La gran ventaja de esta hoja de ruta es que está pormenorizado, por cada paquete, en qué estado se encuentra el código y qué áreas necesitan más trabajo, de forma que es mucho más fácil para potenciales colaboradores decidir en qué quieren ayudar.

Durante los próximos meses (aunque este proceso se puede demorar todavía algunos años) se seguirán liberando versiones en la rama 0.x donde poco a poco se irán incorporando estos cambios, que incluyen:

  • Limpieza de código antiguo y supresión de duplicidades: scipy.fftpack, scipy.misc, eliminación de weave.
  • Cambios en la API que rompen la compatibilidad hacia atrás: scipy.integrate.
  • Revisión del código para asegurar que produce el resultado correcto: scipy.stats, scipy.signal.
  • Mejora de la documentación y de la cobertura por tests: scipy.stats y en general.
  • Aumento de la eficiencia de los algoritmos.
  • Reescritura en Cython de algunas partes de SciPy: scipy.cluster, scipy.sparse, scipy.spatial.

Cuando estos objetivos se alcancen se liberará SciPy 1.0.

Si estás interesado en empezar a contribuir código a SciPy, deberías también leerte esta guía:

http://docs.scipy.org/doc/scipy-dev/reference/hacking.html

En ella explican la estructura de paquetes de SciPy, cómo organizar el código, cómo disponer un entorno de desarrollo para que sea más fácil probar tus cambios, cómo ejecutar los tests y muchas cosas más.

Desde aquí animamos a todos aquellos lectores que alguna vez hayan pensado en contribuir a SciPy a leer detenidamente el roadmap y la guía para desarrolladores y que nos cuenten qué les parece. ¿Es suficiente con esto? ¿Habéis tenido dificultades para aclararos por dónde empezar? ¿Creéis que hay algo que se podría mejorar? ¿Pensáis que esto puede ser importante para atraer nuevos desarrolladores a SciPy?

Si tenéis alguna sugerencia que hacer podemos discutirlo juntos y transmitírselo al equipo de desarrolladores :)

¡Un saludo!

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

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

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

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

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

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

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

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

¡Saludos!

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

Introducción

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

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

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

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

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

Interpolación

Polinomios no, ¡gracias!

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

Continue reading

Análisis Cluster (II): Clasificación no supervisada mediante clasificación jerárquica aglomerativa

(Este es el segundo 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 el artículo inicial échale un ojo ahora).

Como vimos anteriormente, existen diferentes formas de hacer clustering y, como también comentamos anteriormente, una de las más habituales es el clustering jerárquico.

El clustering jerárquico asociativo pretende, partiendo de m observaciones, ir encontrado agrupaciones de forma jerarquizada. Para ello, primero busca similitudes entre las observaciones y después procura asociar en grupos o 'clusters' las observaciones que se encuentran 'más cercanas' o presentan mayor similitud.

Si os acordáis, en el primer capítulo de esta mini-serie, entrecomillamos la palabra ‘similitud’. Vamos a ver qué significa esto de similitud en general y en nuestro ejemplo concreto (ver capítulo anterior de la serie para conocer el ejemplo). Dependiendo del problema concreto necesitaremos asociar las variables para poder medir como son de similares o a qué distancia se encuentran entre sí. Estas medidas de asociación (similitud o distancia) dependerán del problema concreto con el que nos encontremos y no se usaría lo mismo para variables booleanas, binarias, reales,... Gracias al módulo scipy.spatial.distance podemos ver muchas de estas medidas de asociación en las que no voy a entrar en detalle. Un resumen de la mayoría de ellas lo podéis ver en el siguiente enlace del profesor José Ángel Gallardo San Salvador, nuevamente. En el ejemplo propuesto usaremos la correlación como medida de asociación ya que es una medida invariante aunque las variables se escalen y/o se les sumen/resten parámetros (algo muy útil para el caso que nos ocupa).

Continue reading

Cómo resolver ecuaciones algebraicas en Python con SciPy

Introducción

En este artículo vamos a utilizar las rutinas de búsqueda de raíces ya disponibles en el módulo scipy.optimize para resolver ecuaciones algebraicas con Python. Ya vimos hace tiempo cómo encontrar el mínimo de una función con SciPy, y también cómo implementar los métodos de la bisección y de Newton en Python. Ahora, además, exploraremos el caso de sistemas de ecuaciones.

En esta entrada se ha usado python 2.7.3, numpy 1.6.2 y scipy 0.11.

Ejemplo básico: métodos de Brent y de Newton

Aunque el método de la bisección es muy conocido porque conceptualmente es muy sencillo de entender, nos vamos a olvidar de él y vamos a utilizar en su lugar el método de Brent, que viene implementado en SciPy en la función scipy.optimize.brentq, ya que según la documentación es la mejor opción. El método de Brent es un algoritmo complicado que combina el método de la bisección, el de la secante e interpolación cuadrática inversa: como el método de la bisección, tiene convergencia garantizada siempre que demos un intervalo \([a, b]\) donde \(f(a) f(b) < 0\); es decir, que haya un cambio de signo en el intervalo.

Continue reading