一.思路
一维元胞自动机
一维元胞自动机的意思是,当前元胞的状态,只能被左右元胞的状态影响,而不能被上下元胞影响(能上下就是二维了)
为什么普通一维元胞自动机(只有当前元胞和左右两个邻居)规则有256个?
如上图所示,假如第一个版本有12个元胞(这个数目可以在代码开头设置),为了让每个元胞都有邻居,可以把这12个元胞头尾相连,如上图所示,一个版本就是1行,下一个版本就是下一行。而为了得到下一行的版本,我们需要规则(rule)来确定元胞下次是变白还是变黑。
在上图中,1号元胞的下一个版本是由自身(2)还有左右邻居(a和2)所决定,也就是说是这三个元胞的状态决定了这个元胞(1号元胞)下一个版本的状态,那么规则就可以定8个,2的3次方。另外,每个规则可以决定该规则决定下一时刻该元胞是黑的还是白的,有两种方案。也就是说,最后有256种策略,2的8次方。
这样我们就大致知道,这个一维元胞自动机到底是个什么东西了。
juge函数
用来判断当前三元组是否在rule中是1,是1就返回1,不是1就返回0
这就是juge函数传入的两个参数,左边那个是 1 1 0 ,也就是黑黑白,预迭代的是中间这个1,边上两个是他的邻居
右边那个对应rule30,也就是 0 0 0 1 1 1 1 0 的规则。也就是在这些规则下,我下一步迭代能到1
直观上看,我们知道,左边这样的状态在右边找不到,那么下一步就是迭代到0
具体算法实现思路如下:
通过把左边这个行数组扩大到右边这个矩阵的规模,然后相减,得到结果矩阵A
那么如果我们能在A矩阵中找到一行是 0 0 0的 ,那就能迭代到1,不然就迭代到0
显然右边没有
在程序中,可以先把A中所有元素取绝对值,再转置,原因如下
matlab中有个sum函数,可以得到每一列的和,先给A取绝对值,不然 1 -1 0这样的相加也是0
abs(A)',这里的'是转置的意思,先函数abs(取绝对值函数)发生作用,然后转置
再通过 sum(abs(A)')' 给每一列求和,再转置回来,得到一个列向量
再通过find函数找到这个列向量中为0的下标序列
狗=find(sum(abs(A)')')==0)
然后通过size函数,size(狗,1)就是得到狗的行数,如果狗中有为0的元素,那么就得到一个大于1的数字猫
再通过logical()函数,logical(猫)的意思就是,猫如果是0,返回0,不是0就返回1(负数也算),以达到最后的判断
1.需要确定的参数
lineLength,即第一个版本的元胞数量
rule,即规则,wolfram最钟爱rule30,也就是说rule设置成 0 0 0 1 1 1 1 0
binaryTable, 01234567..的排列方式
注:wolfram钟爱的rule30 设置成0 0 0 1 1 1 1 0 ,顺序是按正常的二进制从右往左,从小到大的顺序
也就对应着76543210这样的顺序,程序内也默认这样的顺序
这里的数字的解释是二进制的解释,比如7就是111,就是3个格子都是黑色
6就是110,就是当前元胞是黑的,左边邻居是黑的,右边邻居是白的
0 - 7每个数组代表一种 3个元胞的组合 然后rule 确定的就是0 1也就是黑白
2.初始元胞序列最中间两个元胞赋1,其余全0
3.进入版本迭代,一个版本算一行
#1.先对预迭代版本的元胞序列进行首尾连接
#2.取得所有的3元组,3个一行的形式保存在矩阵中
#3.用juge函数判断,生成下一个迭代版本,然后用下一个版本生成下下个版本
二.代码+解析
1.第一部分,参数确认,初始化数据
function cellularAutomato
%以下是3个需要调整的内容,包括一行的长度,规则,规则对应的01234567或者76543210的顺序
lineLength=100; //初始元胞序列长度设置成100,可调整
rule=[0 0 0 1 1 1 1 0]; //据说是wolfram最喜欢的rule30,要对应上面 7-0的排列 (可以自己设置256种规则)
binaryTable=[1 1 1;1 1 0;1 0 1;1 0 0;0 1 1;0 1 0;0 0 1;0 0 0]
%binaryTable=[0 0 0;0 0 1;0 1 0;0 1 1;1 0 0;1 0 1;1 1 0;1 1 1] //给出两种0-7的序列排列方式,默认使用上面一种
onePosition=find(rule==1); //通过find函数找到rule中的1的下标
initial=zeros(1,lineLength); //初始化最初的100个元胞,全0
initial([floor(lineLength/2),floor(lineLength/2)+1])=1; //把中间两个元胞设置为1
lenInitial=length(initial); //初始化元胞数组长度
L=floor(lineLength/2); //迭代数目L,元胞总数的一半
currentGeneration=initial; //当前迭代版本
spaceTimeDiagrams=initial; //最终产生灰度图像的矩阵先放入最初的那一行
2.第二部分,进入版本迭代
for i=1:L
currentGeneration=[currentGeneration(end),currentGeneration,currentGeneration(1)]; //首尾相连
currentGeneration
nextGeneration=[];
r=1;
len=2*r+1; //这里由于底下的规则设置都是针对两个邻居和自己,一共三个,所以也可以直接赋3
sequence=[]; //当前行每3个元素为一组存储在这里
for j=1:lenInitial
sequence=[sequence;currentGeneration(j:j+len-1)]; //把三元组按3个一行,100列的方式存储其中
end
for j=1:lenInitial //每个元胞判断下一次的迭代
if judge(sequence(j,:),binaryTable(onePosition,:))
nextGeneration=[nextGeneration,1];
else
nextGeneration=[nextGeneration,0];
end
end
nextGeneration
currentGeneration=nextGeneration;
spaceTimeDiagrams=[spaceTimeDiagrams;nextGeneration]; //把下一次版本添加到矩阵的下一行
end
figure;
subplot(1,2,1),imshow(spaceTimeDiagrams) //subplot(1,2,1)表示把当前窗口分成1行两列,画在第一列
subplot(1,2,2),imshow(~spaceTimeDiagrams) //黑白相反再画一次
end
function value=judge(now_three_element,A) //juge判断函数
mySize=size(A,1);
p=now_three_element;
validate_matrix=ones(mySize,1)*p-A; //矩阵相减
value= logical(size(find(sum(abs(validate_matrix)')'==0),1)); //在上面介绍过
end
三.源码
function cellularAutomato
%以下是3个需要调整的内容,包括一行的长度,规则,规则对应的01234567或者76543210的顺序
lineLength=100;
rule=[0 0 0 1 1 1 1 0]; %据说是wolfram最喜欢的rule30,要对应上面 7-0的排列
binaryTable=[1 1 1;1 1 0;1 0 1;1 0 0;0 1 1;0 1 0;0 0 1;0 0 0]
%binaryTable=[0 0 0;0 0 1;0 1 0;0 1 1;1 0 0;1 0 1;1 1 0;1 1 1]
onePosition=find(rule==1);
initial=zeros(1,lineLength);
initial([floor(lineLength/2),floor(lineLength/2)+1])=1;
lenInitial=length(initial);
L=floor(lineLength/2);
currentGeneration=initial;
spaceTimeDiagrams=initial;
for i=1:L
currentGeneration=[currentGeneration(end),currentGeneration,currentGeneration(1)];
currentGeneration
nextGeneration=[];
r=1;
len=2*r+1; %这里由于底下的规则设置都是针对两个邻居和自己,一共三个,所以也可以直接赋3
sequence=[]; %对当前行每3个元素存储在这里
for j=1:lenInitial
j
sequence=[sequence;currentGeneration(j:j+len-1)];
end
for j=1:lenInitial
if judge(sequence(j,:),binaryTable(onePosition,:))
nextGeneration=[nextGeneration,1];
else
nextGeneration=[nextGeneration,0];
end
end
nextGeneration
currentGeneration=nextGeneration;
spaceTimeDiagrams=[spaceTimeDiagrams;nextGeneration];
end
figure;
subplot(1,2,1),imshow(spaceTimeDiagrams)
subplot(1,2,2),imshow(~spaceTimeDiagrams)
end
function value=judge(now_three_element,A)
mySize=size(A,1);
p=now_three_element;
validate_matrix=ones(mySize,1)*p-A;
value= logical(size(find(sum(abs(validate_matrix)')'==0),1));
end
四.小结
以下是wolfram最钟爱的rule30