Cellular Automata in Matlab

Cellular Automata in MatlabIntroductionCellular 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 considerationsThe following considerations will be illustrated using a Matlab program which computes Conways life. Parts of the code are explained below. 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 equal axis 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 Matlabs 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 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);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 endExamples1. Conways life.The rule is: o Sum the 8 nearest neighborso If the sum=2 then the state does not change o If the sum=3 then the state=1 o Otherwise the state=0The 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); 2. Surface Tension The rule is:o Sum the 8 nearest neighbors and the cell itselfo If the sum 4 or sum=5 then the state=0 o Otherwise the state=1The 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 0 (at least one neighbor) and the rthreshold, and the cell has never had a neighbor then cell=1. o If the sum 0 set the visited flag for the cell to record that the cell has a nonzero neighbor. 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) ;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. Excitable media (BZ reaction or heart) The rule: o Cells can be in 10 different states. State 0 is resting. States 1-5 are active, and states 6-9 are refractory. o Count the 8 nearest neighbors of each cell which are in one of the active states. o If the sum is greater or equal to 3 (at least three active neighbors) then cell=1. o 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. The update code: x = 2:n-1;y = 2:n-1; sum(x,y) = (cells(x,y-1) 0)&(cells(x,y-1) 0)&(cells(x,y+1) 0)&(cells(x-1, y) 0)&(cells(x+1,y) 0)&(cells(x-1,y-1) 0)&(cells(x-1,y+1) 0)&(cells(x+1,y-1) 0)&(cells(x+1,y+1)=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.5. Forest Fire The rule: o Cells can be in 3 different states. State=0 is empty, state=1 is burning and state=2 is forest. o 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). o There is a low probablity (say 0.000005) of a forest cell (state=2) starting to burn on its own (from lightning). o A cell which is burning (state=1) becomes empty (state=0). o There is a low probability (say, 0.01) of an empty cell becoming forest to simulate growth. o 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. The update code: 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. 6. Gas dynamicsThis 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 evenevencelloddoddodd7. The rule: o This is called a HPP-GAS rule. See reference 1 Chapter 12. o Cells have 2 states. State=0 is empty, state=1 represents a particle in motion. o 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.00goes0000to00o10goes0000to01o10goes1001to01o10goes0110to01o11goes0110to11o11goes1111to11o 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: 10goes0101to10o To make a particle collide with a wall, simply leave its state unchanged. This causes a reflection. The update code: 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 dont 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;8. Diffusion limited aggregation 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.The rule: o 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. o 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. The update code: 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) 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. 9. Sand pile The cross-section of a pile of sand can be modeled using a Margolus neighborhood to propagate cells, but with a different motion ruleThe rule: o See reference 2 Chapter 2.2.6. o Cells have 2 states. State=0 is empty, state=1 represents agrain of sand. o 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. o 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.00goes0000to00o10goes0000to10o01goes0000to01o10goes0010to11o10goes0001to11o01goes0010to11o01goes0001to11o11goes1010to11o11goes0101to11o The update code: o p=mod(i,2); %margolis neighborhoodo sand(nx/2,n