用python演绎神奇的生命游戏，在游戏中学习numpy和matplotlib动画

Posted 2022-08-05 天元浪子

在一个二维网格中，假定每一个方格代表一个细胞，每个细胞有存活和死亡两种状态，其初始生存状态随机确定。每隔一段时间检查一次细胞的生存状态，每个细胞的生存状态由其周围的8个细胞的生存状态决定，具体规则如下：

1. 如果一个细胞周围的活细胞数量超过3个或少于2个，则该细胞死亡；
2. 如果一个细胞周围的活细胞数量等于2个，则该细胞生存状态不变；
3. 如果一个细胞周围的活细胞数量等于3个，则该细胞复活或继续存活；

如果用黑色表示死亡的细胞，用白色表示存活的细胞，这个二维网格最初看起来是杂乱无章的，但是经过一段时间的演化之后，就会产生稳定而规律的变化，甚至是持续的、有规律的移动，类似物质交换或生命运动。

这就是英国数学家约翰·何顿·康威（John Horton Conway）于1970年发明的生命游戏，也被称作康威生命游戏。上面这几张图是生命游戏中的经典图案，转自康威生命游戏官方网站。仅凭3条细胞繁衍和死亡的简单规则，生命游戏就可以在计算机上模拟出丰富的生命演化过程，甚至可以模拟出与真实生命相当的复杂度——只要计算机内存足够大、计算能力足够强。

数学家康威证明了生命游戏具有图灵完备性，允许在生命游戏中模拟任何其他生命游戏规则。康威生命游戏是人工生命的经典研究，推动了元胞自动机（Cellular Automaton）理论的发展。元胞自动机作为一种仿真算法在近两年的数学建模竞赛中经常出现，可谓数学建模竞赛的万金油。康威生命游戏就是一种在二维网格上定义的元胞自动机。

如果用python演绎这个游戏的话，似乎不难。比如， 可以用一个二维列表来模拟游戏中的二维表格，用两层嵌套的循环结构来检查每一个细胞的生存状态。不过，这不是一个好的主意：如果二维表格稍微大一点，这样的代码会慢到无法忍受。下面这段代码使用numpy数组来模拟游戏中的二维表格，不需要循环就可以完成一次细胞演化，速度直追c语言。用空间换时间，这是应用numpy时避免显式循环的最常用的手段之一。

import numpy as np

def evolve(bio):
    """细胞演化，bio是用来模拟二维表格的二维数组"""
    
    # bio右下角元素 + bio尾行 + bio左下角元素，水平堆叠为一行
    top = np.hstack((bio[-1,-1], bio[-1], bio[-1,0])) 
    
    # bio右上角元素 + bio首行 + bio左上角元素，水平堆叠为一行
    bottom = np.hstack((bio[0,-1], bio[0], bio[0,0])) 
    
    # bio尾列 + bio + bio首列，水平堆叠为二维数组
    center = np.hstack((bio[:,-1:], bio, bio[:,0:1])) 
    
    # top + center + bottom，垂直堆叠为二维数组u，u比bio“胖了一圈”
    u = np.vstack((top, center, bottom))  
    
    s1 = u[:-2,:-2]     # 每个元素左上的元素构成的二维数组
    s2 = u[:-2,1:-1]    # 每个元素上方的元素构成的二维数组
    s3 = u[:-2,2:]      # 每个元素右上的元素构成的二维数组
    s4 = u[1:-1,:-2]    # 每个元素左侧的元素构成的二维数组 
    s5 = u[1:-1,2:]     # 每个元素右侧的元素构成的二维数组
    s6 = u[2:,:-2]      # 每个元素左下的元素构成的二维数组
    s7 = u[2:,1:-1]     # 每个元素下方的元素构成的二维数组
    s8 = u[2:,2:]       # 每个元素右下的元素构成的二维数组
    
    # s是和bio同样大小的二维数组，记录了bio每个元素周围值为1的元素个数
    s = s1 + s2 + s3 + s4 + s5 + s6 + s7 + s8 
    
    # bio中每个元素：如果周边元素为1的个数大于3或小于2，则该元素为0
    bio[np.where((s > 3) | (s < 2))] = 0
    
    # bio中每个元素：如果周边元素为1的个数等于3，则该元素为1
    bio[np.where(s == 3)] = 1
    
    return bio

接下来测试细胞演化函数evolve是否符合设计预期。

rows, cols = 8, 8 # 用8行8列的表格测试
n = int(0.5 * rows * cols) # 初始状态50%的细胞是活的
bio = np.zeros((rows, cols), dtype=np.uint8) # rows行cols列全0数组，模拟生命空间
bio[(np.random.randint(0, rows, n), np.random.randint(0, cols, n))] = 1 # 随机填写n个1

print('初始状态：')
print(bio)

bio = evolve(bio)
print('演化一次：')
print(bio)

测试结果看起来是正确的。

初始状态：
[[0 1 1 0 0 0 0 0]
 [1 0 0 1 1 0 0 0]
 [0 1 1 0 1 1 1 1]
 [1 1 0 0 1 0 0 1]
 [1 0 0 1 0 0 1 1]
 [1 0 1 0 0 0 1 0]
 [0 0 0 0 0 1 0 0]
 [1 0 1 0 0 0 0 1]]
演化一次：
[[0 0 1 0 0 0 0 1]
 [1 0 0 0 1 0 1 1]
 [0 0 1 0 0 0 1 0]
 [0 0 0 0 1 0 0 0]
 [0 0 1 1 0 1 1 0]
 [1 1 0 0 0 1 1 0]
 [1 0 0 0 0 0 1 0]
 [1 0 1 0 0 0 0 0]]

请按任意键继续. . .

怎么能够直观地看到细胞连续演化的过程呢？matplotlib.animation提供了一个便捷的手段，这就是FuncAnimation函数。FuncAnimation函数用起来稍显繁琐，不过下面的代码展示了一种逻辑清晰、通俗易懂的用法，可以用来生成复杂的动画，唯一的缺点就是刷新频率较低，且难以提升。

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation

def evolve():
    """细胞演化，bio是用来模拟二维表格的二维数组"""
        
    top = np.hstack((bio[-1,-1], bio[-1], bio[-1,0]))
    bottom = np.hstack((bio[0,-1], bio[0], bio[0,0])) 
    center = np.hstack((bio[:,-1:], bio, bio[:,0:1])) 
    
    u = np.vstack((top, center, bottom))  
    s = u[:-2,:-2] + u[:-2,1:-1] + u[:-2,2:]  + u[1:-1,:-2] + u[1:-1,2:] + u[2:,:-2] + u[2:,1:-1] + u[2:,2:] 
    
    bio[np.where((s > 3) | (s < 2))] = 0
    bio[np.where(s == 3)] = 1

def animate(frame):
    """动画函数"""

    plt.cla() # 清空画布
    plt.imshow(bio, alpha=0.8) # 增加参数cmap='gray'，显示黑白效果
    evolve() # 细胞演化

rows, cols = 16, 16 # 16行16列
n = int(0.5 * rows * cols) # 初始状态50%的细胞是活的
bio = np.zeros((rows, cols), dtype=np.uint8) # rows行cols列全0数组，模拟生命空间
bio[(np.random.randint(0, rows, n), np.random.randint(0, cols, n))] = 1 # 随机填写n个1

plt.figure(figsize=(8, 8)) # 设置画布
anim = FuncAnimation(plt.gcf(), animate) # 参数1是画布，参数2是动画函数，默认刷新周期200毫秒
plt.show() # 显示画布

若想直接生成gif文件或者mp4文件，只需要将上面代码的最后两行做如下改动。

anim = FuncAnimation(plt.gcf(), animate, frames=100) # frames用来设置生成的动画或视频帧数
anim.save(r'd:\\bio.gif', fps=100)
#anim.save(r'd:\\bio.mp4', writer='ffmpeg', fps=100)

这是在16行16列的二维空间中进行100次演化的演示动画。