语言程序计算n重贝努利试验模拟_在Processing和Taichi中模拟《我的世界》里的奇妙水格格...

本项目使用了python模式的Processing，关于交互程序编程软件Processing的介绍请自行搜索或者查看本专栏的早期文章。

在Processing教学名书《the nature of code》中，有提到过细胞自动机和“生命游戏”，非常有趣。还有人改写相关代码，并做成了线上游戏：https://sandspiel.club/ （强烈推荐游玩一下）：

sandspiel

其实细胞自动机的这种离散格子的形式，是现在主要的一种物理仿真方法（有限元分析法）的简略版本。比如做机械仿真的时候，就会把机械分为好多个小格子，然后格子之间按照设定好的规则进行相互作用，就能非常接近实际地计算出机械在各种情景下的表现。

而且除了在科研领域，在游戏领域也大量使用到的这种技术。尤其是对于游戏场景中“水”的运动的仿真，也是和这种方法密切相关的。比如说我的世界里面的水，就是一个很好的典型。

而且大部分我的世界玩家，相信都观察到了，“我的世界”里面的水等流体的体积都是“不守恒”的，通过这种特性就有可能源源不断地获得水这种材料。这种不守恒，是计算机程序的速度与精度之间的妥协。

编程思路

将水抽象为一个类Grid。Grid类中最重要的数据就是一个二维数组，实时记录了各个格子中的“水量”。Grid还有一个方法，能够把自己在Processing中画出来。

而Grid最关键的方法是更新方法。因为二维自动机中的格子一般都只和相邻位置的四个或者8个格子相互作用，所以我们也遵循这样的原则（这样能够大大提升速度）。那么剩下的就是，具体该怎么设定一个格子和周围格子的交互规则了。

如果按照一般的思路，一个格子如果是只出水的话，那么周围格子所能得到的总水量一定和他给出的总水量是相等的。如果想要采取这种方式，就需要有一个统计周围需求，然后计算出水总量，然后再将出的水分别划分给周围的格子。

沿着这种思路，很容易得到一种编程思路，就是“每一个格子的下一个值都是周围格子的平均值”的方法。这是一种经典的方法和思路，我第一次尝试的就是这种方法。然而发现这种模拟的效果并不理想（行为非常单调，像是在模拟一团光，而没有水的特征）。

这种方法，其实无法保证水总量是严格恒定的，而是只能维持总量在一个大致稳定的状态。而且这种方法并不是很实际，水作为一种物质，显然无法掌握“求均值”这种高级的运算。

而如果我们舍弃了这种“总分总”的麻烦逻辑，而是采用一种更加“分布式”的方法： 只要两个相邻的格子a和b，a比b的水少，那么a就直接从b得到一部分水，这部分水是ab水量差的一部分。然后将所有的相邻的格子都两两比较，就能得到结果。这种方式更接近物理实际（有点微分方程的感觉了），而且计算也更加高效。

需要注意的是，在使用上一个时间的数据计算下一个时间的数据，需要对所有格子遍历，在遍历的过程中，要把所有的结果都暂存在一个新的数组中。等到所有的格子数据都计算完毕，再用新数组替换掉旧数组，这样能保证各个格子是“同步”地更新的。

程序效果

知乎视频​www.zhihu.com

交互：鼠标左键拖动加水，右键拖动除水。

效果：可以看出液体格子的行为非常逼真，不仅能从水多的地方流到水少的地方，甚至还有一点“荡漾”。

展望：这种模拟往往计算量非常大，所以一般需要使用专业的软件，所以后续我应该不会对这个项目继续进行迭代了。但是相信文本还是有教育和科普意义的，顺便还能带大家熟悉一下Procesing的编程环境。

注意：代码中有一个关键参数（在代码注释中标出），决定了水的“荡漾程度”，它小的时候虽然水不太荡漾了，但是水流的速度也慢了。而如果它大得超过一个临界的话，水就会变成“无药可救狗眼闪瞎土嗨水”，所及程序中的这个参数其实已经是调了很多次之后的很理想的状态了。

代码

格子类：

zhous = [(x, y) for x in (-1, 0, 1) for y in (-1, 0, 1)]

from copy import deepcopy

class Grid():
    def __init__(self, xn, yn, wid, heig):
        self.data = [[0.0] * wid for i in range(heig)]
        # self.data[3][3] = 125*xn*yn
        self.xn = xn
        self.yn = yn
        self.wid = wid
        self.heig = heig
    
    def update(self):
        "更新状态"
        tmp = 0
        tmp_data = deepcopy(self.data)
        for i in range(self.xn):
            for k in range(self.yn):
                for dx, dy in zhous:
                    tx, ty = i+dx, k+dy
                    if 0<= tx < self.xn and 0<= ty < self.yn:
                        tmp = self.data[i][k]
                        tmp2 = self.data[tx][ty]
                        d = abs(tmp - tmp2) * 0.165  # 这里的这个值0.165，设的越大水越“荡漾”
                        if tmp > tmp2:
                            tmp_data[i][k] -= d
                        else:
                            tmp_data[i][k] += d
                        
        self.data = tmp_data   
                
        if frameCount % 10 == 0:  # 定时输出帧率
            print(frameRate)
    
    def get_touch(self, x, y, p=True):
        "对鼠标做出响应"
        m = int(float(x)/self.wid*self.xn)
        n = int(float(y)/self.heig*self.yn)
        print(m, n)
        if p:
            self.data[m][n] += 500
        else:
            self.data[m][n] -= 500
    
    def display(self):
       "把自己画出来"
        w = self.wid / self.xn
        h = self.heig / self.yn
        for i in range(self.xn):
            for k in range(self.yn):
                pushMatrix()
                translate(w*i, h*k)
                fill(self.data[i][k])
                stroke(100)
                rect(0, 0, w, h)
                popMatrix()

主程序：

from Grid import Grid

def setup():
    global gr
    size(800, 600)
    gr = Grid(30, 30, width, height)
    
def draw():
    background(200)
    gr.update()
    gr.display()


def mouseDragged():
    gr.get_touch(mouseX, mouseY, mouseButton == LEFT)

突然硬核了起来：使用Taichi加速

Taichi是一种用于计算机图形应用的高性能编程语言。其设计目标有：

• 生产力
• 性能
• 可移植性
• 空间稀疏计算
• 可微编程
• 元编程

我记得大约是半年还是一年前，taichi刚出来的时候，真是火遍了知乎，堪称是计算机领域技术“出圈”的典范，甚至连“高中生看了taichi觉得自己一辈子也做不出这样的项目，不想努力了怎么办”这种知乎提问都大火了一把。（这一段的意思就是taichi很火啦）

我对Taichi也是门外汉，粗略地去理解的话，感觉上有点像是更容易使用的着色器语言（GLSL），底层向量和矩阵的使用感觉又有点像numpy，还支持GPU加速或者物理仿真啥的？总之就是各种又快又好用吧~

下面程序使用了Python3环境中的taichi，果然加速了很多：（仿真逻辑和Processing中的是一致的，我从taichi官方的一个生命游戏的例子中尝试改改，没想到还真改出来了）

taichi官方的生命游戏程序截图

知乎视频​www.zhihu.com

（鼠标拖动交互的时候仿真会暂停，需要按下空格键之后继续仿真）

可以看出，在同样的机器上，将仿真数量从30*30扩展到50*50之后，帧率已经从原始的10以内帧数每秒的龟速提升到“任意流畅”的地步，也就是50帧左右，并且经过后续测试，继续提升仿真数量到100×100之后依旧能非常流畅。这就是taichi的魅力吧。

至于写taichi的时候，有一点心得就是：尽量不要在 ti.kernel装饰的函数里面写复杂的循环，不然会编译出错。

代码：

import taichi as ti
# import numpy as np

ti.init()

n = 50
cell_size = 12
img_size = n * cell_size
alive = ti.var(ti.f32, shape=(n, n))  # alive = 1, dead = 0
img = ti.var(ti.f32, shape=(img_size, img_size))  # image to be displayed


zhous = [(x, y) for x in (-1, 0, 1) for y in (-1, 0, 1)]


@ti.func
def calc_rule(i, k):
    return alive[i, k] + 0.18 * (sum(alive[i+a, k+b] for a, b in zhous) - 9 * alive[i, k])


@ti.kernel
def run():
    for i, j in alive:
        alive[i, j] = calc_rule(i, j)


@ti.kernel
def init():
    for i, j in alive:
        alive[i, j] = ((float(i)/n + float(j)/n)-1)*200
        # alive[i, j] = 0


@ti.kernel
def render():
    for i, j in alive:
        for u, v in ti.static(ti.ndrange(cell_size, cell_size)):
            img[i * cell_size + u, j * cell_size + v] = alive[i, j] / 100 + 0.5


gui = ti.GUI('Amazing Water Sim', (img_size, img_size))

print('[Hint] Press `r` to reset')
print('[Hint] Press SPACE to pause')
print('[Hint] Click LMB, RMB and drag to add alive / dead cells')

init()
paused = False
while gui.running:
    for e in gui.get_events(gui.PRESS, gui.MOTION):
        if e.key == gui.ESCAPE:
            gui.running = False
        elif e.key == gui.SPACE:
            paused = not paused
        elif e.key == 'r':
            alive.fill(0)

    if gui.is_pressed(gui.LMB, gui.RMB):
        mx, my = gui.get_cursor_pos()
        alive[int(mx * n), int(my * n)] = (gui.is_pressed(gui.LMB) - 0.5) * 2000
        paused = True

    if not paused and gui.frame % 4 == 0:
        run()

    render()

    gui.set_image(img)
    gui.show()

最后再加两个Taichi官方例子中，两个特别像是用GLSL写出来的项目的效果动图吧：