元胞自动机-附代码注释

元胞自动机（Cellular Automata CA）是一种时空离散的局部动力学模型，是研究复杂系统的一种典型方法，适用于空间复杂系统的时空动态模拟研究。
元胞自动机不是由严格定义的物理方程或函数确定，而是用一系列模型构造的规则构成。满足这些规则的模型都可算作是元胞自动机模型。故元胞自动机是一类模型的总称，或者说是一个方法框架。
在CA模型中，散布在规则格网(Lattice Grid)中的每一元胞(Cell)取有限的离散状态，遵循同样的作用规则，依据确定的局部规则作同步更新。大量元胞通过简单的相互作用而构成动态系统的演化。

元胞自动机的组成部分有元胞、元胞空间、邻居及规则四部分。
简单讲，元胞自动机可以视为由一个元胞空间和定义于该空间的变换函数所组成。

元胞自动机的元胞

元胞是元胞自动机的最基本组成部分。它分布在离散的一维、二维或多维欧几里德空间的晶格点上。
状态，状态可以是{0,1}的二进制形式，或是整数形式的离散集。严格意义上，元胞只能有一个状态变量，但在实际应用中，往往将其进行扩展。

元胞自动机的元胞空间

元胞空间元胞所分布在的空间网点集合就是元胞空间。
对于一维元胞自动机，元胞空间划分只有一种。而高维的元胞自动机，元胞空间的划分则可能有多种形式。对于常见的二维自动机，元胞空间通常可按三角形、四边形或六边形三种网格排列。

• 三角网格：拥有较少的邻居数目，这在某些时候很有用。缺点是计算机的表达与显示不方便。
• 四边形网格 （最常用）：直观简单,特别适合于计算机环境下进行表达显示。
• 六边形网格：能较好的模拟各向同性的现象，因此，模型能更加自然而真实。其缺点同正三角网格一样，在表达显示上较为困难和复杂。

在理论上，元胞空间的边界条件在各个维向上是无限延展的，但实际应用过程中，无法在计算机上实现这一理想条件。

元胞自动机的邻居

在数学上，它通常可以表示为一个多维的整数矩阵。
邻居、元胞和元胞空间只表示了系统的静态成分，为了将动态引入系统，必须加入演化(规则。这些规则是定义在局部空间范围内的，即一个元胞下一时刻的状态决定于本身的状态和它的邻居无胞的状态。因此，在指定规则之前，必须定义一定的邻居规则，明确哪些元胞属于该元胞的邻居。
一维元胞自动机中，通常以半径r来确定邻居，距离一个元胞r内的所有元胞都属于该元胞的邻居。
二维元胞自动机的邻居定义较为复杂，但通常有以下几种（以正方形网格为例)

元胞自动机的规则

规则根据元胞当前状态及其邻居状况确定下一时刻该元胞状态的动力学函数。简单讲，就是状态转移函数。这个函数构造了一种简单的离散的时间和空间范围的局部物理成分。状态的变化可以由状态转移函数表示。t时刻的邻居状态组合时间元胞自动机是一个动态系统，它在时间维上的变化是离散的，即时间t是一个整数值，而且连续且等间距。在转换函数中，一个元胞在t+1时刻的状态只直接决定于t时刻的该元胞及其邻居的状态。

元胞自动机的特性

把一个空间划分成网络，每一个点表示一个元胞，它们的状态赋值，在网格中用颜色的变化来表示，在事先设定的规则下，元胞的演化就用网格颜色的变化来描述，这样的模型就是元胞自动机。
通过对元胞自动机这些网络中的格点的不同定义，以及初始条件的不同，可以模拟出不同的现象和过程。
离散性：元胞自动机是高度离散的。它不仅仅空间离散时间离散，而且在函数值，即元胞的状态值也是离散的。
动力学演化的同步性：元胞自动机具有利用简单的，局部规则的和离散的方法，描述复杂的，全部的和连续系统的能力。
相互作用的局部性：元胞自动机的规则是局部的，而动力学行为规则是全局的，在模拟的过程中，具体的演化过程也是局部的，即仅同周围的元胞有关系。
时间、空间、状态都离散，每个变量只取有限多个状态，且其状态改变的规则在时间和空间上都是局部的。

元胞自动机的应用思想

复杂系统又称为非线性系统。如城市的发展与演化、城市人流与交通流以及交通堵塞的形成、自然环境下的动物的空间分布、河网的形成、疾病的传播等。
传统的自顶向下的分析方法是把系统分割成几个部分，对每一个部分逐个进行研究。而目前提出来的分析复杂动态系统的思想：自底向上的研究方法。

元胞自动机的演化行为的统计特征

Wolfram将元胞自动机的演化行为归纳为四大类：
平稳型(homogeneous) ：自任何初始状态开始经过一定时间演化后，经过若干步运算便停留在一个固定的状态。
周期型(periodic)：经过一定时间演化后，在几种状态之间周期循环。
混沌型(chaos)：自任何初始状态开始，经过一定时间演化后，处于一种完全无序随机的状态，几乎找不到任何规律。
复杂型(edge of chaos)：在演化的过程中可能产生复杂的结构，这种结构既不是完全的随机混乱，又没有固定的周期和状态。
固定点->周期->复杂->混沌

其实元胞自动机更多的是相当于一种思想、一种框架，通过设置一套规则，将时间离散化后，对元胞进行规则下的随时间的演化。

至今元胞自动机已经有了非常广泛的应用，典型应用如交通流模型、人员疏散模型、森林火灾模型、传染病模型、土地利用演化模型、城市发展扩张、交通状况模拟、博物馆逃生模拟等等。

森林火灾模拟

%% 森林火灾
m = 300; n = 300; % 表示森林的矩阵行列 m x n

Plight = 5e-6;    % 闪电概率
Pgrowth = 1e-2;   % 生长概率

% 邻居方位 d 和点燃概率 p
d = {[1,0], [0,1], [-1,0], [0,-1]};
p = [    1,     1,      1,      1];

% % 改进元胞自动机
% d = {[1,0], [0,1], [-1,0], [0,-1], [1,1], [-1,1], [-1,-1], [1,-1]};
% p = [ones(1,4), ones(1,4)*(sqrt(1/2)-1/2)];

% % 考虑风的情况
% d = {[1,0], [0,1], [-1,0], [0,-1], [1,1], [-1,1], [-1,-1], [1,-1], [0,-2]};
% p = [ 0.80,  0.30,   0.80,   1.00,  0.12,   0.12,    0.30,   0.30,    0.8];

% 空=0, 火=1, 树=2
E = 0; F = 1; T = 2;

S = T * (rand(m, n) < 0.5);

imh = image(cat(3, S==F, S==T, zeros(m,n)));
axis image;

for t = 1:3000

    % 计算邻居中能传播着火的个数
    sum = zeros(size(S));
    for j = 1:length(d) % 计算四个方向的着火邻居
        sum = sum + p(j) * (circshift(S,d{j}) == F);
    end

    isE = (S==E); isF = (S==F); isT = (S==T);     % 找出三种不同的状态

    ignite = rand(m,n) < sum | (rand(m,n) < Plight);  % 着火条件

    % 规则 1: 着火
    Rule1 = T * (isT & ~ignite) + F * (isT & ignite);
    % 规则 2: 烧尽
    Rule2 = F*isF - F*isF;
    % 规则 3: 新生
    Rule3 = T * (isE & rand(m,n) < Pgrowth);

    S = Rule1 + Rule2 + Rule3;

    set(imh, 'cdata', cat(3, isF, isT, zeros(m,n)))
    drawnow
end

生命游戏模拟

%% 生命游戏
%% 界面控制
startbutton = uicontrol('style','pushbutton','string','start','fontsize',12,'position',[100,400,50,20],'callback','start=1;');
freezebutton = uicontrol('style','pushbutton','string','freeze','fontsize',12,'position',[320,400,50,20],'callback','freeze=1;');
quitebutton = uicontrol('style','pushbutton','string','quit','fontsize',12,'position',[400,400,50,20],'callback','stop=1;close;');
number = uicontrol('style','text','string','1','fontsize',12,'position',[20,400,50,20]);
freeze_at_text = uicontrol('style','text','string','freeze at：','fontsize',12,'position',[180,400,80,20],'callback','start=1;');
freeze_at = uicontrol('style','edit','string','1000','fontsize',12,'position',[250,400,50,20]);
%% 初始化设置
n = 400; % 元胞数目
z = zeros(n, n);
cells = (rand(n, n)) < 0.6; % 初始值随机分布
big_matrix = zeros(n+2, n+2); % 增补矩阵
big_matrix(2:end-1, 2:end-1) = cells;
imh = image(cat(3,z, z, cells)); % 建立图像句柄
figure(gcf);
set(gcf, 'doublebuffer', 'on'); % 开启双缓冲
%% 主循环
start = 0; stop = 0; freeze = 0;
while stop == 0
    if findobj == 0 % 用来监视窗口是否已经被关闭 防止报错
        break;
    end
    
    if start == 1
        z = zeros(n, n);
        sum = sum_neighbour(big_matrix); % 计算邻居状态和
        z(cells == 1 & (sum == 2 | sum == 3)) = 1; % 按规则更新状态
        z(cells == 0 & sum == 3) = 1;
        big_matrix(2:end-1,2:end-1) = z; % 将新状态填充回增补矩阵
        cells = z; % 将状态更新
        set(imh, 'cdata', cat(3, zeros(n,n), zeros(n,n), cells)); % 画图
        stepnumber = 1 + str2double(get(number, 'string')); % 更新计数
        set(number, 'string', num2str(stepnumber));
    end
    if (freeze == 1 || stepnumber == str2double(get(freeze_at, 'string'))) % 判断是否暂停
        start = 0;
        freeze = 0;
    end
    drawnow % 更新绘图
end
%% 计算邻居和函数
function [A] = sum_neighbour(B)
% 元胞自动机 计算周围邻居的状态和
    A = B(1:end-2,1:end-2)+B(1:end-2,2:end-1)+B(1:end-2,3:end)+B(2:end-1,1:end-2)+B(2:end-1,3:end)+B(3:end,1:end-2)+B(3:end,2:end-1)+B(3:end,3:end);
end