📄 ca.htm
字号:
<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.0 Transitional//EN">
<!-- saved from url=(0053)http://instruct1.cit.cornell.edu/courses/bionb441/CA/ -->
<HTML><HEAD><TITLE>CA</TITLE>
<META http-equiv=Content-Type content="text/html; charset=iso-8859-1">
<META content="MSHTML 6.00.2900.2180" name=GENERATOR></HEAD>
<BODY bgColor=#ffffff>
<P align=center><FONT size=+2>Cornell University<BR>BioNB 441<BR>Cellular
Automata in Matlab</FONT></P>
<P align=left><B>Introduction</B></P>
<P align=left>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. </P>
<HR>
<P><B>Matlab code considerations</B></P>
<P>The following considerations will be illustrated using a <A
href="http://instruct1.cit.cornell.edu/courses/bionb441/CA/calife.m">Matlab
program</A> which computes Conway's life. Parts of the code are explained
below.</P>
<UL>
<LI>Matrixes and images can be converted to one-another, so display is
straighforward. If the array <CODE>cells</CODE> contains the binary state of
each cell and the array <CODE>z</CODE> contains zeros, then the
<CODE>cat</CODE> command builds an RGB image, which is displayed by the
<CODE>image</CODE> command. The <CODE>image</CODE> command also returns a
handle. <PRE>imh = image(cat(3,cells,z,z));
set(imh, 'erasemode', 'none')
axis equal
axis tight
</PRE>
<LI>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. <PRE>z = zeros(n,n);
cells = z;
cells(n/2,.25*n:.75*n) = 1;
cells(.25*n:.75*n,n/2) = 1;
</PRE>
<LI>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. <PRE>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);
</PRE>
<LI>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. <PRE>%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]);
</PRE>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. <PRE>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
</PRE></LI></UL>
<HR>
<P><B>Examples</B></P>
<OL>
<LI><A
href="http://instruct1.cit.cornell.edu/courses/bionb441/CA/calife.m">Conway's
life</A>.<BR>The rule is:
<UL>
<LI>Sum the 8 nearest neighbors<BR>
<LI>If the sum=2 then the state does not change
<LI>If the sum=3 then the state=1
<LI>Otherwise the state=0<BR></LI></UL>The code: <PRE>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);
</PRE>
<LI><A
href="http://instruct1.cit.cornell.edu/courses/bionb441/CA/CAsurfaceT.m">Surface
Tension</A> <BR>The rule is:<BR>
<UL>
<LI>Sum the 8 nearest neighbors and the cell itself<BR>
<LI>If the sum< 4 or sum=5 then the state=0
<LI>Otherwise the state=1<BR></LI></UL>The code: <PRE>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));
</PRE>
<LI><A
href="http://instruct1.cit.cornell.edu/courses/bionb441/CA/CAvine.m">Percolation
Cluster</A> <BR>The rule:
<UL>
<LI>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.
<LI>Compute a random number r between 0 and 1.
<LI>If the sum> 0 (at least one neighbor) and the r>threshold, and the
cell has never had a neighbor then cell=1.
<LI>If the sum> 0 set the "visited" flag for the cell to record that the
cell has a nonzero neighbor. </LI></UL>The update code:<BR><PRE> 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) ;
</PRE>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 <CODE>getframe</CODE> function.. <PRE>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;
</PRE>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.<BR><A
href="http://instruct1.cit.cornell.edu/courses/bionb441/CA/percolation.png"><IMG
height=225 src="CA.files/percolationSmall.png" width=300 border=0></A>
<P></P>
<LI><A
href="http://instruct1.cit.cornell.edu/courses/bionb441/CA/excitable.m">Excitable
media</A> (BZ reaction or heart) <BR>The rule:
<UL>
<LI>Cells can be in 10 different states. State 0 is resting. States 1-5 are
active, and states 6-9 are refractory.
<LI>Count the 8 nearest neighbors of each cell which are in one of the
active states.
<LI>If the sum is greater or equal to 3 (at least three active neighbors)
then cell=1.
<LI>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. </LI></UL>The update code: <PRE>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);
</PRE>An image from this CA after spiral waves develop is shown
below.<BR><IMG height=251 src="CA.files/excitable1.png" width=250>
<P></P>
<LI><A
href="http://instruct1.cit.cornell.edu/courses/bionb441/CA/forest.m">Forest
Fire</A> <BR>The rule:
<UL>
<LI>Cells can be in 3 different states. State=0 is empty, state=1 is burning
and state=2 is forest.
<LI>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).
<LI>There is a low probablity (say 0.000005) of a forest cell (state=2)
starting to burn on its own (from lightning).
<LI>A cell which is burning (state=1) becomes empty (state=0).
<LI>There is a low probability (say, 0.01) of an empty cell becoming forest
to simulate growth.
<LI>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. </LI></UL>The update code: <BR><PRE>
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) ;
</PRE>Note that the toroidal connection is implemented by the ordering of
subscripts.
<P></P>
<LI><A href="http://instruct1.cit.cornell.edu/courses/bionb441/CA/gas2.m">Gas
dynamics</A><BR>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.
<P>
<TABLE height=116 width="17%" border=2>
<TBODY>
<TR>
<TD width="33%">even </TD>
<TD width="33%">even</TD>
<TD width="34%"> </TD></TR>
<TR>
<TD width="33%">even</TD>
<TD width="33%">cell</TD>
<TD width="34%">odd</TD></TR>
<TR>
<TD width="33%"> </TD>
<TD width="33%">odd</TD>
<TD width="34%">odd</TD></TR></TBODY></TABLE>
<P>The rule: <BR>
<UL>
<LI>This is called a HPP-GAS rule. See reference [1] Chapter 12..
<LI>Cells have 2 states. State=0 is empty, state=1 represents a particle in
motion.
<LI>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.<BR>
<TABLE width="23%" border=2>
<TBODY>
<TR>
<TD width="18%">0</TD>
<TD width="18%">0</TD>
<TD width="32%">goes</TD>
<TD width="15%">0</TD>
<TD width="17%">0</TD></TR>
<TR>
<TD width="18%">0</TD>
<TD width="18%">0</TD>
<TD width="32%">to</TD>
<TD width="15%">0</TD>
<TD width="17%">0</TD></TR></TBODY></TABLE>
<P>
<TABLE width="23%" border=2>
<TBODY>
<TR>
<TD width="18%">1</TD>
<TD width="18%">0</TD>
<TD width="32%">goes</TD>
<TD width="15%">0</TD>
<TD width="17%">0</TD></TR>
<TR>
<TD width="18%">0</TD>
<TD width="18%">0</TD>
<TD width="32%">to</TD>
<TD width="15%">0</TD>
<TD width="17%">1</TD></TR></TBODY></TABLE>
<P>
<TABLE width="23%" border=2>
<TBODY>
<TR>
<TD width="18%">1</TD>
<TD width="18%">0</TD>
<TD width="32%">goes</TD>
<TD width="15%">1</TD>
<TD width="17%">0</TD></TR>
<TR>
<TD width="18%">0</TD>
<TD width="18%">1</TD>
<TD width="32%">to</TD>
<TD width="15%">0</TD>
<TD width="17%">1</TD></TR></TBODY></TABLE>
<P>
<TABLE width="23%" border=2>
<TBODY>
<TR>
<TD width="18%">1</TD>
<TD width="18%">0</TD>
<TD width="32%">goes</TD>
<TD width="15%">0</TD>
<TD width="17%">1</TD></TR>
<TR>
<TD width="18%">1</TD>
<TD width="18%">0</TD>
<TD width="32%">to</TD>
<TD width="15%">0</TD>
<TD width="17%">1</TD></TR></TBODY></TABLE>
<P>
<TABLE width="23%" border=2>
<TBODY>
<TR>
<TD width="18%">1</TD>
<TD width="18%">1</TD>
<TD width="32%">goes</TD>
<TD width="15%">0</TD>
<TD width="17%">1</TD></TR>
<TR>
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -