（CA学习资料）美国康奈尔大学BioNB441在Matlab中的元胞自动机

网站翻译内容 美国康奈尔大学BioNB441在Matlab中的元胞自动机 介绍 元胞自动机（CA）是用于计算利用当地的规则和本地通信。普遍CA定义一个网格，网格上的每个点代表一个有限数量的状态中的细胞。过渡规则同时应用到每一个细胞。 典型的转换规则依赖于细胞和它的（4个或8个）近邻的状态，虽然临近的细胞也同样使用。 CA的应用在并行计算研究，物理模拟和生物模拟。这个页面将考虑如何写出高效的MATLAB代码的CA的实施和看一些有趣的规则。 Matlab代码注意事项 以下注意事项将说明使用Matlab程序计算康威的生命。部分的代码的解释如下。 矩阵和图像可以被转换为一个另一个，所以显示是为straighforward的的。如果阵列细胞包含二进制的每一个电池单体的状态和数组Z含有零，那么cat命令建立一个RGB图像，显示图像命令。图像命令也返回一个句柄。 imh = image(cat(3,cells,z,z)); set(imh, 'erasemode', 'none') axis equal axis tight 矩阵和图像可以被转换为一个另一个，那么初始条件，可以计算矩阵或图形命令。下面的代码生成一个元素为零的数组，细胞状态初始化为零，然后使细胞的交叉状态=1的数组的中心。使用（渗流群集）下面的例子中的一个图形命令来初始化CA阵列。 z = zeros(n,n); cells = z; cells(n/2,.25*n:.75*n) = 1; cells(.25*n:.75*n,n/2) = 1; Matlab代码进行矢量化，以减少开销。下面的两个语句计算的总和近邻的计算CA规则。代码利用Matlab的非常灵活的索引指定的近邻。 x = 2:n-1; y = 2:n-1; sum(x,y) = cells(x,y-1) + cells(x,y+1) + ... cells(x-1, y) + cells(x+1,y) + ... cells(x-1,y-1) + cells(x-1,y+1) + ... cells(x+1,y-1) + cells(x+1,y+1); cells = (sum==3) | (sum==2 & cells); 添加一个简单的GUI是很容易的。在这个例子中，三个按钮和一个文本字段implmented。三个按钮允许用户运行，停止，或选择退出。文本字段中显示模拟执行的步骤数。 %build the GUI %define the plot button plotbutton=uicontrol('style','pushbutton',... 'string','Run', ... 'fontsize',12, ... 'position',[100,400,50,20], ... 'callback', 'run=1;'); %define the stop button erasebutton=uicontrol('style','pushbutton',... 'string','Stop', ... 'fontsize',12, ... 'position',[200,400,50,20], ... 'callback','freeze=1;'); %define the Quit button quitbutton=uicontrol('style','pushbutton',... 'string','Quit', ... 'fontsize',12, ... 'position',[300,400,50,20], ... 'callback','stop=1;close;'); number = uicontrol('style','text', ... 'string','1', ... 'fontsize',12, ... 'position',[20,400,50,20]); 后控制（和CA）初始化，程序进入一个循环的测试中设置的每个按钮的回调函数的变量的状态下降。就目前而言，只是看在while循环，if语句的嵌套结构。程序循环，直到按下退出按钮。另外两个按钮导致当推一个if语句来执行。 stop= 0; %wait for a quit button push run = 0; %wait for a draw freeze = 0; %wait for a freeze while (stop==0) if (run==1) %nearest neighbor sum sum(x,y) = cells(x,y-1) + cells(x,y+1) + ... cells(x-1, y) + cells(x+1,y) + ... cells(x-1,y-1) + cells(x-1,y+1) + ... cells(3:n,y-1) + cells(x+1,y+1); % The CA rule cells = (sum==3) | (sum==2 & cells); %draw the new image set(imh, 'cdata', cat(3,cells,z,z) ) %update the step number diaplay stepnumber = 1 + str2num(get(number,'string')); set(number,'string',num2str(stepnumber)) end if (freeze==1) run = 0; freeze = 0; end drawnow %need this in the loop for controls to work End 例子 1.康威的生命。 其规则是： 总结8近邻 如果总和=2 那么状态不改变 如果总和=3，然后在状态= 1 否则状态= 0 代码： x = 2:n-1; y = 2:n-1; %nearest neighbor sum sum(x,y) = cells(x,y-1) + cells(x,y+1) + ... cells(x-1, y) + cells(x+1,y) + ... cells(x-1,y-1) + cells(x-1,y+1) + ... cells(3:n,y-1) + cells(x+1,y+1); % The CA rule cells = (sum==3) | (sum==2 & cells); 2.表面张力 其规则是： 总结8近邻和细胞本身 如果总和<4或总和= 5，则状态= 0 否则状态= 1 代码： x = 2:n-1; y = 2:n-1; sum(x,y) = cells(x,y-1) + cells(x,y+1) + ... cells(x-1, y) + cells(x+1,y) + ... cells(x-1,y-1) + cells(x-1,y+1) + ... cells(3:n,y-1) + cells(x+1,y+1)+... cells(x,y); % The CA rule cells = ~((sum< 4) | (sum==5)); 3.渗流集群 规则： 每个单元格的8个最近邻的求和（细胞是二进制值，0/1）。细胞也有一个独立的状态变量（称为'参观'）记录，如果他们曾经有过一个非零的前邻居。 计算0和1之间的随机数r。 如果总和> 0（至少一个相邻）的r>阈值，和小区从未有过的邻居小区= 1。 如果总和>0记录细胞，该细胞具有非零邻居“访问”的标志。 更新代码： The update code: sum(2:a-1,2:b-1) = cells(2:a-1,1:b-2) + cells(2:a-1,3:b) + ... cells(1:a-2, 2:b-1) + cells(3:a,2:b-1) + ... cells(1:a-2,1:b-2) + cells(1:a-2,3:b) + ... cells(3:a,1:b-2) + cells(3:a,3:b); pick = rand(a,b); cells = cells | ((sum>=1) & (pick>=threshold) & (visit==0)) ; visit = (sum>=1) ; 变量a和b的图像的大小。 >初始的图像是由图形操作。设立轴大小是固定的，下面的语句写文字转换成轴，然后抢轴的内容，并把它们放回一个数组，用的getFrame功能. ax = axes('units','pixels','position',[1 1 500 400],'color','k'); text('units', 'pixels', 'position', [50,255,0],... 'string','BioNB','color','w','fontname','helvetica','fontsize',100) text('units', 'pixels', 'position', [120,120,0],... 'string','441','color','w','fontname','helvetica','fontsize',100) initial = getframe(gca); [a,b,c]=size(initial.cdata); z=zeros(a,b); cells = double(initial.cdata(:,:,1)==255); visit = z ; sum = z; 4.激发介质（BZ反应或心脏） 规则： 细胞可以在10个不同的状态。休息状态0。国家1-5活跃，状态6-9难治。 计数8个最近邻的每个单元格是在一个激活状态。 如果总和大于或等于3（至少三个活性邻居的），那么细胞=1。 状态1〜9中，没有更多的输入逐步发生。如果状态=1，那么下一个状态=2。如果状态= 2那么下一状态= 3，同样的所有的州多达9个。如果状态= 9，那么下一个状态=0和细胞是在休息。 更新代码： x = [2:n-1]; y = [2:n-1]; sum(x,y) = ((cells(x,y-1)> 0)&(cells(x,y-1)< t)) + ((cells(x,y+1)> 0)&(cells(x,y+1)< t)) + ... ((cells(x-1, y)> 0)&(cells(x-1, y)< t)) + ((cells(x+1,y)> 0)&(cells(x+1,y)< t)) + ... ((cells(x-1,y-1)> 0)&(cells(x-1,y-1)< t)) + ((cells(x-1,y+1)> 0)&(cells(x-1,y+1)< t)) + ... ((cells(x+1,y-1)> 0)&(cells(x+1,y-1)< t)) + ((cells(x+1,y+1)> 0)&(cells(x+1,y+1)< t)); cells = ((cells==0) & (sum>=t1)) + ... 2*(cells==1) + ... 3*(cells==2) + ... 4*(cells==3) + ... 5*(cells==4) + ... 6*(cells==5) +... 7*(cells==6) +... 8*(cells==7) +... 9*(cells==8) +... 0*(cells==9); 5.森林火灾 规则： 细胞可以在3个不同的状态。状态= 0是空的，state = 1时燃烧，状态=2是森林。 如果一个或更多的4个邻居，如果细胞在燃烧，这是森林（状态= 2），那么新的状态（状态 = 1时）燃烧。 有森林细胞（状态=2），开始自身的燃烧（闪电）（0.000005）是一种低概率。 燃烧（state = 1时）后一个细胞变成空的（状态= 0）。 有一个空单元格成为森林模拟增长是一种低概率（比如说，0.01）。 该阵列被视为toroidly连接的，使燃烧的火焰，以左侧将启动火灾在右边。类似地连接的顶部和底部。 更新代码： sum = (veg(1:n,[n 1:n-1])==1) + (veg(1:n,[2:n 1])==1) + ... (veg([n 1:n-1], 1:n)==1) + (veg([2:n 1],1:n)==1) ; veg = ... 2*(veg==2) - ((veg==2) & (sum> 0 | (rand(n,n)< Plightning))) + ... 2*((veg==0) & rand(n,n)< Pgrowth) ; 请注意的是，环形连接的顺序实施的下标。 6.气体动力学 此CA（及之后的两个）被用来计算粒子的运动。此应用程序需要不同种类的邻里。瓦里斯的小区附近的每一个时间步长，使运动在一定的方向上，将继续在同一方向。换句话说，规则节省势头，这是根据动力计算。所用的邻域被称为一个的马戈利斯附近，由重叠的2×2块的细胞。在下面的图中，左上角的4细胞是邻近地区的偶数时间步长期间，在奇数时间步长，右下4个细胞分数。一个给定的单元有3个邻居在每个时间步长，但具体的细胞构成的邻居来回翻转。 even even even cell odd odd odd 规则： 这就是所谓的HPP-GAS规则。见参考文献[1]第12章.. 细胞有2个状态。状态= 0是空的，状态=1表示运动中的粒子。 在任何一个步骤中，粒子被假设刚才输入的2x2块对角线路径。通过它的中心，因此，它必须跨块。因此，在任何时间步长，每个单元的内容交换与细胞斜对面。对块如下所示。左边的一个是前一段时间后步的权利。六例示，但他们每个人的所有旋转版本相同的方式处理。两种特殊的情况，颗粒与颗粒的碰撞和粒子壁被认为是如下。 0 0 goes 0 0 0 0 to 0 0 1 0 goes 0 0 0 0 to 0 1 1 0 goes 1 0 0 1 to 0 1 1 0 goes 0 1 1 0 to 0 1 1 1 goes 0 1 1 0 to 1 1 1 1 goes 1 1 1 1 to 1 1 为了使粒子碰撞（节省动量和能量），治疗的情况下，如果他们正好一个对角线上的两个粒子击中对方90度偏转。这是通过将一条对角线上的其他时间步长。您可以实现由一个细胞的细胞逆时针旋转。上述第三条规则变为： 1 0 goes 0 1 0 1 to 1 0 为了使粒子碰撞与墙，只需将其状态不变。这导致的一种体现。 更新代码： p=mod(i,2); %margolis neighborhood, where i is the time step %upper left cell update xind = [1+p:2:nx-2+p]; yind = [1+p:2:ny-2+p]; %See if exactly one diagonal is ones %only (at most) one of the following can be true! diag1(xind,yind) = (sand(xind,yind)==1) & (sand(xind+1,yind+1)==1) & ... (sand(xind+1,yind)==0) & (sand(xind,yind+1)==0); diag2(xind,yind) = (sand(xind+1,yind)==1) & (sand(xind,yind+1)==1) & ... (sand(xind,yind)==0) & (sand(xind+1,yind+1)==0); %The diagonals both not occupied by two particles and12(xind,yind) = (diag1(xind,yind)==0) & (diag2(xind,yind)==0); %One diagonal is occupied by two particles or12(xind,yind) = diag1(xind,yind) | diag2(xind,yind); %for every gas particle see if it near the boundary sums(xind,yind) = gnd(xind,yind) | gnd(xind+1,yind) | ... gnd(xind,yind+1) | gnd(xind+1,yind+1) ; % cell layout: % x,y x+1,y % x,y+1 x+1,y+1 %If (no walls) and (diagonals are both not occupied) %then there is no collision, so move opposite cell to current cell %If (no walls) and (only one diagonal is occupied) %then there is a collision so move ccw cell to the current cell %If (a wall) %then don't change the cell (causes a reflection) sandNew(xind,yind) = ... (and12(xind,yind) & ~sums(xind,yind) & sand(xind+1,yind+1)) + ... (or12(xind,yind) & ~sums(xind,yind) & sand(xind,yind+1)) + ... (sums(xind,yind) & sand(xind,yind)); sandNew(xind+1,yind) = ... (and12(xind,yind) & ~sums(xind,yind) & sand(xind,yind+1)) + ... (or12(xind,yind) & ~sums(xind,yind) & sand(xind,yind))+ ... (sums(xind,yind) & sand(xind+1,yind)); sandNew(xind,yind+1) = ... (and12(xind,yind) & ~sums(xind,yind) & sand(xind+1,yind)) + ... (or12(xind,yind) & ~sums(xind,yind) & sand(xind+1,yind+1))+ ... (sums(xind,yind) & sand(xind,yind+1)); sandNew(xind+1,yind+1) = ... (and12(xind,yind) & ~sums(xind,yind) & sand(xind,yind)) + ... (or12(xind,yind) & ~sums(xind,yind) & sand(xind+1,yind))+ ... (sums(xind,yind) & sand(xind+1,yind+1)); sand = sandNew; 7.扩散限制聚集 该系统模拟粘性颗粒聚集形成的分形结构的。粒子的运动发生，例如6 HPP-气体规则类似的规则。所不同的是该颗粒被认为是在某些致密的（但不可见）的液体蹦跳着。其效果是随机的每一个粒子的运动方向在每个时间步长。换句话说，每一个时间步长的碰撞的步骤。仿真也接种与固定在阵列中心的粒子。任何接触它的散射粒子粘到它，本身成为一个不可移动，粘性颗粒。 规则： 使用附近Margolus。在每一个时间步长，旋转4顺时针或逆时针以相等的概率由一个细胞的细胞。旋转随机速度。 移动后，如果一个或多个的八个最接近的neighboors的是一个固定的，粘性颗粒，然后冻结粒子，并使其粘。 更新代码： p=mod(i,2); %margolis neighborhood %upper left cell update xind = [1+p:2:nx-2+p]; yind = [1+p:2:ny-2+p]; %random velocity choice vary = rand(nx,ny)< .5 ; vary1 = 1-vary; %diffusion rule -- margolus neighborhood %rotate the 4 cells to randomize velocity sandNew(xind,yind) = ... vary(xind,yind).*sand(xind+1,yind) + ... %cw vary1(xind,yind).*sand(xind,yind+1) ; %ccw sandNew(xind+1,yind) = ... vary(xind,yind).*sand(xind+1,yind+1) + ... vary1(xind,yind).*sand(xind,yind) ; sandNew(xind,yind+1) = ... vary(xind,yind).*sand(xind,yind) + ... vary1(xind,yind).*sand(xind+1,yind+1) ; sandNew(xind+1,yind+1) = ... vary(xind,yind).*sand(xind,yind+1) + ... vary1(xind,yind).*sand(xind+1,yind) ; sand = sandNew; %for every sand grain see if it near the fixed, sticky cluster sum(2:nx-1,2:ny-1) = gnd(2:nx-1,1:ny-2) + gnd(2:nx-1,3:ny) + ... gnd(1:nx-2, 2:ny-1) + gnd(3:nx,2:ny-1) + ... gnd(1:nx-2,1:ny-2) + gnd(1:nx-2,3:ny) + ... gnd(3:nx,1:ny-2) + gnd(3:nx,3:ny); %add to the cluster gnd = ((sum> 0) & (sand==1)) | gnd ; %and eliminate the moving particle sand(find(gnd==1)) = 0; 8.砂桩 一堆沙子的横截面可以仿照使用一个Margolus的邻里传播细胞，但具有不同的运动规律 规则： 见参考文献[2]第2.2.6章.. 细胞有2个状态。状态= 0是空的，状态=1表示沙agrain。 在任何步骤中，粒子可以落在朝下方的2×2的块。可能的转换如下所示。墙壁和地板停止运动。 为了使议案略有随机，我添加了一个规则，有时反转的两个较低的细胞状态，所有的动作都完成后。 0 0 goes 0 0 0 0 to 0 0 1 0 goes 0 0 0 0 to 1 0 0 1 goes 0 0 0 0 to 0 1 1 0 goes 0 0 1 0 to 1 1 1 0 goes 0 0 0 1 to 1 1 0 1 goes 0 0 1 0 to 1 1 0 1 goes 0 0 0 1 to 1 1 1 1 goes 1 0 1 0 to 1 1 1 1 goes 0 1 0 1 to 1 1 更新代码： p=mod(i,2); %margolis neighborhood sand(nx/2,ny/2) = 1; %add a grain at the top %upper left cell update xind = [1+p:2:nx-2+p]; yind = [1+p:2:ny-2+p]; %randomize the flow -- 10% of the time vary = rand(nx,ny)< .9 ; vary1 = 1-vary; sandNew(xind,yind) = ... gnd(xind,yind).*sand(xind,yind) + ... (1-gnd(xind,yind)).*sand(xind,yind).*sand(xind,yind+1) .* ... (sand(xind+1,yind+1)+(1-sand(xind+1,yind+1)).*sand(xind+1,yind)); sandNew(xind+1,yind) = ... gnd(xind+1,yind).*sand(xind+1,yind) + ... (1-gnd(xind+1,yind)).*sand(xind+1,yind).*sand(xind+1,yind+1) .* ... (sand(xind,yind+1)+(1-sand(xind,yind+1)).*sand(xind,yind)); sandNew(xind,yind+1) = ... sand(xind,yind+1) + ... (1-sand(xind,yind+1)) .* ... ( sand(xind,yind).*(1-gnd(xind,yind)) + ... (1-sand(xind,yind)).*sand(xind+1,yind).*(1-gnd(xind+1,yind)).*sand(xind+1,yind+1)); sandNew(xind+1,yind+1) = ... sand(xind+1,yind+1) + ... (1-sand(xind+1,yind+1)) .* ... ( sand(xind+1,yind).*(1-gnd(xind+1,yind)) + ... (1-sand(xind+1,yind)).*sand(xind,yind).*(1-gnd(xind,yind)).*sand(xind,yind+1)); %scramble the sites to make it look better temp1 = sandNew(xind,yind+1).*vary(xind,yind+1) + ... sandNew(xind+1,yind+1).*vary1(xind,yind+1); temp2 = sandNew(xind+1,yind+1).*vary(xind,yind+1) + ... sandNew(xind,yind+1).*vary1(xind,yind+1); sandNew(xind,yind+1) = temp1; sandNew(xind+1,yind+1) = temp2; sand = sandNew;