cellular automaton

Cellular automata (cellular automata,CA) It's a kind of time 、 Space 、 All States are discrete , Spatial interaction and temporal causality are local grid dynamics models , It has the ability to simulate the spatio-temporal evolution of complex systems .

In fact, I tried to teach myself when I was preparing for the national games last summer vacation , But limited by , There are few easy to understand and efficient documents or video teaching on the Internet , And personal learning ability and time are limited , I was confused at that time , So I gave up studying at that time .

During the recent vacation preparing for the mathematical modeling competition , Modeling teammates let me learn about cellular automata simulation . Spent some time studying and thinking , Sort out some things .

In my opinion, cellular automata , The most puzzling thing is ：

1. How to set the specific color of cells

2. How to apply cellular automata to specific areas （ Such as the realization of forest fire in Australia , That is, import the green plant data from Australia ）

3. How to adapt the state matrix to more complex rules

Recently, the author watched a large number of teaching videos and papers , According to the above questions, I have done some sorting , This process makes me understand the running process of cellular automata again , And a set of programming framework suitable for most simulation processes is summarized . Here we will not repeat the basic knowledge of cellular automata , Let's talk about our opinions directly in combination with several specific examples of the above three problems . （ All of them are used in this paper MATLAB To achieve ）

1. How to set the specific color of cells

First , In my personal understanding , The process of cell changing with time is essentially the transformation of matrix . The values in the state matrix are taken from a set of finite elements （ For example, in forest fire cellular automata , We make the forest 0, Land is 1, The ignition point is 3）. And cellular automata “ Color ” It is just for the convenience of us to observe the change law “ visualization ” The process , We make the forest （0） The position of is green , land （1） The position of is yellow , Ignition point （3） The position of is red .

In color image , The color of each point is affected by three primary color vectors [R,G,B] To control ,R,G,B The value of is floating on 0 To 1 Within the range of . And the two extremes of color —— Black and white RGB The vector is [0,0,0] and [1,1,1]. in other words , In order to make the value of a certain size in the matrix show the color we need , We need to set a specific RGB vector . It's not difficult to understand. , If one m n Each point of the size matrix shows the color we want to see , We need to give a mn3 To control its color , namely R,G,B From a vector to a matrix .

                

Here we can use cat Function to combine its RGB matrix

When encountering a concrete implementation , I see a lot B standing UP Or the document is directly transformed in the final loop RGB matrix , Here we can choose to write a color conversion function , Input state matrix , Turn it into , In this way, the state matrix can be directly transformed in each cycle , Let each function of the program have a clear division of labor , It is convenient for us to better debug the logic of the program . The following is the color transformation function of forest fire ：

function [R,G,B] = transformcolor(b1)
[m,n]=size(b1);    
[R,G]=deal(ones(m,n));
B=zeros(m,n);
R(b1==0)=0;
B(b1==2)=1;
G(b1==3)=0;
end

[R,G,B]=transformcolor(b1);
ch=imagesc(cat(3,R,G,B));

cat The function will have three R,G,B The matrix is synthesized into a high-dimensional matrix , It is convenient for us to visualize the current state . 

MATLAB Graphic R,G,B The color comparison table can be found in the following links ：

https://blog.csdn.net/ha_____ha/article/details/103683988

2. How to apply cellular automata to specific areas

  Actually , Importing environment specific data can be a simple process of image processing . Take forest fires for example , We can find the vegetation distribution map of each country on the Internet , The author looked for the national vegetation map of the United States during his own practice , as follows ：

A simple way to do this is to gray the color map , About to image （mn3）  The 3 d RGB Convert to matrix （mn） Two dimensional matrix of .（ As I said before RGB The matrix is slightly different , When importing an image RGB Values are 0-255 The numerical , Not at all 0-1）

The author chose the red R matrix （ selection RGB The one with the most obvious color difference in the matrix ）, use imshow The observations are as follows

here , We can select a threshold according to the gray image ： The values of the parts circled in green are mostly distributed in 0-160 Between （ The forest part ）, The values of the continental part pointed by the yellow arrow are mostly distributed in 160-245 Between , The rest are distributed in 245-255 Between . The blank space does not belong to the territory of the United States , Convenient for post transformation processing , The author set the part outside the national boundary as 2.

Here the entire matrix is iterated using a loop Zoning ：

a=imread("E:\USEfulEEE\Mathmodel\CA\China\America.jpg"); % Read images 
b3=a(:,:,3);
b2=a(:,:,2);
b1=a(:,:,1);      % Observed separately RGB matrix , Select the one with the most obvious value difference 
imshow(b1)  %imshow(b2);imshow(b3)
[m,n]=size(b1);
for i=1:m
    for j=1:n
        if b1(i,j)<=160
            b1(i,j)=0;   % green （ The forest ） by 0
        elseif b1(i,j)>160&&b1(i,j)<245
            b1(i,j)=1;   % land （ An area that can grow into a forest ） by 1
        else
            b1(i,j)=2;   % The boundary is 2
        end
    end
end
[R,G,B]=transformcolor(b1);   % Color transformation of state matrix 
ch0=imagesc(cat(3,R,G,B));    % A map of the United States vegetation divided by visual areas 

  The matrix divided by regions is the state matrix that we use for transformation in the later stage , This process is equivalent to the initialization of the state .

 3. How to adapt the state matrix to more complex rules

The process of thinking about the computer realization of state transformation is a complicated logical sorting process . Here, the author also uses a single encapsulated function to realize each transformation process .

The change rule of forest fire cellular machine is ：

1. Around the forest 8 Neighbors, if any 1 Is on fire , The cell has a high probability of fire in the next state

2. The fire state will go out in the next state , Into land

3. Each forest cell has a very small probability of being struck by lightning , And start a fire

4. Each land cell has a very small probability of growing a new forest

5. Cells outside the U. S. territory are not affected by any cells , Keep yourself in shape

We are at the beginning of every new state , The state transition function captures the forest of old States 、 land 、 To catch fire 、 Status outside the territory

function new = generatenew(old)
[m,n]=size(old);
treeold=old==0;            %old The location of the forest 
landold=old==1;            %old The location of the middle land 
fireold=old==3;            %old Where the fire broke out 
lightning=0.000002;        % Probability of fire 
growth=0.001;              % Probability of new tree growth 

With the old state , We need to generate new states based on rules and old states , That is, which locations in the new state are forests , To catch fire , Beyond land and territory .

The author lists several important relations according to the five rules ：

1. The location of the new state forest = The location of the old forest - Where the fire broke out in the new state + Where new forests are growing

2. The location of the new state land = The location of land in the old state - Where new forests are growing + Where the fire broke out in the old state

3. Where the fire broke out in the new state = The position in the lightning strike + The fire location associated with the fire in the old state

It's not hard to see. , The total amount of several cells in the old state and the total amount of several cells in the new state are always conserved .

countfire=zeros(m,n);
countfire(2:m-1,2:n-1)=fireold(1:m-2,2:n-1)+fireold(3:m,2:n-1)+fireold(2:m-1,1:n-2)+fireold(2:m-1,3:n);

[r,r1]=deal(rand(m,n));
grow=r.*landold;                        %grow It has value where it can grow 
light=r.*treeold;                       %light It's worth it where there might be a fire 
liandai=(countfire>=1).*r1;
liandai=(liandai>=0.35)&treeold;

grownew=(grow>=1-growth);               % The position of the newly grown tree is 1, Others are 0 
firenew=(light>=1-lightning)|liandai;   % The new ignition point is the position in the lightning strike + The location of the old ignition point 
land=landold-grownew+fireold;           % The location of the land is the old land + Where the flame burned out 
tree=treeold+grownew -firenew;

new=2*ones(m,n);                        % New initialization is all about borders 
new(tree==1)=0;
new(land==1)=1;
new(firenew==1)=3;
end

Realized a state transition , We add a in the main program for Loop or while loop , To observe the continuous change process of the fire system .

PS： Pay attention to parallelism , We cannot use loops to change the state of cells in a cell space one by one , This will cause the previous changes to affect the subsequent operation . And the use of MATLAB The advantages of the matrix can be changed to ensure the speed .（ I've also tried to use loops to change each cell in turn , The time cost of running like this is really very high ）

hold on;
axis equal
axis tight
axis off
for i=1:100
    b1=generatenew(b1);
    [R,G,B]=transformcolor(b1);
    ch0=imagesc(cat(3,R,G,B));
    pause(0.0001)
end

The change process is demonstrated as follows （ conversion gif After that, it was a little vague ）：

Example 2 of cellular automata —— Desertification protection forest

The author is in bilibili See a more interesting cellular automata simulation on , About the density of the shelter forest 、 Study on the relationship between width and desertification prevention capacity , That is to simulate the protective effect of the shelter forest .UP The LORD made such a rule description for it ：

Links are as follows ：https://www.bilibili.com/video/BV1nK4y1p7BD

  Follow the template mentioned before , I programmed .

First initialize the area , One 300300 Region ,100-150 The area is shelterbelt , On its left is the desert （ Not to be used by others ）, On the right is the open space （ Can be used by people ）, For the convenience of handling , We also set the boundaries of cellular automata .

Here we make the open space 0, Desert is 1, Forest as 2, The boundary is 3.

The planting density of the forest area can be changed , With variable percentage To express . It's not hard to see. , The larger the variable , The denser the green part .

Encapsulate the following color matching function ：

function [R,G,B] = transformcolor(b1)
[m,n]=size(b1);    
[R,G]=deal(ones(m,n));
B=zeros(m,n);
R(b1==0)=0;
B(b1==2)=1;
G(b1==3)=0;
end

Then according to the rules of state change , List the following formula ：

1. The location of the new state desert = The location of the old desert - The location of the desert to the open space + The open space turns to the desert

2. Location of the new state vacant lot = The location of the old vacant lot - The open space turns to the desert + The location of the desert to the open space + The location of the forest to the open space

3. The location of the new state forest = The location of the old forest - The location of the forest to the open space

The state transition function is as follows ：

function outlinenew = transferarray(outlineold,m,n)
spaceold=outlineold==0;
desertold=outlineold==1;
forestold=outlineold==2;
[count1,count2]=deal(zeros(m,n));
c1=desertold(1:m-2,2:n-1)+desertold(3:m,2:n-1)+desertold(2:m-1,1:n-2)+desertold(2:m-1,3:n);
count1(2:m-1,2:n-1)=c1;     %count1 For the number of deserts around each point 
c2=forestold(1:m-2,2:n-1)+forestold(3:m,2:n-1)+forestold(2:m-1,1:n-2)+forestold(2:m-1,3:n);
count2(2:m-1,2:n-1)=c2;     %count2 For the number of forests around each point 
r=rand(m,n);
r1=r>0.7;   
%%%%% Open space turns into desert %%%%%%
desertadd=spaceold & (count1>=1) & r1 ;
%%%%% Desert to open space %%%%%%
desertreduce=desertold & (count2>=3) &r1;
%%%%% Forest turns to open space %%%%%%
spacenew=forestold & (count1>=3) & r1;
%%%%% Overall change process %%%%%%
outlinenew=3*ones(m,n);
outlinenew(desertold+desertadd-desertreduce==1)=1;
outlinenew(forestold-spacenew==1)=2;
outlinenew(spaceold+spacenew-desertadd+desertreduce==1)=0;
end

The main function is as follows ：

clear;clc;
percentage=0.45;    % Adjustable forest density 
in=zeros(300);
in(:,1:99)=1;
r=rand(300,51);
r(r>=1-percentage)=2;
r(r<1-percentage)=0;
in(:,100:150)=r;
outline=3*ones(302);
outline(2:301,2:301)=in;

[R,G,B]=transformcolor(outline);
ch=imagesc(cat(3,R,G,B));
[m,n]=size(outline);
hold on
for i=1:300
    outline=transferarray(outline,m,n);
    [R,G,B]=transformcolor(outline);
    ch=imagesc(cat(3,R,G,B));
    pause(0.08)
end

Running effect ：

1. The forest density is less than 50%（ Set to 35%）：

It can be seen that , The desert quickly broke through the shelter forest , Eroded the available land of human beings .

 2. The forest density is greater than 50%（ Set to 65%）：

It can be seen that , The protective effect of the shelter forest is very good .

This passage records how an idle man spent his boring Spring Festival , In the process of thinking, I realized a lot , It also brought me a little sense of achievement and self-confidence .

Strongly push the cellular automata tutorial of the lecture on large mathematical modeling ：

https://www.bilibili.com/video/BV1fV411q7SA?from=search&seid=6908856970843696115&spm_id_from=333.337.0.0

Cellular automata （ Simulated maple leaf corrosion ）_ Bili, Bili _bilibili

Keep going .