[size=+2]Cornell University7 V# i8 w) d1 ?- ^" f
BioNB 441
9 B' o# {9 x0 h% M" V! S4 q6 CCellular Automata in Matlab Introduction Cellular Automata (CA) are a scheme for computing using local rules and local communication. Typcally a CA is defined on a grid, with each point on the grid representing a cell with a finite number of states. A transition rule is applied to each cell simultaneously. Typical transition rules depend on the state of the cell and its (4 or 8) nearest neighbors, although other neighborhoods are used. CAs have applications in parallel computing research, physical simulations, and biological simulations. This page will consider how to write efficient matlab code to implememt CAs and look at some interesting rules.
Matlab code considerations, e8 _1 j3 ~ A' H3 x9 `
The following considerations will be illustrated using a Matlab program which computes Conway's life. Parts of the code are explained below.. v _/ X8 F( U+ {
- Matrixes and images can be converted to one-another, so display is straighforward. If the array cells contains the binary state of each cell and the array z contains zeros, then the cat command builds an RGB image, which is displayed by the image command. The image command also returns a handle. imh = image(cat(3,cells,z,z));set(imh, 'erasemode', 'none')axis equalaxis tight
- Matrixes and images can be converted to one-another, so initial conditions may be computed by matrix or graphics commands. The following code generates an array of zeros, initializes the cell state to zero, then makes a cross of cells in state=1 in the center of the array. One of the examples below (percolation cluster) uses graphics commands to initialize the CA array. z = zeros(n,n);cells = z;cells(n/2,.25*n:.75*n) = 1;cells(.25*n:.75*n,n/2) = 1;
- The Matlab code should be vectorized to minimize overhead. The following two statements compute the sum of nearest neighbors and compute the CA rule. The code makes use of Matlab's very flexible indexing to specify the nearest neighbors. 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);
- Adding a simple GUI is easy. In this example three buttons and a text field were implmented. The three buttons allow a user to Run, Stop, or Quit. The text field displays the number of simulation steps executed. %build the GUI%define the plot buttonplotbutton=uicontrol('style','pushbutton',... 'string','Run', ... 'fontsize',12, ... 'position',[100,400,50,20], ... 'callback', 'run=1;');%define the stop buttonerasebutton=uicontrol('style','pushbutton',... 'string','Stop', ... 'fontsize',12, ... 'position',[200,400,50,20], ... 'callback','freeze=1;');%define the Quit buttonquitbutton=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]); After the controls (and CA) are initialized, the program drops into a loop which tests the state of the variables which are set in the callback functions of each button. For the moment, just look at the nested structure of the while loop and if statements. The program loops until the Quit button is pushed. The other two buttons cause an if statement to execute when pushed. stop= 0; %wait for a quit button pushrun = 0; %wait for a draw freeze = 0; %wait for a freezewhile (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 B6 d2 H2 @" H
Examples
5 S( s8 m, G$ h1 A3 I- Conway's life.
0 X% R9 `5 `2 v, }The rule is:
+ B0 u: f7 Z1 n7 q) B- Sum the 8 nearest neighbors
- If the sum=2 then the state does not change
- If the sum=3 then the state=1
- Otherwise the state=0+ u& ~0 ?* h% D+ w
The code: x = 2:n-1;y = 2:n-1; %nearest neighbor sumsum(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 rulecells = (sum==3) | (sum==2 & cells); - Surface Tension / V/ N# v* D6 V1 i+ H1 B) X
The rule is:/ m1 K# b6 T& n% [& B2 o! @, n; R
- Sum the 8 nearest neighbors and the cell itself
- If the sum< 4 or sum=5 then the state=0
- Otherwise the state=15 H! ]: X/ o1 \/ k8 f2 k
The code: 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)); - Percolation Cluster 9 ]- ?% a( k/ \2 t
The rule:
" S9 d% i5 H; r E$ L- Sum the 8 nearest neighbors of each cell (cells are binary valued, 0/1). Cells also have a separate state variable (called 'visited') which records if they have ever had a nonzero neighbor before.
- Compute a random number r between 0 and 1.
- If the sum> 0 (at least one neighbor) and the r>threshold, and the cell has never had a neighbor then cell=1.
- If the sum> 0 set the "visited" flag for the cell to record that the cell has a nonzero neighbor.
6 p. m( g6 P6 o3 O; u% C) A1 x
The update code:# u8 e& O( i, {1 C1 z% d
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) ; The variables a and b are the size of the image. The inital image is determined by graphics operations. The following statements set up axes of a fixed size, write text into the axes, then grab the contents of the axes and put them back into an array using the getframe function.. 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; After a few dozen time steps (starting with the BioNB 441 image) we get the following image. Click on it to see a full size image.4 A/ x" A1 `3 Q2 \7 K/ z$ C/ y
 - Excitable media (BZ reaction or heart) : d/ _, ^/ k w: ?0 q1 d
The rule:6 x2 i# m- C- J; O, o
- Cells can be in 10 different states. State 0 is resting. States 1-5 are active, and states 6-9 are refractory.
- Count the 8 nearest neighbors of each cell which are in one of the active states.
- If the sum is greater or equal to 3 (at least three active neighbors) then cell=1.
- States 1 to 9 occur stepwise with no more input. If state=1 then the next state=2. If state=2 then the next state=3, and similarly of all the states up to 9. If state=9 then the next state=0 and the cell is back at rest.
9 H/ M6 {7 S; g" I) z
The update code: 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); An image from this CA after spiral waves develop is shown below.
; K' [, C$ F9 Y+ R; X+ j% x - Forest Fire
& m. w4 k) I% |2 Y: [ q* {5 pThe rule:
" N1 Z: R1 O2 p) M/ O- Cells can be in 3 different states. State=0 is empty, state=1 is burning and state=2 is forest.
- If one or more of the 4 neighbors if a cell is burning and it is forest (state=2) then the new state is burning (state=1).
- There is a low probablity (say 0.000005) of a forest cell (state=2) starting to burn on its own (from lightning).
- A cell which is burning (state=1) becomes empty (state=0).
- There is a low probability (say, 0.01) of an empty cell becoming forest to simulate growth.
- The array is considered to be toroidly connected, so that fire which burns to left side will start fires on the right. The top and bottom are similarly connected.
! q/ }/ P, {2 R
The update code: h% S+ d4 x2 ~1 J. P+ s h( r
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) ; Note that the toroidal connection is implemented by the ordering of subscripts. - Gas dynamics. w* x6 Z: X* p/ I, {0 x
This CA (and the next two) are used to compute motion of particles. This application requires a different kind of neighborhood. The neighborhood of a cell varys on each time step so that motion in a certain direction will continue in the same direction. In other words, the rule conserves momentum, which is the basis of dynamical calculations. The neighborhood used is called a Margolis neighborhood and consists of overlapping 2x2 blocks of cells. In the following diagram, the upper-left 4 cells are the neighborhood during an even time-step, the bottom-right 4 cells are the neighborhood during an odd time-step. A given cell has 3 neighbors at each time-step, but the specific cells which constitute the neighbors flips back and forth. | even | even | | | even | cell | odd | | | odd | odd |
The rule:
/ }& Y- J# [1 W+ E- This is called a HPP-GAS rule. See reference [1] Chapter 12..
- Cells have 2 states. State=0 is empty, state=1 represents a particle in motion.
- On any step, a particle is assumed to have just entered the 2x2 block on a diagonal path. It must therefore cross the block through its center. So on any time-step, swap the contents of each cell with that of the cell diagonally across from it. Pairs of blocks are shown below. The left one is before and the right after a time-step. Six cases are shown, but all rotational versions of each of them are treated the same way. Two special cases, particle-particle collision and particle-wall are considered below.% P; G7 K$ O* f; {6 O4 D) R
" v0 a: ^7 `3 M' `9 w/ V
- To make a particle collision (which conserves momentum and energy), treat the case of exactly two particles on a diagonal as if they hit and deflected each other 90 degrees. This is done by converting one diagonal to the other on the time-step. You can implement this by rotating the 4 cells counterclockwise by one cell. The third rule above becomes:
]5 b& C) H1 P( o, v% o - To make a particle collide with a wall, simply leave its state unchanged. This causes a reflection.
$ k9 Y8 e/ Z: y& D. L9 e/ [" L, d5 F, {- v5 \; p- ]* P) `; Y
The update code: - d, E0 Z$ |: Q
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; - Diffusion limited aggregation 8 h2 w$ y2 U! o8 r( ^, \
This system simulates sticky particles aggregating to form a fractal structure. particle motion takes place with a rule similar to the HPP-GAS rule in example 6. The difference is that particles are assumed to be bouncing around in some dense (but invisible) liquid. The effect is to randomize the direction of motion of every particle at every time step. Put differently, every time-step is a collision step. The simulation is also seeded with one fixed particle in the center of the array. Any diffusing particle which touches it sticks to it, and itself becomes a non-moving, sticky particle.
" s* P7 D6 t( V- W' J$ D; ?; cThe rule:3 P5 m7 @4 I% ^5 S! J( o" Y1 w ^1 p
- Use a Margolus neighborhood. At every time step, rotate the 4 cells either clockwise or counterclockwise by one cell with equal probability. The rotation randomizes the velocities.
- After the move, if one or more of the eight nearest neighboors is a fixed, sticky particle, then freeze the particle and make it sticky.: y% T% r" ]2 l7 I* P
The update code:
, j/ j# Y- c# Y2 n. ^; k0 ?' | 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; The following image shows the fixed cluster after many time steps.  - Sand pile
7 r2 B( ?: a( ?$ f( u1 P- o- YThe cross-section of a pile of sand can be modeled using a Margolus neighborhood to propagate cells, but with a different motion rule+ c3 s& h1 U6 v
The rule:4 k4 U7 P; @; e4 V: K! x
- See reference [2] Chapter 2.2.6..
- Cells have 2 states. State=0 is empty, state=1 represents agrain of sand.
- On any step, a particle can fall toward the bottom of the the 2x2 block. The possible transitions are shown below. Walls and floors stop motion.
- To make the motion slightly random, I added a rule which sometimes reverses the state of the two lower cells after all the moves are complete." _; V6 K5 X8 [% ]1 o0 f
The update code:
) O' x; A$ Y9 l4 h7 o: F 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; & y3 D* t! X3 I* @, e4 F$ \6 T/ W3 z
References
8 E$ M0 G# l; k7 W* v[1] Cellular Automata Machines by Tommaso Toffoli and Norman Margolus, MIT Press, 1987.
5 M0 S( {" b& b) p: d[2] Cellular Automata Modeling of Physical Systems by Bastien Chopard and Michel Droz, Cambridge University Press, 1998." |; c5 [2 x) [# I# u: z
$ g; Z8 p+ d* `3 Q, o& S- j
https://courses.cit.cornell.edu/bionb441/CA/+ Y9 H: Y8 ]1 @7 n9 }$ I
|
发表于 2014-8-5 12:33
自动把其余的细胞变成黑色或保留白色。上左图就是根据规则110运行了前20步的情况，在这里似乎看不出什么有趣的东西。但运行到几百步后，就出现了一些有趣的特征，一些结构开始既不是周期性地也不是完全随机地出现在画面上。上右图就是按规则110运行到700步的情况。1984年沃尔夫勒姆把256种规则分成了四类：第一类只生成简单重复的图案，比如全黑、全白、或黑白相间如国际象棋棋盘等等；第二类规则产生一些自相似的分形图案，形成稳定的嵌套结构；第三类规则产生的图案具有明显的随机性；第四类规则产生复杂的图案，但不简单重复。这些图案既不是规则的也不是完全随机的。它们呈现出某种有序性，但却不能被预言。 
errier，J.Y.,1996)。同时。它为动力学系统理论中有关秩序 (Ordering)、紊动 (Turbulence)、混沌 (Chaos)、非对称(Symmetry-Breaking)、分形(Fractality)等系统整体行为与复杂现象的研究提供了一个有效的模型工具 (Vichhac。G，1984; Bennett，C,1985)。 