朝花夕拾

轻飘飘的旧时光就这么溜走 转头回去看看时已匆匆数年

碉堡的小程序:用 Python 制作演示迷宫算法的 gif 动画

2017-10-02


本文要介绍的是我写的一个有趣的小程序,一个脱离了低级趣味的程序,一个有益于广大人民学习算法的程序。代码在 Github 上。

先上效果图:

这个图初看起来让人眼花缭乱,不知所云,但是看完以后应该不难明白,它演示的是一个迷宫生成算法和一个迷宫求解算法 (只不过嵌入到了一张背景图片中)。实际上它演示的是概率论中的 Wilson 一致生成树算法 (UST 意为 uniform spanning tree) 和宽度优先搜索算法运行在一个 2D 网格图上的效果。还记得吗?一个图的生成树可以看作是这个图上的一个迷宫,图上的任何两个顶点之间沿着生成树有唯一的一条路径相连。

至于什么是 Wilson 一致生成树算法 (这里的重点在 "一致" 上),本文后面会介绍。

这个程序是怎么来的?

许多年前我在网上闲逛的时候,偶然发现了一个网站:https://bl.ocks.org/mbostock,当时我对上面展示的各种丰富炫酷的动态效果惊羡不已,尤其是其中对 Wilson 算法的演示,让我对此算法有了更直观和深入的理解。我立刻萌发了用 Python 制作一个 gif 版的动画演示的想法,但是思考了许久也不知道从何入手。这里困难的地方在于 Wilson 算法是一个随机算法,其运行时间是不确定的,一个动画里面可能包含数千帧,如果采用把每一帧保存为图像再合并到一起的话,最终得到的文件会非常庞大。而且这种纯暴力的做法逼格太 low,自诩为 geek 的我辈岂屑于为之?限于能力不足,这个想法只好被暂时压在心底,但是一直念念不忘。过了几年后,在一个偶然的机会我接触到了 gif 图像的编码协议,豁然开朗:为什么不直接把动画过程编码为字节流呢?通过精确定位每一帧的位置,控制 LZW 压缩过程的编码长度,文件过大的问题是可以解决的!前后捣鼓了半个月,反复研究协议细节,debug 了无数次后,这才作出了上面的效果。

后来我把这个程序作了改进,可以用来演示更多的 (理论上是全部) 迷宫生成/求解算法,比如:

  1. Prim 最小生成树算法:

  2. 随机深度优先搜索生成树算法:

  3. Kruskal 最小生成树算法:

    (免责声明:上图是个正经图,没有别的意思,不要把我批判一番)

其中染成白色的房间表示生成树,染成黑色的房间表示墙,任何两个白色房间之间有且只有一条路径相连 (根据树的定义)。

这个程序有如下特点:(容我厚颜无耻自吹自擂一番)

  1. 所有代码全部由纯 Python 写成,没有用到任何第三方库或者外部软件,也没有任何绘图/染色命令,仅使用了内置的 struct, random 模块和一些内置函数。后来的版本中为了显示进度条引入了 tqdm;为了把整个动画嵌入一张背景图片引入了 pillow,这些都属于特效,本质不需要。

  2. 实现了一个小型但高效的 gif 编码器,可以在数秒之内生成高度优化的动态图。这是这个程序最让人意外的一点:Python 生成图像的慢是出了名的,它居然能在几秒内生成一张包含几百乃至几千帧的 gif 动图?这是个大新闻啊!

  3. 严格遵循 GIF89a 协议,生成的图片在 chrome, firefox, IE 和 Eog 下都可以正常显示。

程序运行的相当快,生成一副 600x400 像素,演示 Wilson 算法的动图只要数秒,得到的文件包含 1000~3000 帧,但大小不超过 1M 左右。没想到吧?

其实程序最初的版本比现在要慢 5x~10x 倍左右,而且绝大部分功能都放在一个庞大的类里面 (20+ 个方法),现在的代码是经过不计其数次的修改,反复绞尽脑汁优化以后才得到的。

程序的核心部分是实现 Wilson 一致生成树算法和一个高效的 gif 编码器,所以除了 Wilson 一致生成树算法,还需要你理解 GIF89a 编码协议以及其中用到的 LZW 压缩算法。这些完全靠自己钻研的话还是很有难度的,要踩的坑非常多。由于篇幅所限这里我只简要介绍一下 Wilson 一致生成树算法,至于 GIF89a 协议你需要且只需要一份参考:What’s in a gif,我当时写这个程序的时候只参考了它。

什么是 Wilson 一致生成树算法

考虑如下问题:

问题:给定一个有限无向的连通图 \(G\),怎样在 \(G\) 的所有生成树中等概率地任选一种?(这样选出的生成树由于服从所有生成树上的一致分布,故而得名一致生成树)

你也许会想,这有何难?把所有的生成树都列举出来,然后从中随便挑一个不行吗?让我们以 \(n\) 阶的完全图 \(K_n\) 为例,根据 Cayley 公式它总共有 \(n^{n-2}\) 个不同的生成树,当 \(n=100\) 时这个数字是 \(100^{98}\),远远大于人类已知的宇宙中所有粒子的总数!

目前已知的算法里面,最有效的是 Wilson 在 1997 年的论文

Generating random spanning trees more quickly than the cover time

中提出来的,它是一个随机算法,这里随机算法的意思是说算法的运行时间是不确定的,有的时候你可能很幸运很快就找到了一个生成树,也有可能会一直等下去。不过算法以概率 1 会在有限时间内结束,而且大多数情况下它的表现还不错。

理解 Wilson 算法的关键概念是所谓的擦圈的随机游动 (loop erased random walk):一个图 \(G\) 上的随机游动就是从某个初始顶点 \(v_0\) 出发,每次从当前顶点 \(v_i\) 的邻居中完全随机地选择一个 (设为 \(v_{i+1}\)),然后移动到 \(v_{i+1}\)。这个随机游动的路径记作 \(\{v_0,v_1,\ldots,v_n,\ldots\}\)。图 \(G\) 上的擦圈的随机游动就是在随机游动的过程中,每当遇到一个先前访问过的顶点 \(u\) 时,立刻擦除路径中这两次访问 \(u\) 所形成的回路,然后从 \(u\) 继续作随机游动。

由定义可见擦圈的随机游动得到的路径中总不包含任何回路。

其实只要你看看开头的 gif 动态图,就很容易直观上理解擦圈的随机游动是怎么回事。

下面我们来介绍 Wilson 的一致生成树算法:

Wilson 一致生成树算法

  1. 选择任一顶点 \(v\) 作为根节点,维护一个树 \(T\),初始时刻 \(T=\{v\}\)
  2. \(z\) 是一个不属于 \(T\) 的顶点,从 \(z\) 出发作擦圈的随机游动,直到这个随机游动撞到 \(T\) 为止,然后把随机游动经过的路径加入 \(T\) 中。
  3. 重复步骤 2 直到 \(G\) 的所有顶点都属于 \(T\),这时得到的 \(T\) 就是一个一致生成树。

Wilson 算法的证明有相当的技巧性,本文就没法展开讲了,读者可以参考下面的文献:

Probability on Trees and Networks, by Russell Lyons and Yuval Peres.

或者我以前的一篇笔记

关于代码

理解代码的关键有这么几个:

  1. 由于 gif 图像的每一帧占据的是整个图像窗口的一个矩形子区域,在一个包含很多帧的动图中,相邻的两帧之间的变动可能很小,没有必要每次都将整个图像所占区域全部编码。我们只需要记录帧和帧之间的变化情况,得出每一帧所占的矩形子区域,每次编码时只针对这个子区域编码即可。这样就大大减小了生成的文件体积。

  2. 采用变长的 LZW 压缩算法。GIF89a 协议允许每个打包的数据块指明其所使用的最小码字的长度 (仅适用于本块数据),如果你确定这一帧图像用到的颜色很少,比如只有 4 种颜色,那么 2 个比特就足以表示这 4 种颜色,从而最小编码长度可以设置为 2。这样根据具体情况采用不同的编码长度能有效减少文件体积。

  3. 因为要频繁的进行字节流的操作,所以每次将编码后的数据先写入一个 BytesIO 对象中,放在内存里,最后一次性输出到硬盘。由于优化的好,动图大小普遍不超过 3M,所以放在内存里不是问题。

代码的组织结构是简单的三层论:顶层是抽象的 Maze 类,其本质就是一个 2D 网格图,用来跑各种图算法,它不关心动图的任何细节。底层是 GIFSurface 类,负责维护 gif 图片的全局信息,比如大小,循环次数,背景颜色,全局调色板等。中间层是 Animation 类,用来控制帧的信息,在算法运行过程中它按照一定的频率将 Maze 染色并编码写入 GIFSurface

2018/11/18 更新

加入了一个新脚本,演示二维的 Hilbert 空间填充曲线,方法是利用 Gray 码生成曲线上每个点对应的像素位置,结果如下:

此外还加入了一个演示 Langton 蚂蚁的例子: