生命游戏中的经典涌现“滑翔机”
Coway的生命游戏在先前的文章中介绍过，是一种元胞自动机模型。把一个具有类似于棋盘网格的平面看成是一个“世界”，其中每一个格子都看成是一个“生命”。这个生命基于一些简单的法则改变自己生存或死亡的状态。之前已经说过，这个看似“简陋的”世界，已经被MIT的数学家Wolfram证明其计算能力与通用计算机相同。简单来说，就是给定合适的规则和初值，你可以得到任何可计算的结果。今天介绍用Python和C++两种语言共四种不同的实现。在了解生命模型的同时，对比一下两种语言的性能。由于时间和设备等条件的限制，本文的对比仅作为一个粗略直观的参考。
下面给出要实现的生命游戏规则：
1、一个生命为生存状态时，若其周围8个邻居有且仅有2或者3个活着，则其继续存活；
2、一个生命为生存状态时，若其周围8个邻居中活着的数量超过3个，则拥挤而死；
3、一个生命为生存状态时，若其周围8个邻居中活着的数量少于2个，则寂寞而死；
4、一个生命为死亡状态时，若其周围8个邻居中或者数量仅有3个，则获得重生。
为了展示动态效果，我们这里使用的Python中的pygame游戏引擎，而为了测试其性能上限，利用一个像素来代表一个生命。因而，接下来代码中所有的操作都是针对屏幕像素的操作。Talk is cheap, show me the code.话不多说，直接上代码和结果。
import time
import pygame, sys
from pygame.locals import *
import random
HEIGHT = 400
WIDTH =400
WHITE = (255, 255, 255)
BLACK = (0,0,0)
FPS = 60
pygame.init()
fpsClock = pygame.time.Clock()
screen = pygame.display.set_mode((HEIGHT, WIDTH))
screen.fill(BLACK)
pxarray_orig = pygame.PixelArray(screen.copy())
for i in range(1,WIDTH-1):
for j in range(1,HEIGHT-1):
if random.randrange(2) < 1:
screen.fill(0,((i,j),(1,1)))
else:
screen.fill(16777215,((i,j),(1,1)))
screen_new = screen.copy()
pygame.display.update()
x = 0
startTime = time.time()
while x<100: # main game loop
pxarray_orig = pygame.PixelArray(screen.copy())
init_color_value = 16777215
for i in range(1,WIDTH-1):
for j in range(1,HEIGHT-1):
my_state = pxarray_orig[i,j]
neib_position = [(i,j+1),(i,j-1),(i+1,j),(i-1,j),(i+1,j+1),(i-1,j-1),(i+1,j-1),(i-1,j+1)]
neighbor_state = 0
for position in neib_position:
neighbor_state += pxarray_orig[position[0],position[1]]
if my_state == 16777215 and (neighbor_state == init_color_value*2 or neighbor_state==init_color_value*3):
screen_new.fill(16777215,((i,j),(1,1)))
elif my_state == 16777215 and neighbor_state > init_color_value *3:
screen_new.fill(0,((i,j),(1,1)))
elif my_state == 16777215 and neighbor_state < init_color_value*2:
screen_new.fill(0,((i,j),(1,1)))
elif my_state == 0 and neighbor_state==init_color_value*3:
screen_new.fill(16777215,((i,j),(1,1)))
screen.blit(screen_new,screen.get_rect().topleft)
pygame.display.flip()
for event in pygame.event.get():
if event.type == QUIT:
pygame.quit()
sys.exit()
#fpsClock.tick(FPS)
x += 1
executionTime = (time.time() - startTime)
print('Execution time in seconds: ' + str(executionTime))
pygame.quit()
输出结果为：
以上代码设定了400*400=16000个像素，即16000个生命的世界。每一轮模拟中，系统都要询问16000个生命的状态及其周围8个邻居的状态，设定模拟100轮后停止，并输出所用时间。可以看到画面的流畅度非常差，感觉像在放ppt。
众所周知，Python是一门解释的语言，在运行过程中，其解释器会一边对变量类型进行自动判定，而在嵌套的for循环中，判定次数陡增，进而造成性能急速下降。因此，Pyhon的确不适合执行大规模的仿真。但需要了解是，Python的许多数据分析库是用C语言写的，例如numpy，因此其数据分析的性能依然表现不俗。C语言的执行效率比Python快得多。作为一门接近机器而非人类的语言，其编写的复杂程度也比Python高很多。为了综合Python的“人性”和C语言的高效率（机器运行的效率），Cython横空出世。Cython是Python的一个库，底层是由C语言写的，并且有自己的编译器。其最大的好处就是可以在Python中无缝添加Cython代码，其编译器会将Python代码翻译成C代码，进行高效率执行。事实上，纯Python的代码在经过Cython的翻译加工后也能提升不少性能。在Jupyter Notebook中仅需要在上面代码顶端加上“%%cython”这样一个简单的操作就能提升20%的速度。
为了将进一步利用Cython的优势，我们将其中可以定义出明确类型的变量全部显示地给出类型，如下面的代码所示。
%%cython
import time
import pygame, sys
from pygame.locals import *
import random
import numpy
cdef int HEIGHT = 400
cdef int WIDTH =400
WHITE = (255, 255, 255)
BLACK = (0,0,0)
cdef int FPS = 100
#cdef numpy.ndarray pxarray_orig
pygame.init()
fpsClock = pygame.time.Clock()
screen = pygame.display.set_mode((HEIGHT, WIDTH))
screen.fill(BLACK)
pxarray_orig = numpy.array(pygame.PixelArray(screen.copy()))
pxarray_orig = pygame.PixelArray(screen.copy())
cdef int i
cdef int j
cdef my_state
cdef init_color_value
cdef list neib_position
cdef int neighbor_state
for i in range(1,WIDTH-1):
for j in range(1,HEIGHT-1):
if random.randrange(2) < 1:
screen.fill(0,((i,j),(1,1)))
else:
screen.fill(16777215,((i,j),(1,1)))
screen_new = screen.copy()
pygame.display.update()
cdef x = 0
startTime = time.time()
while x<100: # main game loop
pxarray_orig = pygame.PixelArray(screen.copy())
pxarray_new = pygame.PixelArray(screen.copy())
init_color_value = 16777215
for i in range(1,WIDTH-1):
for j in range(1,HEIGHT-1):
my_state = pxarray_orig[i,j]
neib_position = [[i,j+1],[i,j-1],[i+1,j],[i-1,j],[i+1,j+1],[i-1,j-1],[i+1,j-1],[i-1,j+1]]
neighbor_state = 0
for position in neib_position:
neighbor_state += pxarray_orig[position[0],position[1]]
if my_state == 16777215 and (neighbor_state == init_color_value*2 or neighbor_state==init_color_value*3):
pxarray_new[i,j] = 16777215
elif my_state == 16777215 and neighbor_state > init_color_value *3:
pxarray_new[i,j] = 0
elif my_state == 16777215 and neighbor_state < init_color_value*2:
pxarray_new[i,j] = 0
elif my_state == 0 and neighbor_state==init_color_value*3:
pxarray_new[i,j] = 16777215
screen_new = pygame.PixelArray.make_surface(pxarray_new).copy()
screen.blit(screen_new,screen.get_rect().topleft)
pygame.display.flip()
for event in pygame.event.get():
if event.type == QUIT:
pygame.quit()
sys.exit()
fpsClock.tick(FPS)
x += 1
executionTime = (time.time() - startTime)
print('Execution time in seconds: ' + str(executionTime))
pygame.quit()
输出为：
经过测试发现，上述代码的速度提升了约100%!当然，这对画面的流畅性提升明显，但仍然有显著的卡顿感。然而，上面的代码与纯Python代码相比只是增加了变量类型的定义。按照Cython的宣传，其提供的不仅仅是C语言一半的性能，而是几乎一模一样的性能。由于Pygame的内部的数据类型我们无法直接定义，因此目前来无法充分发挥Cython的全部优势。
为了获得更加流畅的画面，可以尝试用C++来实现生命游戏。C++是许多大型游戏引擎的开发语言，也是许多物理仿真的引擎，在进行流体仿真、软体仿真方面具有巨大优势。许多大型多主体仿真（agent-based modeling）需要用到超级计算机，其模型的实现语言就要求是C++，而不是Python或者JAVA等。
本文给出两种C++的实现，一种是利用向量加指针的方式来控制像素，一种是通过结构体的引用。经过测试，两种方式的画面流畅感差不多。
向量加指针的实现代码如下：
#include
const
int WIDTH = 600, HEIGHT=600;
//template
const sf::Color WHITE = sf::Color(255, 255, 255, 255);
const sf::Color BLACK = sf::Color(0, 0, 0, 255);
class ParticleSystem : public sf::Drawable, public sf::Transformable
{
public:
ParticleSystem(unsigned int count) :
m_vertices(sf::Points, count)
{
}
void generate_particle_matrix(int count)
{
std::vector<:color> m_color_row;
std::vector<:vertex> m_particle_row;
for (int j = 0; j < HEIGHT; ++j) {
m_particle_row.clear();
m_color_row.clear();
for
(int k = 0; k < WIDTH; ++k) {
sf::Vertex *p = &m_vertices[j * WIDTH + k];
sf::Vector2f position = { float(j),float(k) };
(*p).position = position;
if (std::rand() % 100 < 50) {
(*p).color = WHITE;
}
else {
(*p).color = BLACK;
}
m_particle_row.push_back(p);
m_color_row.push_back((*p).color);
}
m_particle_matrix.push_back(m_particle_row);
m_particle_colors.push_back(m_color_row);
}
}
void update()
{
for (int j = 1; j < HEIGHT - 1; ++j) {
for (int k = 1; k < WIDTH - 1; ++k) {
m_particle_colors[j][k] = (*m_particle_matrix[j][k]).color;
}
}
for (int j = 1; j < HEIGHT-1; ++j) {
for (int k = 1
; k < WIDTH-1; ++k) {
int neighbor_position[8][2] = { {j,k + 1}, {j,k - 1}, {j + 1,k}, {j + 1,k + 1}, {j + 1,k - 1}, {j - 1,k}, {j - 1,k + 1}, {j - 1,k - 1} };
sf::Color& my_color = (*m_particle_matrix[j][k]).color;
int color_count = 0;
for (int l = 0; l < 8;++l) {
if (m_particle_colors[neighbor_position[l][0]][neighbor_position[l][1]] == WHITE){
color_count++;
}
}
if (my_color == WHITE && (color_count == 2 || color_count==3)) {
my_color = WHITE;
}
else if (my_color == WHITE && color_count>3) {
my_color = BLACK;
}
else if (my_color == WHITE && color_count < 2) {
my_color = BLACK;
}
else if (my_color == BLACK && color_count == 3) {
my_color = WHITE;
}
}
}
}
private:
virtual void draw(sf::RenderTarget& target, sf::RenderStates states) const
{
states.transform *= getTransform();
states.texture = NULL;
target.draw(m_vertices, states);
}
private:
struct Particle
{
sf::Vector2f velocity;
sf::Time lifetime;
};
sf::VertexArray m_vertices;
std::vector<std::vector<:vertex>> m_particle_matrix;
std::vector<std::vector<:color>> m_particle_colors;
};
int main()
{
sf::
RenderWindow window(sf::VideoMode(WIDTH, HEIGHT), "Life Game");
ParticleSystem particles(WIDTH*HEIGHT);
particles.generate_particle_matrix(WIDTH*HEIGHT);
sf::Clock clock;
while (window.isOpen())
{
sf::Event event;
while (window.pollEvent(event))
{
if (event.type == sf::Event::Closed)
window.close();
}
particles.update();
window.clear();
window.draw(particles);
window.display();
}
return 0;
}
结构体加引用的实现方式：
#include
const int WIDTH = 500, HEIGHT = 500;
const sf::Color WHITE = sf::Color(255, 255, 255, 255);
const sf::Color BLACK = sf::Color(0, 0
, 0, 255);
class ParticleSystem : public sf::Drawable, public sf::Transformable
{
public:
ParticleSystem(unsigned int count) :
m_particles(count),
m_vertices(sf::Points, count)
{
}
void generate_particle_matrix(int count)
{
for (int j = 0; j < HEIGHT; ++j) {
for (int k = 0; k < WIDTH; ++k) {
sf::Vector2f position = { float(j),float(k) };
Particle& p = m_particles[j * WIDTH + k];
sf::Vertex& v = m_vertices[j * WIDTH + k];
p.position = position;
v.position = p.position;
if (std::rand() % 100 < 50) {
p.my_color = WHITE;
}
else {
p.my_color = BLACK;
}
v.color = p.my_color;
}
}
}
void update()
{
for (int j = 0; j < HEIGHT; ++j) {
for (int k = 0; k < WIDTH; ++k) {
if (j - 1 > 0 && j + 1 < HEIGHT && k - 1 > 0 && k + 1 < WIDTH) {
Particle& p = m_particles[j * (WIDTH)+k];
sf::Vertex& v = m_vertices[j * (WIDTH)+k];
int neighbor_position[8] = { j*WIDTH+k + 1, j*WIDTH+k - 1, (j + 1)*WIDTH+k,
(j + 1)*WIDTH+k + 1, (j + 1)*WIDTH+k - 1, (j - 1)*WIDTH+k, (j - 1)*WIDTH+k + 1, (j - 1)*WIDTH+k - 1};
int color_count = 0;
for (int l = 0; l < 8; ++l) {
if (m_vertices[neighbor_position[l]].color == WHITE) {
color_count++;
}
}
if (v.color == WHITE && (color_count == 2 || color_count == 3)) {
p.my_color = WHITE;
}
else if (v.color == WHITE && color_count > 3) {
p.my_color = BLACK;
}
else if (v.color == WHITE && color_count < 2) {
p.my_color = BLACK;
}
else if (v.color == BLACK && color_count == 3) {
p.my_color = WHITE;
}
}
}
}
for (int i = 0; i < WIDTH * HEIGHT; ++i) {
m_vertices[i].color = m_particles[i].my_color;
}
}
private:
virtual void draw(sf::RenderTarget& target, sf::RenderStates states) const
{
states.transform *= getTransform();
states.texture = NULL;
target.draw(m_vertices, states);
}
private:
struct
Particle
{
sf::Vector2f position;
sf::Color my_color;
};
std::vector m_particles;
sf::VertexArray m_vertices;
std::vector<std::vector<:vertex>> m_particle_matrix;
std::vector<std::vector<:color>> m_particle_colors;
};
int main()
{
sf::RenderWindow window(sf::VideoMode(WIDTH, HEIGHT), "Life Game");
ParticleSystem particles(WIDTH * HEIGHT);
particles.generate_particle_matrix(WIDTH * HEIGHT);
sf::Clock clock;
while (window.isOpen())
{
sf::Event event;
while (window.pollEvent(event))
{
if (event.type == sf::Event::Closed)
window.close();
}
particles.update();
window.clear();
window.draw(particles);
window.display();
}
return 0;
}
为了对比C++和Python的输出效果，我们将两者的输出并列放置，其中Python的代码已经经过Cython编译。下面是不同像素的下两者效果的比较。虽然动图截图也有其FPS限制，但仍然也可以看出即使是经过Cython编译过，利用Python的流畅性仍然远低于C++。
C++400*400 VS Python400*400
C++500*500 VS Python 400*400
C++600*600 VS Python 400*400
C++ 600*600 VS Python 600*600
可以看出，即使处理600*600的像素，C++的流畅度仍然超过Python处理的400*400。而Python在600*600的像素下可以明显看出力不从心。
结论：
C++不愧是支持超算的语言，在处理大规模仿真实体的条件下，将Python远远抛开。但是Python在数据分析和人工智能领域的主导地位仍然让其具有强大的生态优势。因此，如果模型要处理的问题规模很大，且对动态可视化要求较高，则可以使用C++进行实现；但可以用Python作为数据分析的主要工具。同时，Python还可以用来做模型原型的开发和探索，其开发效率和难度要远低于C++。因而，两种语言对于建模来说都有其独特优势，在建模的过程中应取长补短。
注：全文转载自科学建模DrZ的头条专栏，详情点击阅读原文。
