渗流模拟 — 把全书零件组装成一个程序Copy Anchor Link

这是全书的收官章:一个从头写到尾的完整程序。 题目叫渗流(percolation)—— 大意是「森林火灾模拟」:在棋盘上随机种树,沿着一边点火,看火能不能烧穿到对面。 它的迷人之处在于存在一个临界概率:树的密度超过它,火几乎必然蔓延;低于它,火几乎必然熄灭。 读完你会看到:前 18 章的每一件工具,都在这里找到了岗位。

1. 这一章讲什么Copy Anchor Link

问题一句话就能说清:¹

拿一块棋盘,每个格子随机种上树或留成沙漠;沿着一边把所有树点燃,然后看着火蔓延。 如果种树的概率超过某个临界值,火会烧到另一边;低于它,火会先熄灭。

物理背景(书里没展开,补充一句):渗流是统计物理的经典模型,用来研究「随机介质中连通性何时突然涌现」—— 森林火灾、多孔材料渗水、流行病传播,背后是同一套数学。它最出名的结论是相变: 宏观行为在某个人为无法解析求出的临界密度处突变。

  程序的组装顺序(也是任何仿真的标准顺序):
  ① 定义世界的词汇(enum:四种地块)
  ② 生成世界(随机撒树)
  ③ 定义局部规则(邻居、点火、烧尽)
  ④ 主循环(只要火还烧着,就推进一步)
  ⑤ 输出(把地块数值化成颜色,画出来)
  图说:每一步都只依赖前面已经建好的东西——这就是「自底向上组装」。

2. 核心原理:一步步把世界搭出来Copy Anchor Link

2.1 词汇:enum 定义四种地块Copy Anchor Link

解决什么问题: 棋盘的每个格子是四种状态之一:沙漠(不可燃)、树(可燃)、燃烧中、灰烬(已烧尽)。 用数字 0–3 记也行,但读者看到 map[i][j] == 2 要回去查 2 是什么——魔法数字,第 06 章刚判过反模式。

怎么做: 用 enum(第 06 章的防守工具):²

from enum import Enum

class Field(Enum):
    DESERT = 0
    TREE = 1
    BURNING = 2
    ASH = 3

之后代码里出现 Field.BURNING,语义自带;想拼错成 Field.BURN 会当场报错。 书里明说:「这让代码比用魔法数字好读得多」——第 06 章的预言兑现。

2.2 生成世界:随机撒树Copy Anchor Link

怎么做: 密度参数 probability_tree(书的主程序用 0.55,50×50 的棋盘), 每个格子独立地 random.random() < probability_tree 决定种树还是留沙漠(第 09 章的均匀随机数):³

def initialize_map(probability_tree, size):
    map = []
    for _ in range(size):
        row = []
        for _ in range(size):
            if random.random() < probability_tree:
                row.append(Field.TREE)
            else:
                row.append(Field.DESERT)
        map.append(row)
    return map

嵌套列表 = 二维矩阵(第 04 章),外层循环走行、内层走列——第 08 章的「嵌套循环预算 = 维度数」的教科书场景:两层循环,没有额外业务逻辑。

2.3 局部规则:邻居、点火、蔓延Copy Anchor Link

坐标打包成类。 x 是行号(上→下),y 是列号(左→右);坐标总是成对出现, 按第 11 章的原则打包成一个小类:⁴

class Point:
    def __init__(self, x, y):
        self.x = x
        self.y = y

邻居与周期边界。 periodic_neighbors 找上下左右四个邻居。特别的是左右方向用了周期边界条件 (periodic boundary conditions):棋盘在水平方向「卷成圆筒」——最左列的左邻居是最右列。 实现用的是取余运算 %:y=0 的左邻居是 y=49,y=49 的右邻居是 y=0。⁵

  普通边界:          周期边界(左右):
  □ □ □             □ □ □
  □ ■ □   ← 左邻居不存在  ■ 的最右邻居是最左 □ ←
  □ □ □             (上下保持普通边界:顶行无邻居)
  图说:% 运算把 y-1 = -1 折回 size-1,y+1 = size 折回 0——一行代码换来的拓扑。

三条蔓延规则,各是一个小函数,职责单一:⁶

第三步值得停一拍:它造一个新棋盘,而不是原地改——因为「正在烧的树既要点燃邻居、又要变灰」, 如果边遍历边改,规则会互相污染。这是仿真代码的经典纪律(我们的评注;书里只展示了「返回新矩阵」的做法)。 书里也坦白:「真要做计算密集的应用,这步可以优化;对我们这里,够用。」

2.4 主循环:while 是关键词Copy Anchor Link

def still_burning(map):
    # 棋盘上还有任何 BURNING 格子吗?
    ...

def perculate(map):
    while still_burning(map):        # while:条件成立就一直做(第 08 章没教,这里补上)
        map = perculation_step(map)
    return map

while 循环:先检查条件、再执行体,条件假了才停——「让火烧尽」这种事先不知道要几步的场景, 正是 for(已知步数)与 while(已知条件)的分工线。点火只需 ignite_front_row:把第一行的树全置 BURNING。⁷

程序入口收在标准姿势里:

def main():
    map = initialize_map(0.55, 50)
    ignite_front_row(map)
    map = perculate(map)

if __name__ == "__main__":
    main()

if __name__ == "__main__": 的意思是「只有直接运行这个文件时才执行 main;被别的文件 import 时不执行」—— 第 10 章 sys.argv 一节出现过同一用法,这是它的正式登场。⁷

2.5 输出:把世界画出来Copy Anchor Link

怎么做: 两种表示的桥接:enum → 数值(plotable_map 把每格换成 .value)→ 颜色。 颜色用 ListedColormap 自定义四色:⁸

colors = ["yellow", "green", "red", "black"]    # 沙漠黄、树绿、火红、灰烬黑
cmap = matplotlib.colors.ListedColormap(colors)

plt.imshow(plotable_map(map), cmap=cmap)        # imshow:把二维数值表当图显示
plt.show()

主程序点火前先画一张、烧完再画一张:黄绿棋盘 → 烧出一道黑色火廊。 (书里的小补丁 fix_color_schema 故意把 [0][0] 改成灰烬,保证四种颜色都在色标里出现——演示性技巧,不必学。)

3. 走一遍:一次完整运行Copy Anchor Link

把主线例子从头走到尾,每步旁边是具体状态(密度 0.55、棋盘 50×50 取自书的主程序):⁷

  ① initialize_map(0.55, 50)   → 约 55% 的格子是 TREE,其余 DESERT(0.55 是概率,实际 1375 棵左右,随种子浮动)
  ② ignite_front_row           → 第一行的树全部 BURNING(约 27 棵,按种子的随机结果)
  ③ while still_burning:
        perculation_step       → 每轮:邻火树变 BURNING、旧火树变 ASH
                                 火沿树带向对面推进;沙漠是不可逾越的隔离带
  ④ 火无处蔓延 → still_burning 返回 False → 循环结束
  ⑤ imshow 画出:黑色灰烬区从顶边延伸——0.55 高于渗流阈值,烧穿棋盘
  图说:把 0.55 改成 0.4 再跑,火几乎总在一小片就灭——临界现象亲眼可见。

临界概率是多少? 对这种方形网格近邻渗流,理论值约 0.5927(补充,不在书里)—— 0.55 其实低于它,严格说「烧穿」不是必然;但在有限 50×50 棋盘上,靠近阈值就能频繁观察到穿与不穿两种结局。 这正是蒙特卡洛实验的素材:重复跑 N 次、统计烧穿比例——把本章和第 09 章的 π 估算接起来, 就是一道完整的科研入门作业(书没布置,留给你)。

4. 作者的判断与证据Copy Anchor Link

书里给了证据的: 全部函数可直接运行;两张图(点火前/烧尽后)是程序的直接产物; 「临界概率以上蔓延、以下熄灭」是渗流理论的标准结论,书里以定性方式陈述。¹

作者的判断:

• 「这个算法可以优化,但够用」——第 07 章「复杂度 = 问题复杂度」的得体应用;
• enum 替代魔法数字——他的贯穿立场,最后一次示范;
• 「动手算一遍 (y−1+l) % l 就会懂,没有比这更好的解释」——他一贯的「自己算一遍」教学法。⁵

判断(我们的,不是书里的): 这个例子真正的教学价值不在渗流,而在程序的形状: 词汇(enum)→ 世界(数据结构)→ 规则(纯函数)→ 循环(推进)→ 可视化(输出), 五步骨架适用于绝大多数仿真(元胞自动机——格点上按局部规则同步更新的模型,渗流就是它的亲戚、分子动力学、排队模型……)。 学会看「程序的形状」之后,读陌生仿真代码的速度会快一个量级。 如果错,会错在: 交互式程序、事件驱动系统(游戏、GUI)的主循环不是「推进世界」而是「响应事件」, 骨架不同——把本章的五步当成所有程序的形状,会在那些领域读歪。

5. 边界与局限Copy Anchor Link

• 没有自动化的实验:临界概率是「说」出来的,没有「跑 1000 次统计烧穿率」的代码(我们在第 3 节给了思路);
• 邻居只算四近邻(上下左右),不含对角——这是模型选择,书里未讨论其他选法的差别;
• 周期边界只在一个方向,不对称的原因书里没解释(实验设计的选择);
• 代码里有几处转印级小错(如 neighbors.append(... % l) 的变量 l 未定义、periodic_neighbors 的拼写), 照抄前要自己修——也算一个免费的「读代码找茬」练习。⁵

6. 可带走的Copy Anchor Link

本章主线例子,一行写完: 0.55 密度撒树 → 顶行点火 → while still_burning: step() → 黑色灰烬烧穿棋盘;改 0.4 再跑,火一片即灭。⁷

1. 仿真的五步骨架:词汇(enum)→ 世界(嵌套列表)→ 规则(纯函数)→ 主循环(while)→ 可视化(imshow)。
2. enum 又一次替代魔法数字:四种地块语义自带、拼写错误当场报错。
3. 坐标成对出现就该打包成类(Point)。
4. 周期边界 = 取余运算:一行代码把棋盘卷成圆筒。
5. 逐步推进要造新状态,别边遍历边改(防规则互相污染)。
6. while 管「不知要几步、只知停止条件」的活;for 管已知步数的活。
7. if __name__ == "__main__": 是脚本入口的标准姿势,import 时不执行。
8. imshow + 自定义色表 = 把二维数值表变成世界地图的最短路径。
9. 临界现象是科研素材:重复跑、统计比例,就是蒙特卡洛实验(衔接第 09 章)。

7. 原文地图Copy Anchor Link