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

(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;