Hace meses publicábamos en Pybonacci una entrada sobre cómo implementar el Juego de la vida de Conway en Python utilizando arrays de NumPy. Recientemente el gran Jake VanderPlas hizo lo mismo, dando un código mucho más compacto y centrándose más en la creación de animaciones. A raíz de ello, Chema Cortés preguntó en nuestro blog si no existiría una manera más eficiente de calcular cada paso. En esta entrada vamos a explicar nuestros resultados, y ya os adelantamos que son muy buenos :)
Aviso a navegantes: hoy vamos a ver más matemáticas que Python, porque creo que es valioso explicar cómo hallé este nuevo método. Pero ¿a quién no le gustan las matemáticas?
(Si no te interesa esa parte, puedes saltar directamente a cómo construir la matriz de adyacencia).
En esta entrada se han utilizado numpy 1.7.1, scipy 0.12.0 y matplotlib 1.3.
Supongamos que tenemos un tablero muy pequeño, de 5 filas y 5 columnas nada más, y que hemos etiquetado las celdas. Las vamos a llamar $c_{ij}$, donde $i$ será la fila y $j$ será la columna (empezando por cero). Veamos cómo quedaría, y fijémonos en la célula del centro. Por un momento olvidémonos de células muertas y vivas:
In [1]:
from IPython.display import HTML
from jinja2 import Template
def make_table(M, N, id):
table = Template("""<table id="{{ id }}">
{% for row in rows %}
<tr>
{% for col in cols %}<td id="c{{ row }}{{ col }}">$c_{ {{ row }}{{ col }} }$</td>{% endfor %}
</tr>
{% endfor %}
</table>""").render(rows=range(M), cols=range(N), id=id)
return table
style = """<style>
table#t1 #c22 { background-color: #bcd; border-width: 2px; }
table#t1 #c11,
table#t1 #c12,
table#t1 #c13,
table#t1 #c21,
table#t1 #c23,
table#t1 #c31,
table#t1 #c32,
table#t1 #c33 { background-color: #def; }
</style>
"""
HTML(style + make_table(5, 5, "t1"))
Out[1]:
Todas las células tienen 8 vecinos, también las de los bordes (imaginad el tablero como un nivel de Asteroids). Si ahora las células pueden estar vivas o muertas, el número de vecinos vivos de $c_{22}$ sería:
$$ v_{22} = c_{11} + c_{12} + c_{13} + c_{21} + c_{23} + c_{31} + c_{32} + c_{33} $$donde $c_{ij}$ puede valer 1 o 0, dependiendo de si está viva o no, y el número $v_{22}$ es el número de vecinos vivos de $c_{22}$. También podemos escribir esto así, aunque ahora parezca tonto:
$$ v_{22} = 1 · c_{11} + 1 · c_{12} + 1 · c_{13} + 1 · c_{21} + 1 · c_{23} + 1 · c_{31} + 1 · c_{32} + 1 · c_{33} $$o, incluyendo todas las células del tablero:
$$ v_{22} = 0 · c_{00} + 0 · c_{01} +{}...{}+ 0 · c_{10} + 1 · c_{11} + 1 · c_{12} +{}...{}+ 0 · c_{44} $$Por ejemplo, si la única célula viva del tablero fuese $c_{12}$, obviamente esta operación resulta 1:
$$ v_{22} = 0 · 0 + 0 · 0 +{}...{}+ 0 · 0 + 1 · 0 + 1 · 1 +{}+ 0 · 0 = 1 $$Vemos que podemos expresar esto como un producto escalar:
$$ v_{22} = \begin{pmatrix}0 & 0 &{}...{}& 0 & 1 & 1 &{}...{}& 0\end{pmatrix} · \begin{pmatrix}c_{00} \\ c_{01} \\ \vdots \\ c_{10} \\ c_{11} \\ c_{12} \\ \vdots \\ c_{44}\end{pmatrix} = R_{22} · T$$siendo $R_{22}$ el vector de unos y ceros que indica cuáles células cuentan en la suma y cuáles no y $T$ un vector columna donde tenemos todas las células de tablero. Esto se pone interesante ;)
Veamos cómo es el vector $R_{22}$ completo:
$$ R_{22} = \small \begin{pmatrix}0 & 0 & 0 & 0 & 0 & | & 0 & 1 & 1 & 1 & 0 & | & 0 & 1 & 0 & 1 & 0 & | & 0 & 1 & 1 & 1 & 0 & | & 0 & 0 & 0 & 0 & 0\end{pmatrix}$$Esto es importante. He escrito un separador cada cinco elementos para que te imagines la situación por filas. Los primeros cinco son todo ceros, porque ninguna de las celdas de la primera fila es vecina de la celda $v_{22}$. La segunda fila tiene los tres del centro, la tercera fila tiene dos (¡la celda $v_{22}$ no es vecina de sí misma!), la cuarta vuelve a tener las tres del centro y la quinta de nuevo todo ceros.
La pregunta es, ¿qué pinta tienen el resto de vectores $R$? ¿Se parecen entre sí? Veamos otro ejemplo:
In [2]:
style = """<style>
table#t2 #c21 { background-color: #bcd; border-width: 2px; }
table#t2 #c10,
table#t2 #c11,
table#t2 #c12,
table#t2 #c20,
table#t2 #c22,
table#t2 #c30,
table#t2 #c31,
table#t2 #c32 { background-color: #def; }
</style>
"""
HTML(style + make_table(5, 5, "t2"))
Out[2]:
¿Ven bien mis ojos? ¿Es el mismo vector de antes, desplazado una posición a la izquierda? ¡Sí! :D Si repetimos este proceso con todas las celdas y apilamos todos esos vectores por filas, obtendremos una matriz que visualizada en matplotlib quedaría así:
In [3]:
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import cm
A55 = np.array([
[0, 1, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 1],
[1, 0, 1, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0],
[0, 1, 0, 1, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0],
[0, 0, 1, 0, 1, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1],
[1, 0, 0, 1, 0, 1, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1],
[1, 1, 0, 0, 1, 0, 1, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1],
[1, 1, 1, 0, 0, 1, 0, 1, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 1, 1, 1, 0, 0, 1, 0, 1, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 1, 1, 1, 0, 0, 1, 0, 1, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 1, 1, 1, 0, 0, 1, 0, 1, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 1, 1, 1, 0, 0, 1, 0, 1, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 1, 0, 1, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 1, 0, 1, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 1, 0, 1, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 1, 0, 1, 0, 0, 1, 1, 1, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 1, 0, 1, 0, 0, 1, 1, 1, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 1, 0, 1, 0, 0, 1, 1, 1, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 1, 0, 1, 0, 0, 1, 1, 1, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 1, 0, 1, 0, 0, 1, 1, 1],
[1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 1, 0, 1, 0, 0, 1, 1],
[1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 1, 0, 1, 0, 0, 1],
[1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 1, 0, 1, 0, 0],
[0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 1, 0, 1, 0],
[0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 1, 0, 1],
[1, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 1, 0]
])
"""
$$ A = \tiny \begin{pmatrix}
0 & 1 & 0 & 0 & 1 & 1 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 1 & 1 & 0 & 0 & 1 \\
1 & 0 & 1 & 0 & 0 & 1 & 1 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 1 & 1 & 0 & 0 \\
0 & 1 & 0 & 1 & 0 & 0 & 1 & 1 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 1 & 1 & 0 \\
0 & 0 & 1 & 0 & 1 & 0 & 0 & 1 & 1 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 1 & 1 \\
1 & 0 & 0 & 1 & 0 & 1 & 0 & 0 & 1 & 1 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 1 \\
1 & 1 & 0 & 0 & 1 & 0 & 1 & 0 & 0 & 1 & 1 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 \\
1 & 1 & 1 & 0 & 0 & 1 & 0 & 1 & 0 & 0 & 1 & 1 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
0 & 1 & 1 & 1 & 0 & 0 & 1 & 0 & 1 & 0 & 0 & 1 & 1 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
0 & 0 & 1 & 1 & 1 & 0 & 0 & 1 & 0 & 1 & 0 & 0 & 1 & 1 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
0 & 0 & 0 & 1 & 1 & 1 & 0 & 0 & 1 & 0 & 1 & 0 & 0 & 1 & 1 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 & 1 & 1 & 1 & 0 & 0 & 1 & 0 & 1 & 0 & 0 & 1 & 1 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 & 0 & 1 & 1 & 1 & 0 & 0 & 1 & 0 & 1 & 0 & 0 & 1 & 1 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 & 0 & 0 & 1 & 1 & 1 & 0 & 0 & 1 & \normalsize{0} & 1 & 0 & 0 & 1 & 1 & 1 & 0 & 0 & 0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 1 & 1 & 0 & 0 & 1 & 0 & 1 & 0 & 0 & 1 & 1 & 1 & 0 & 0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 1 & 1 & 0 & 0 & 1 & 0 & 1 & 0 & 0 & 1 & 1 & 1 & 0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 1 & 1 & 0 & 0 & 1 & 0 & 1 & 0 & 0 & 1 & 1 & 1 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 1 & 1 & 0 & 0 & 1 & 0 & 1 & 0 & 0 & 1 & 1 & 1 & 0 & 0 \\
0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 1 & 1 & 0 & 0 & 1 & 0 & 1 & 0 & 0 & 1 & 1 & 1 & 0 \\
0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 1 & 1 & 0 & 0 & 1 & 0 & 1 & 0 & 0 & 1 & 1 & 1 \\
1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 1 & 1 & 0 & 0 & 1 & 0 & 1 & 0 & 0 & 1 & 1 \\
1 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 1 & 1 & 0 & 0 & 1 & 0 & 1 & 0 & 0 & 1 \\
1 & 1 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 1 & 1 & 0 & 0 & 1 & 0 & 1 & 0 & 0 \\
0 & 1 & 1 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 1 & 1 & 0 & 0 & 1 & 0 & 1 & 0 \\
0 & 0 & 1 & 1 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 1 & 1 & 0 & 0 & 1 & 0 & 1 \\
1 & 0 & 0 & 1 & 1 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 1 & 1 & 0 & 0 & 1 & 0
\end{pmatrix} $$
"""
plt.matshow(A55, cmap=cm.gray_r)
plt.grid(False)
Lo especial de esta matriz es que si la multiplicamos por el vector $T$, nos da un vector $V$ que cuenta cuántas vecinas vivas tiene cada celda.
$$ V = A T $$¡Ya tenemos el misterio resuelto! Podemos calcular cada paso en el Juego de la vida con una simple multiplicación matricial. Esta matriz no es más que la matriz de adyacencia de un grafo tal que cada vértice es una célula del tablero y cada arista indica dos células adyacentes. Hemos pasado de autómatas celulares a teoría de grafos, casi sin despeinarnos :)
Sin embargo, no vamos a utilizar el grafo del tablero porque nos complicaríamos en exceso. En su lugar, vamos a construir la matriz directamente usando matrices dispersas.
Ya tenemos un artículo en Pybonacci sobre cómo construir matrices tridiagonales en Python. Visto por la gráfica anterior que la matriz de adyacencia tiene una estructura diagonal, vamos a hacer algo muy similar a lo que ya hicimos entonces: solo tenemos que crear el vector $R$, y después vamos a crear la matriz por diagonales utilizando el módulo scipy.sparse.
In [4]:
M, N = 5, 5 # Dimensiones del tablero
R_22 = np.zeros((5, 5)) # Idea de Chema Cortés
R_22[1:4, 1:4] = np.ones((3, 3))
R_22[2, 2] = 0
print(R_22)
R_22 = R_22.flatten()
Nos interesa $R_{00}$:
In [5]:
R_00 = np.roll(R_22, 2 * (-N - 1)) # Trasladamos el vector
R_00
Out[5]:
Y ahora viene el truco: como sólo nos interesa almacenar los unos, hay que decir a la función spdiags cuántos hay y en qué posiciones se encuentran.
In [6]:
cnt = np.count_nonzero(R_00)
nonzero_idx = R_00.nonzero()[0]
print(cnt, nonzero_idx)
In [7]:
import scipy.sparse as ss
np.set_printoptions(linewidth=100)
ss.diags(np.ones(cnt), # Número de unos
nonzero_idx, # Posiciones
shape=(M * N, M * N), # Tamaño de la matriz
format="csr", # Formato (CSR es mejor para productos matriz-vector)
dtype=int # Tipo entero
).todense()
Out[7]:
El problema es que los índices positivos solo representan el triángulo superior de la matriz. Para rellenar también el triángulo inferior podemos escribir esto:
In [8]:
ss.diags(np.ones(cnt * 2), # Número de unos
list(-nonzero_idx) + list(nonzero_idx), # Posiciones
shape=(M * N, M * N), # Tamaño de la matriz
format="csr", # Formato (CSR es mejor para productos matriz-vector)
dtype=int # Tipo entero
).todense()
Out[8]:
¡Es el mismo resultado que hemos obtenido antes! Si generalizamos lo que hemos escrito a una función:
In [9]:
def adjacency_matrix(M, N):
"""Matriz de adyacencia de un tablero M x N.
Reglas:
* 8 vecinos por celda
* Geometría toroidal
No he implementado el caso N < 3 o M < 3.
"""
if M < 3 or N < 3:
raise NotImplementedError
R_11 = np.zeros((M, N))
R_11[0:3, 0:3] = np.ones((3, 3))
R_11[1, 1] = 0
R_00 = np.roll(R_11.flatten(), -N - 1)
mn = M * N
cnt = np.count_nonzero(R_00)
nonzero_idx = R_00.nonzero()[0]
A = ss.diags(np.ones(cnt * 2, dtype=int), list(-nonzero_idx) + list(nonzero_idx),
shape=(mn, mn), format="csr", dtype=int)
return A
In [10]:
def plot_board(A, axis=True):
plt.matshow(A, cmap=cm.gray_r)
plt.grid(False)
if not axis:
plt.axis('off')
plot_board(adjacency_matrix(5, 5).todense())
plot_board(adjacency_matrix(4, 12).todense())
Vemos que calcular la matriz de adyacencia es un tanto costoso para dimensiones grandes:
In [11]:
%timeit adjacency_matrix(100, 200)
Pero no importa porque este cálculo sólo lo vamos a hacer una vez en toda la simulación, puesto que las dimensiones del tablero no cambian. Podemos cachear nuestros resultados:
In [12]:
adj_matrices = {}
def get_adj_matrix(M, N):
try:
A = adj_matrices[(M, N)]
except KeyError:
A = adjacency_matrix(M, N)
adj_matrices[(M, N)] = A
return A
In [13]:
%timeit -n1 -r1 get_adj_matrix(100, 200)
In [14]:
%timeit get_adj_matrix(100, 200)
Y ahora, los resultados :)
Jake VanderPlas proponía estos dos métodos para calcular el paso:
In [15]:
from scipy.signal import convolve2d
def life_step_1(X):
"""Game of life step using generator expressions"""
nbrs_count = sum(np.roll(np.roll(X, i, 0), j, 1)
for i in (-1, 0, 1) for j in (-1, 0, 1)
if (i != 0 or j != 0))
return (nbrs_count == 3) | (X & (nbrs_count == 2))
def life_step_2(X):
"""Game of life step using scipy tools"""
nbrs_count = convolve2d(X, np.ones((3, 3)), mode='same', boundary='wrap') - X
return (nbrs_count == 3) | (X & (nbrs_count == 2))
Veamos cómo funcionan frente al nuestro:
In [16]:
def life_step_3(X):
"""Calcula el paso usando la matriz de adyacencia."""
A = get_adj_matrix(*X.shape)
nbrs_count = A.dot(X.flatten()).reshape(X.shape)
return (nbrs_count == 3) | (X & (nbrs_count == 2))
In [17]:
tab = np.round(np.random.rand(100, 200)).astype(int) # Tablero aleatorio
get_adj_matrix(*tab.shape) # Cacheamos la matriz correspondiente
Out[17]:
In [18]:
%timeit life_step_1(tab)
In [19]:
%timeit life_step_2(tab)
In [20]:
%timeit life_step_3(tab)
Vencedor: matriz de adyacencia :D ¡Cinco veces más rápido! (En mi ordenador)
Y para terminar:
In [21]:
from functools import reduce
for ii in range(0, 18, 6):
plot_board(reduce(lambda x, _: life_step_3(x), range(ii), tab), False)
plt.title("{} iteraciones".format(ii))
¡Un saludo! ;)