使用元胞自动机借助MATLAB预测大黄蜂种群发展趋势
clear;clc;
dbstop if error
doc = '2021MCMProblemC_DataSet.xlsx';
data = readtable(doc,'VariableNamingRule','preserve');
clear doc;
shp = shaperead('usastatehi.shp');
n = size(data,1);
posid = [];
lat = zeros(n,1);
lon = zeros(n,1);
for i = 1:n
if strcmp(cell2mat(data.("Lab Status")(i)),'Positive ID')
posid(end+1) = i;
end
lat(i) = data.Latitude(i);
lon(i) = data.Longitude(i);
end
%% 创建地图
M.ylim = [min(lat),max(lat)];
M.xlim = [min(lon),max(lon)];
M.pos = [lat(posid),lon(posid)];
map = creatworldmap(M);
%% 抽象出对应的网格
G.size = 100;%地图被划分为的栅格数
fakemap = zeros(G.size);
temprowid = linspace(M.ylim(1),M.ylim(2),G.size);
tempcolid = linspace(M.xlim(1),M.xlim(2),G.size);
%找到正类所在坐标索引
latid = zeros(size(posid,2),1);
lonid = zeros(size(posid,2),1);
for i = 1:size(latid,1)
temp = 1;
while lat(posid(i))>temprowid(temp)
temp = temp+1;
end
latid(i) = temp;
end
for i = 1:size(lonid,1)
temp = 1;
while lon(posid(i))>tempcolid(temp)
temp = temp+1;
end
lonid(i) = temp;
end
ind = [latid,lonid];
clear latid lonid temp
% x = ind(:,2);
% y = 100-ind(:,1);
% fakemap(x,y)=1;
figure(2)
map2 = imshow(fakemap);
for i = 1:size(ind,1)
fakemap(G.size-ind(i,1),ind(i,2)) = 1;
map2.CData(G.size-ind(i,1),ind(i,2))=1;
end
%% 开始迭代
ITER = 500;
%计算每个格子的概率
point = [10,10;20,20;30,30;40,40];
probmap = getprob(fakemap,point);
% probmap = zeros(G.size);
deadperid = 3;
r = cell(deadperid,1);
c = cell(deadperid,1);
for Times = 1:ITER
if Times <=deadperid
if Times == deadperid
[r{deadperid},c{deadperid}] = find(map2.CData==1);
else
[r{mod(Times,deadperid)},c{mod(Times,deadperid)}] = find(map2.CData==1);
end
else
r(1:deadperid-1) = r(2:deadperid);c(1:deadperid-1) = c(2:deadperid);
[r{end},c{end}] = find(map2.CData==1);
end
%密度检测>4则死亡
mapexd = padarray(map2.CData,[1,1]);
nummap = mapexd(1:end-2,1:end-2)+mapexd(1:end-2,2:end-1)+mapexd(1:end-2,3:end)+...
mapexd(2:end-1,1:end-2)+mapexd(2:end-1,3:end)+...
mapexd(3:end,1:end-2)+mapexd(3:end,2:end-1)+mapexd(3:end,3:end);
map2.CData(nummap>4)=0;
clear mapexd nummap
%周期检测>?死亡
if Times>deadperid
deadgroup = [r{1},c{1}];
map2.CData(deadgroup(:,1),deadgroup(:,2)) = 0;
end
clear deadgroup
%依概率选择下一个要转移到哪里去
livegroup = [r{deadperid},c{deadperid}];
for i = 1:size(livegroup,1)
p = livegroup(i,:);
mapexd = padarray(probmap,[1,1]);
neighbors = [p(1)-1,p(2)-1;p(1),p(2)-1;p(1)+1,p(2)-1;p(1),p(2)-1;p(1),p(2)+1;p(1)+1,p(2)-1;...
p(1)+1,p(2);p(1)+1,p(2)+1];
%剔除格子上已有蜂群的
for j = 8:-1:1
if neighbors(j,1)>G.size || neighbors(j,2)>G.size || neighbors(j,1)<1 || neighbors(j,2)<1
neighbors(j,:) = [];
end
end
for j = size(neighbors,1):-1:1
if map2.CData(neighbors(j,1),neighbors(j,2))==1
neighbors(j,:) = [];
end
end
if ~isempty(neighbors)
neigh = getneigh(neighbors,mapexd);
map2.CData(neigh(1),neigh(2))=1;
end
end
% pause(0.5);
%在map上改变热力图
[rowid,colid] = find(map2.CData==1);
m = length(rowid);
templat = zeros(m,1);
templon = zeros(m,1);
for i = 1:m
if rowid(i)==100 || colid(i) == 100
continue
else
templat(i) = temprowid(G.size-rowid(i))+(rand(1)*2-1)/8000;
templon(i) = tempcolid(colid(i))+(rand(1)*2-1)/8000;
end
end
weight = ones(size(templat,1),1);
map.LatitudeData = templat';
map.LongitudeData = templon';
map.WeightData = weight';
pause(0.1)
disp(['第',num2str(Times),'次迭代'])
end
%%
function map = creatworldmap(M)
weights = ones(size(M.pos,1),1);
figure(1)
map = geodensityplot(M.pos(:,1),M.pos(:,2),weights);
geolimits([M.ylim],[M.xlim])
end
function map = getprob(fakemap,point)
map = fakemap;
for i = 1:size(map,1)
for j = 1:size(map,2)
for k = 1:size(point,1)
map(i,j) = map(i,j)+exp(-((i-point(k,1))^2+(j-point(k,2))^2)/1000);
end
end
end
end
function neigh = getneigh(neighbors,mapexd)
corprob = zeros(size(neighbors,1),1);
for i = 1:size(corprob,1)
corprob(i) = mapexd(neighbors(i,1),neighbors(i,2));
end
corprob = corprob./sum(corprob);
corprob = cumsum(corprob);
a = rand(1);
t = 1;
while a>corprob(t)
t = t+1;
end
neigh = [neighbors(t,1),neighbors(t,2)];
end