SIRS.m

% SIRS.m
% converted from diffusionmat on 1/4/16
% Susecpt-Infect-Remove-Suscept  plus an innoculation of highest-degree node

clear
close all

T1 = 100;
T2 = 200;

showfig = 1;

Innoc = 1;  % Turn innoculation on

% beta is infection rate
% mu is recovery rate

N = 50;  m = 2;    beta = 0.2;  mu = 0.65;               node = makeSF(N,m);
%N = 50;  m = 2;    beta = 0.2;  mu = 0.4;   p = 0.1;   node = makeSW(N,m,p);   
%N = 50;  p = 0.06; beta = 0.2;  mu = 0.4;               node = makeER(N,p);      

[N,e,avgdegree,maxdegree,mindegree,numclus,meanclus,Lmax,L2,LmaxL2] = clusterstats(node);

disp(' ')
displine('Number of nodes = ',N)
disp(strcat('Number of edges = ',num2str(e)))
disp(strcat('Mean degree = ',num2str(avgdegree)))
displine('Maximum degree = ',maxdegree)
disp(strcat('Number of clusters = ',num2str(numclus)))
disp(strcat('mean cluster coefficient = ',num2str(meanclus)))
disp(' ')
disp(strcat('Lmax = ',num2str(Lmax)))
disp(strcat('L2 = ',num2str(L2)))
disp(strcat('Lmax/L2 = ',num2str(LmaxL2)))
disp(' ')
disp(strcat('fac =',num2str(mu/avgdegree/beta)))
disp(' ')

[A,degree,Lap] = adjacency(node);

[V,D] = eig(Lap);

for loop = 1:N
    eigval(loop) = D(loop,loop);
end

figure(1)
plot(eigval)
title('Eigenvalues')

% initial values
a = zeros(N,1);
[Y,I] = max(degree);
a(I) = 1;

fh2 = figure(2);
drawnet(node)

% The discrete-time approach

c0 = a;
dt = 1;

R = eye(N,N);
Ct = zeros(1,N);
for nloop = 1:N
    deg(nloop) = node(nloop).numlink;
end

tic
c = c0;
fh3 = figure(3);
eps = 1e-6;
Con = zeros(T1,N);
flag = 0; timeloop = 0;
while (flag == 0)&&(timeloop <= T1)
    timeloop = timeloop + 1;
    
    %M = eye(N,N) - beta*Lap*dt.*randbin2(N,N,1-beta).*(ones(N,N)-eye(N,N));
    M = eye(N,N) + beta*A*dt.*randbin2(N,N,1-beta).*(ones(N,N)-eye(N,N));
    
    ctmp = M*c;
    
    ctmp2 = ceil(ctmp);
    
    
    c = maskbilevel(ctmp2,0,1,0,1);
    
    
    ctmp3 = R*c;
    ctmp4 = floor(ctmp3);
    c = maskbilevel(ctmp4,0,1,0,1);
    
    
    Con(timeloop,:) = c';
    
    Pop = sum(Con(timeloop,:));
    if Pop == 0
        flag = 1;
    end
    
    
    for nodeloop = 1:N
        node(nodeloop).value = c(nodeloop);
        Rtmp(nodeloop,nodeloop) = c(nodeloop);
    end
    
    if showfig == 1
        drawnet(node,2)
        pause(0.01)
    end
    
    for nloop = 1:N
        if node(nloop).value == 1
            Ct(nloop) = Ct(nloop) + 1/T1;
        end
    end
    
    R = eye(N,N) - (Rtmp.*mu*dt.*eye(N,N).*randbin2(N,N,1-mu));
    
    %     if timeloop > T1-2
    %         keyboard
    %     end
    
end

if Innoc == 1
    % Innocluate
    % Remove the highest-degree node
    displine('avgdegree = ',avgdegree)
    disp(strcat('fac =',num2str(mu/avgdegree/beta)))
    node = subnode(I,node);
    %snode = removenode(I,node);
    [A,degree,Lap] = adjacency(node);
    
    %keyboard
    
    [N,e,avgdegree,maxdegree,mindegree,numclus,meanclus,Lmax,L2,LmaxL2] = clusterstats(node);
    displine('avgdegree = ',avgdegree)
    disp(strcat('fac =',num2str(mu/avgdegree/beta)))
    
    for timeloop = T1+1:T2
        
        %M = eye(N,N) - beta*Lap*dt.*randbin2(N,N,1-beta).*(ones(N,N)-eye(N,N));
        M = eye(N,N) + beta*A*dt.*randbin2(N,N,1-beta).*(ones(N,N)-eye(N,N));
        
        ctmp = M*c;
        
        ctmp2 = ceil(ctmp);
        c = maskbilevel(ctmp2,0,1.01,0,1);
        
        ctmp3 = R*c;
        ctmp4 = floor(ctmp3);
        c = maskbilevel(ctmp4,-0.01,1.01,0,1);
        
        
        Con(timeloop,:) = c';
        
        
        for nodeloop = 1:N
            node(nodeloop).value = c(nodeloop);
            Rtmp(nodeloop,nodeloop) = c(nodeloop);
        end
        
        if showfig == 1
            drawnet(node,2)
            pause(0.01)
        end
        
        R = eye(N,N) - (Rtmp.*mu*dt.*eye(N,N).*randbin2(N,N,1-mu));
        
        %     if timeloop > T1-2
        %         keyboard
        %     end
        
    end
    toc
    
    
    
    x = 0:T2-1;
    h = colormap(jet);
    figure(4)
    for tloop = 1:T2
        Ssum = 0;
        for nodeloop = 1:N
            Ssum =  Ssum + Con(tloop,nodeloop);
        end
        In(tloop) = Ssum;
    end
    
    figure(4)
    plot(In)
    title('Infected')
    
    
    mn1 = mean(In(10:T1));
    mn2 = mean(In(T1+10:T2));
    
    displine('del pop = ',mn2-mn1)
    displine('rel del pop = ',(mn2-mn1)/mn1)
    
    figure(5)
    plot(deg,Ct,'o')
    xlabel('degree')
    ylabel('avg infection')
    
    
else  % non-Innoc case
    
    
    x = 0:T1-1;
    h = colormap(jet);
    figure(4)
    for tloop = 1:T1
        Ssum = 0;
        for nodeloop = 1:N
            Ssum =  Ssum + Con(tloop,nodeloop);
        end
        In(tloop) = Ssum;
    end
    
    figure(4)
    plot(In)
    title('Infected')
    
end