A raíz de que Jake VanderPlas publicara el otro día su propia aproximación al juego de la vida de Conway, mucho mejor que aquella que publicamos en Pybonacci, se vio que nuestro artículo necesitaba una revisión importante. Chema Cortés además me dio una idea en los comentarios de un método para calcular la matriz vecindario y aquí voy a explorar si la mejora es significativa (con respecto a los otros dos métodos que dio JvdP).
Vamos a utilizar nosotros también una geometría toroidal.
import numpy as np
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"""
from scipy.signal import convolve2d
nbrs_count = convolve2d(X, np.ones((3, 3)), mode='same', boundary='wrap') - X
return (nbrs_count == 3) | (X & (nbrs_count == 2))
Vamos a construir un tablero aleatorio y a comparar los dos métodos.
tab = np.round(np.random.rand(60, 60)).astype(int)
%timeit life_step_1(tab)
1000 loops, best of 3: 451 µs per loop
%timeit life_step_2(tab)
1000 loops, best of 3: 541 µs per loop
Se ve que el segundo es claramente más lento. Veamos ahora cómo funciona el mío.
Si consideramos el tablero del Juego de la vida 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
$$ A v_T = v_V $$siendo $v_T$ el vector columna resultante de aplanar la matriz $T$.
Por ejemplo, para este tablero 3×3:
$$ T = \begin{pmatrix} 1 & 1 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \end{pmatrix} $$%matplotlib inline
import matplotlib.pyplot as plt
from matplotlib import cm
T = np.zeros((3, 3), dtype=int)
T[0, 0:2] = 1
plt.matshow(T, cmap=cm.gray_r)
plt.grid(False)
La matriz de adyacencia (dada solo por la geometría) será:
$ A = \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} $
El vector será:
$ v_T = \begin{pmatrix} 1 \\ 1 \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ \end{pmatrix} $
y por tanto el vecindario:
$ v_V = A v_T = \begin{pmatrix} 1 \\ 1 \\ 2 \\ 2 \\ 2 \\ 2 \\ 2 \\ 2 \\ 2 \\ \end{pmatrix} $
y la matriz correspondiente:
$ V = \begin{pmatrix} 1 & 1 & 2 \\ 2 & 2 & 2 \\ 2 & 2 & 2 \\ \end{pmatrix} $
Tomemos la función para calcular el vecindario de game_step_1
:
def nbrs_1(X):
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
nbrs_1(T)
array([[1, 1, 2], [2, 2, 2], [2, 2, 2]])
Nuestro ejemplo funciona.
La ventaja de este método es que la matriz de adyacencia se construye una sola vez, cuando conocemos las dimensiones del tablero y las condiciones de contorno, y calcular el paso siguiente solo es un producto matriz por vector. Pero ¿cómo construimos la matriz de adyacencia? Seguro que hay métodos más eficientes, pero este es el que he utilizado yo:
def adjacency(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
block_ij = np.ones((3, 3), dtype=int)
block_ij[1, 1] = 0
row_11 = np.vstack([block_ij] + [np.zeros(3, dtype=int)] * (M - 3))
row_11 = np.column_stack([row_11] + [np.zeros(M, dtype=int)] * (N - 3)).flatten()
row_00 = np.roll(row_11, -N - 1)
A = np.vstack(np.roll(row_00, ii) for ii in range(M * N))
return A
from numpy.testing import assert_raises
def test_adjacency_fails_small_board():
assert_raises(NotImplementedError, adjacency, 3, 2)
assert_raises(NotImplementedError, adjacency, 2, 3)
assert_raises(NotImplementedError, adjacency, 2, 2)
def test_adjacency_shape():
M = 50
N = 30
assert adjacency(M, N).shape == (M * N, M * N)
def test_adjacency_int():
M = 50
N = 30
assert adjacency(M, N).dtype == np.dtype(int)
def test_adjacency_diag_zero():
M = 4
N = 3
assert np.all(np.diag(adjacency(M, N)) == 0)
def test_adjacency_symmetric():
M = 42
N = 31
assert np.all(adjacency(M, N).T == adjacency(M, N))
test_adjacency_fails_small_board()
test_adjacency_shape()
test_adjacency_int()
test_adjacency_diag_zero()
test_adjacency_symmetric()
%timeit adjacency(*tab.shape)
1 loops, best of 3: 210 ms per loop
Vemos que construir la matriz de adyacencia es bastante lento. Sin embargo, como hemos dicho, solo vamos a construirla una vez: podemos cachearla.
adjacency_matrices = {}
def get_adj_matrix(M, N):
try:
A = adjacency_matrices[(M, N)]
except KeyError:
A = adjacency(M, N)
adjacency_matrices[(M, N)] = A
return A
%timeit -n1 -r1 get_adj_matrix(*tab.shape) # Sin cachear
1 loops, best of 1: 234 ms per loop
%timeit get_adj_matrix(*tab.shape) # Cacheado
1000000 loops, best of 3: 530 ns per loop
Un millón de veces más rápido. Probemos este método:
def life_step_3(X):
"""Game of life step using adjacency matrices"""
A = get_adj_matrix(*X.shape)
nbrs_count = np.dot(A, X.flatten()).reshape(X.shape)
return (nbrs_count == 3) | (X & (nbrs_count == 2))
%timeit life_step_3(tab)
100 loops, best of 3: 13.7 ms per loop
¡Mal! El método resulta ser 100 veces más lento que los anteriores. Además, vemos que hay problemas cuando las dimensiones de la matriz empiezan a crecer:
# ¡Cuidado! Esto puede congelar tu ordenador unos instantes
#A = adjacency(90, 90)
#print(A.nbytes / 1024. / 1024.) # Megabytes
Pero en realidad, casi todos los elementos de la matriz son nulos:
def sparsity(A):
return np.count_nonzero(A) / A.size
print("Nonzero elements: {:.2f} %".format(sparsity(get_adj_matrix(*tab.shape)) * 100))
plt.figure(1, figsize=(12, 12))
plt.matshow(get_adj_matrix(*tab.shape), fignum=1, cmap=cm.gray_r)
plt.grid(False)
Nonzero elements: 0.22 %
Lo cual sugiere que podemos usar matrices dispersas, construidas a partir de las diagonales. Solo tenemos que reescribir la función adjacency
utilizando scipy.sparse
.
from scipy import sparse as ss
def adjacency2(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
block_ij = np.ones((3, 3), dtype=int)
block_ij[1, 1] = 0
row_11 = np.pad(block_ij, ((0, M - 3), (0, N - 3)), mode='constant').flatten()
row_00 = np.roll(row_11, -N - 1)
mn = M * N
A = ss.diags(np.concatenate([row_00, row_00[1:][::-1]]), np.concatenate([np.arange(mn), -np.arange(mn)[1:]]),
shape=(mn, mn), dtype=int)
return A
adjacency2(*tab.shape)
<3600x3600 sparse matrix of type '<class 'numpy.int64'>' with 12960000 stored elements (7199 diagonals) in DIAgonal format>
Todavía podemos hacerlo mucho mejor. Aunque hemos ganado algo, estamos guardando una cantidad enorme de ceros innecesarios en la matriz. Tenemos que idear un método más inteligente de construirla.
row_00 = np.ones(9, dtype=int)
row_00[0] = 0
row_00
array([0, 1, 1, 1, 1, 1, 1, 1, 1])
np.nonzero(row_00)[0]
array([1, 2, 3, 4, 5, 6, 7, 8])
M, N = tab.shape
block_ij = np.ones((3, 3), dtype=int)
block_ij[1, 1] = 0
row_11 = np.vstack([block_ij] + [np.zeros(3, dtype=int)] * (M - 3))
row_11 = np.column_stack([row_11] + [np.zeros(M, dtype=int)] * (N - 3)).flatten()
row_00 = np.roll(row_11, -N - 1)
row_00
array([0, 1, 0, ..., 0, 0, 1])
list(-row_00.nonzero()[0]) + list(row_00.nonzero()[0])
[-1, -59, -60, -61, -3539, -3540, -3541, -3599, 1, 59, 60, 61, 3539, 3540, 3541, 3599]
def adjacency3(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
block_ij = np.ones((3, 3), dtype=int)
block_ij[1, 1] = 0
row_11 = np.pad(block_ij, ((0, M - 3), (0, N - 3)), mode='constant').flatten()
row_00 = np.roll(row_11, -N - 1)
mn = M * N
cnt = np.count_nonzero(row_00)
nonzero_idx = row_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
adjacency2(100, 20)
<2000x2000 sparse matrix of type '<class 'numpy.int64'>' with 4000000 stored elements (3999 diagonals) in DIAgonal format>
adjacency3(100, 20)
<2000x2000 sparse matrix of type '<class 'numpy.int64'>' with 16000 stored elements in Compressed Sparse Row format>
Ahora parece que sí :)
def test_adjacency3_fails_small_board():
assert_raises(NotImplementedError, adjacency3, 3, 2)
assert_raises(NotImplementedError, adjacency3, 2, 3)
assert_raises(NotImplementedError, adjacency3, 2, 2)
def test_adjacency3_shape():
M = 50
N = 30
assert adjacency3(M, N).shape == (M * N, M * N)
def test_adjacency3_int():
M = 50
N = 30
assert adjacency3(M, N).dtype == np.dtype(int)
def test_adjacency3_diag_zero():
M = 4
N = 3
assert np.all(np.diag(adjacency3(M, N).todense()) == 0)
def test_adjacency3_symmetric():
M = 42
N = 31
assert np.all(adjacency3(M, N).todense().T == adjacency3(M, N).todense())
test_adjacency3_fails_small_board()
test_adjacency3_shape()
test_adjacency3_int()
test_adjacency3_diag_zero()
test_adjacency3_symmetric()
%timeit adjacency(*tab.shape)
1 loops, best of 3: 186 ms per loop
%timeit adjacency3(*tab.shape)
100 loops, best of 3: 6.25 ms per loop
adjacency3_matrices = {}
def get_adj3_matrix(M, N):
try:
A = adjacency3_matrices[(M, N)]
except KeyError:
A = adjacency3(M, N)
adjacency3_matrices[(M, N)] = A
return A
%timeit -n1 -r1 get_adj3_matrix(*tab.shape)
1 loops, best of 1: 22.9 ms per loop
%timeit get_adj3_matrix(*tab.shape)
1000000 loops, best of 3: 507 ns per loop
def life_step_33(X):
"""Game of life step using adjacency matrices"""
A = get_adj3_matrix(*X.shape)
nbrs_count = A.dot(X.flatten()).reshape(X.shape)
return (nbrs_count == 3) | (X & (nbrs_count == 2))
%timeit life_step_1(tab)
1000 loops, best of 3: 437 µs per loop
%timeit life_step_33(tab)
10000 loops, best of 3: 92.6 µs per loop
¡Éxito! El nuevo método es cinco veces más rápido que los anteriores.
from functools import reduce
from time import sleep
def random_board(M, N):
return np.round(np.random.rand(M, N)).astype(int)
def plot_board(bb, num=None, save=False):
plt.matshow(bb, cmap=cm.gray_r)
plt.grid(False)
if save:
plt.savefig("bb{:03d}.png".format(num or 0))
bb = random_board(48, 92)
for ii in range(0, 18, 6):
plot_board(reduce(lambda x, _: life_step_33(x), range(ii), bb), ii)