元胞自动机--大黄蜂预测

使用元胞自动机借助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