Pybonacci

Computación Científica con Python en castellano

Juego de la vida de Conway con Python usando NumPy

with 11 comments

Introducción

En esta breve entrada vamos a ver cómo implementar el juego de la vida de Conway en Python utilizando arrays de NumPy. El juego de la vida es posiblemente el ejemplo más famoso de autómata celular, es decir, un sistema dinámico que evoluciona en pasos discretos. Actualmente están siendo estudiados en profundidad por Stephen Wolfram (creador además del software Mathematica) y tienen aplicaciones en campos tan diversos como la biología o la economía.

Aunque existe una implementación muy inteligente en Rosetta Code utilizando objetos defaultdict, voy a escribir la mía propia porque me parece «visualmente» más intuitiva. Para eso me voy a basar en el algoritmo seguido por Michael Gertelman[I] para escribir el programa en una línea usando APL (!) y en una entrada de Peter Collingridge[II], basada a su vez en la primera, donde aplica la misma idea en Python. Ambos están en las referencias al final de este texto.

¡Vamos allá!

Nota: Jake VanderPlas ofrece una versión todavía más corta del juego de la vida. Básicamente es la misma idea pero resumida en una línea :) La otra versión emplea la función `scipy.signal.convolve2d`.

El juego de la vida de Conway

Aunque hay muchas variantes del juego de la vida, nosotros nos vamos a centrar en la que publicó originalmente Conway en 1970. El juego de la vida es un juego de cero jugadores, lo que quiere decir que su evolución está determinada por el estado inicial y no necesita ninguna entrada de datos posterior (Wikipedia). Las reglas son las siguientes:

  • El tablero es una rejilla de celdas cuadradas que tienen dos posibles estados: viva y muerta.
  • Cada célula tiene ocho células vecinas: se cuentan también las de las diagonales.
  • En cada paso, todas las células se actualizan instantáneamente teniendo en cuenta la siguiente regla:
    • Cada célula viva con 2 o 3 células vecinas vivas sobrevive.
    • Cada célula con 4 o más células vecinas vivas muere por superpoblación. Cada célula con 1 o ninguna célula vecina viva muere por soledad.
    • Cada célula muerta con 3 células vecinas vivas nace.

¡Vamos con la implementación!

Código Python

Para recrear el tablero voy a utilizar un array de NumPy de tipo int, donde el valor 0 corresponderá a célula muerta y el 1 a célula viva. La clave para calcular el número de células vivas va a estar en la función roll, que toma un array de NumPy y desplaza todos los elementos en la dirección indicada. Veamos un ejemplo:

>>> b = np.diag([1, 2, 3])
>>> b
array([[1, 0, 0],
       [0, 2, 0],
       [0, 0, 3]])
>>> np.roll(b, 1, axis=0)  # Eje 0: filas
array([[0, 0, 3],
       [1, 0, 0],
       [0, 2, 0]])
>>> np.roll(b, 2, axis=0)  # Dos posiciones
array([[0, 2, 0],
       [0, 0, 3],
       [1, 0, 0]])
>>> np.roll(b, -1, axis=0)  # Una posición hacia atrás
array([[0, 2, 0],
       [0, 0, 3],
       [1, 0, 0]])
>>> np.roll(b, 1, axis=1)  # Eje 1: columnas
array([[0, 1, 0],
       [0, 0, 2],
       [3, 0, 0]])

Pues bien, vamos a desplazar nuestro tablero una posición en las ocho direcciones posibles (recuerda que se cuenta la diagonal también) y vamos a sumar el resultado. Como cada célula viva tendrá el valor 1, esto me dará, en cada celda, el número de células vecinas vivas.

def vecindario(b):
    """Array de células vivas en el vecindario."""
    vecindario = (
        np.roll(np.roll(b, 1, 1), 1, 0) +  # Arriba-izquierda
        np.roll(b, 1, 0) +  # Arriba
        np.roll(np.roll(b, -1, 1), 1, 0) +  # Arriba-derecha
        np.roll(b, -1, 1) +  # Derecha
        np.roll(np.roll(b, -1, 1), -1, 0) +  # Abajo-derecha
        np.roll(b, -1, 0) +  # Abajo
        np.roll(np.roll(b, 1, 1), -1, 0) +  # Abajo-izquierda
        np.roll(b, 1, 1)  # Izquierda
    )
    return vecindario

date cuenta de que para desplazar en diagonal tenemos que aplicar la función roll dos veces.

Nota: Los elementos dan la vuelta en la frontera del tablero, así que será como si considerásemos que es «esférico». En otras versiones del juego se considera siempre que las células de fuera del tablero están muertas.

Una vez que tenemos el array vecindario, es sencillísimo determinar qué células sobreviven, cuáles mueren y cuáles nacen:

def paso(b):
    """Paso en el juego de la vida de Conway."""
    v = vecindario(b)  # Se calcula el vecindario
    buffer_b = b.copy()  # Hacemos una copia de la matriz
    for i in range(buffer_b.shape[0]):
        for j in range(buffer_b.shape[1]):
            if v[i, j] == 3 or (v[i, j] == 2 and buffer_b[i, j]):
                buffer_b[i, j] = 1
            else:
                buffer_b[i, j] = 0
    return buffer_b

Y el tablero en el nuevo paso será el resultado de la función anterior.

Y esto ya está, casi la parte más complicada del código es conseguir crear una animación a partir de aquí. Como llevo un rato frustrado, y aunque podría arreglar esto de manera limpia leyendo un rato más el tutorial de Jake Vanderplas sobre matplotlib.animation, voy a hacer algo de lo que no me siento orgulloso y nunca, nunca tienes que intentar hacer en tu casa: usar una variable global. Si el gran KvdP lo hizo, yo también puedo :P

Animación del juego de la vida de Conway

Animación del juego de la vida de Conway

Y este es el código:

Espero que os haya gustado, ¡un saludo!

P.D.: Perdonad la ligera inactividad y las prisas en escribir esta entrada, pero estamos muy ocupados organizando la que será la primera PyCon en España. Más información en #PyConES ;)

Referencias

  1. GERTELMAN, Michael. Conway’s Game of Life in one line of APL [en línea]. Disponible en Web: <http://catpad.net/michael/apl/>. [Consulta: 30 de noviembre de 2012]
  2. COLLINGRIDGE, Peter. Game of Life in one line of Python [en línea]. Disponible en Web: <http://www.petercollingridge.co.uk/blog/python-game-of-life-in-one-line>. [Consulta: 30 de noviembre de 2012]
  3. GARDNER, Martin. Mathematical Games – The fantastic combinations of John Conway’s new solitaire game “life” [en línea]. Disponible en Web: <http://ddi.cs.uni-potsdam.de/HyFISCH/Produzieren/lis_projekt/proj_gamelife/ConwayScientificAmerican.htm>. [Consulta: 30 de noviembre de 2012]
About these ads

Written by Juanlu001

30 de noviembre de 2012 a 10:26

11 comentarios

Suscríbete a los comentarios mediante RSS.

  1. Yo te apoyo con la variable global, no hay que ser un taliban del código :P

    El zombi de Schrödinger

    30 de noviembre de 2012 at 12:21

  2. Un artículo muy interesante, sobre la forma de resolverlo con numpy y matplot.

    Pero ya perdonarás si discrepo en el uso de variables globales. No es por el hecho de usarlas, sino porque “mutan”. Si realmente necesitas un objeto que cambie de estado, ¿no estaría todo mejor dentro de una clase?

    Otro detalle: los comentarios de la función vencidario son un poco confusos. Si estás calculando el número de vecinos, con np.roll(b, 1, 0) estarías mirando si hay vecino “arriba”. Creo que quedaría más claro el cálculo si los comentarios fueran acorde con el vecino que estás mirando y no con el movimiento que le das al tablero.

    Y por último proponer un ejercicio: si las operaciones np.roll son operaciones con matrices: ¿no sería más sencillo calcular una matriz generadora ‘G’ que calcule los vecinos multiplicando ‘G*tablero’ y luego aplicar una expresión lógica para obtener el siguiente ciclo? (Lo digo sin pensar mucho. Igual es una tontería :-P)

    Chema Cortés

    14 de agosto de 2013 at 16:29

    • ¡Gracias por el comentario Chema!

      Jeje aún me queda mucho para adquirir tu mentalidad funcional… :P Este código lo escribí en dos patadas y la función animate ahí donde la ves es un petardo, así que la variable global me sacó de apuros rápidamente.

      (http://xkcd.com/292/ s/goto/global/g)

      Con lo de los comentarios tienes toda la razón. En cuanto a lo de la matriz, no lo veo muy directo; se me ocurre que tal vez podría escribirse un array de tres o cuatro dimensiones, hacer un producto y luego contraer una o dos (usando einsum), pero necesitaría un rato para pensarlo (y a estas horas desde luego no me va a salir). Con una matriz normal desde luego es imposible, o eso me parece a mí (una prueba matemática de esto también me requeriría algunas horas de sueño primero).

      Al final, el método de Jake es como comprimir lo mismo que he escrito yo: hacer una lista por comprensión doble y sumar todos sus elementos. El otro, el de hacer una convolución, no entiendo cómo funciona porque no me acuerdo de cómo funcionan las convoluciones :)

      Juanlu001

      14 de agosto de 2013 at 22:22

      • Estoy vago de vacaciones, pero voy a intentar explicarme algo más:

        Agrupemos la operaciones roll en estas dos matrices:

        A = np.roll(np.identity(N), 1, 0) + np.roll(np.identity(N), -1, 0)
        B = np.roll(np.identity(M), 1, 1) + np.roll(np.identity(M), -1, 1)

        La matriz “vencindario” sería: A*T + T*B + A*T*B , o sea,

        V = np.dot(A, tablero) + np.dot(tablero, B) + np.dot(np.dot(A, tablero), B)

        Se puede acortar a dos operandos, pero ya se ve que son menos operaciones que las 8 que tenías en la función vencindario.

        Para calcular el siguiente ciclo, es más “elegante” algo así:

        tablero = np.select([V==3, V+tablero==3], [1,1], 0)

        Ya dirás si funciona :-P

        Chema Cortés

        16 de agosto de 2013 at 13:26

      • Tu idea funciona, ¡pero tengo algo mejor! :D

        Si escribimos el tablero como un grafo donde cada nodo es una célula y cada arista indica células adyacentes según nuestras reglas, su matriz de adyacencia A nos da el número de vecinos vivos, sin más que hacer

        \displaystyle A v_T = v_V

        siendo v_T el vector columna resultante de aplanar el array T.

        Por ejemplo, para este tablero 3×3:

        \displaystyle \begin{pmatrix} 1 & 1 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \end{pmatrix}

        la matriz de adyacencia (dada solo por la geometría) será:

        \begin{pmatrix} 0 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 \\ 1 & 0 & 1 & 1 & 1 & 1 & 1 & 1 & 1 \\ 1 & 1 & 0 & 1 & 1 & 1 & 1 & 1 & 1 \\ 1 & 1 & 1 & 0 & 1 & 1 & 1 & 1 & 1 \\ 1 & 1 & 1 & 1 & 0 & 1 & 1 & 1 & 1 \\ 1 & 1 & 1 & 1 & 1 & 0 & 1 & 1 & 1 \\ 1 & 1 & 1 & 1 & 1 & 1 & 0 & 1 & 1 \\ 1 & 1 & 1 & 1 & 1 & 1 & 1 & 0 & 1 \\ 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 0 \\ \end{pmatrix}

        (teniendo en cuenta que estamos usando una geometría toroidal). El vector será:

        v_T = \begin{pmatrix} 1 \\ 1 \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ \end{pmatrix}

        y por tanto el vecindario:

        A v_T = \begin{pmatrix} 1 \\ 1 \\ 2 \\ 2 \\ 2 \\ 2 \\ 2 \\ 2 \\ 2 \\ \end{pmatrix}

        y por tanto el array correspondiente:

        V = \begin{pmatrix} 1 & 1 & 2 \\ 2 & 2 & 2 \\ 2 & 2 & 2 \\ \end{pmatrix}

        ¡Y fin! :D

        Tan solo me falta escribir la regla para construir la matriz de adyacencia (NetworkX ayuda) y comprobar si este método es más rápido que los dos que propone Jake VanderPlas.

        Juanlu001

        17 de agosto de 2013 at 11:45

  3. En definitiva, este artículo está lleno de errores, omisiones y posibles mejoras. No sé yo si tirarlo a la basura y empezar a escribirlo desde cero xD

    Juanlu001

    17 de agosto de 2013 at 11:50

    • La mejoras siempre son positivas, pero nunca tires a la basura el proceso si no quieres perder el contexto.

      He escrito el código tal como lo había pensado: https://gist.github.com/chemacortes/6257951#file-life-py
      No poca experiencia en animaciones con matplotlib, pero parece que funciona bastante bien.

      Chema Cortés

      17 de agosto de 2013 at 18:48

      • La mejoras siempre son positivas, pero nunca tires a la basura el proceso si no quieres perder el contexto.

        Me lo apunto :)

        Juanlu001

        17 de agosto de 2013 at 18:53

    • Resultados preliminares de mis pesquisas: éxito

      http://nbviewer.ipython.org/6260832

      Juanlu001

      18 de agosto de 2013 at 10:57

      • Muy interesante. Aunque se diría que aumenta la complejidad algoritmica, numpy/scipy se las apaña bastante bien.

        Sólo una cosa más: creo que se puede hacer más claro el cálculo de la matriz de adjaciencia y más sencillo así:

        
        row_11 = np.zeros((M,N),dtype=int)
        row_11[0:3,0:3] = block_ij
        row_00 = np.roll(row_11.flatten(), -N - 1)
        

        Chema Cortés

        19 de agosto de 2013 at 18:03


¡Deja un comentario!

Introduce tus datos o haz clic en un icono para iniciar sesión:

Logo de WordPress.com

Estás comentando usando tu cuenta de WordPress.com. Cerrar sesión / Cambiar )

Imagen de Twitter

Estás comentando usando tu cuenta de Twitter. Cerrar sesión / Cambiar )

Foto de Facebook

Estás comentando usando tu cuenta de Facebook. Cerrar sesión / Cambiar )

Google+ photo

Estás comentando usando tu cuenta de Google+. Cerrar sesión / Cambiar )

Conectando a %s

Seguir

Recibe cada nueva publicación en tu buzón de correo electrónico.

Únete a otros 1.103 seguidores

%d personas les gusta esto: