SIRSMC.m

% SIRSMC.m
% converted from diffusionmat on 1/4/16
% Susecpt-Infect-Suscept  plus an innoculation of highest-degree node
% Monte Carlo to get good statistics

clear
close all
format compact

T1 = 100;
T2 = 300;

showfig = 0;

Innoc = 1;  % Turn innoculation on

Nensemble = 100;
for eloop = 1:Nensemble
    eloop
    
    % beta is infection rate
    % mu is recovery rate
    
    %N = 50;  m = 2;    beta = 0.2;  mu = 0.6;               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;   % Infect the high-degree node
    
    if showfig == 1
        fh2 = figure(2);
        drawnet(node)
    end
    
    % 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(T2,N);
    in = zeros(1,T2);
    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)))
    end
        
        
        while (flag == 0)&&(timeloop <= T2)
            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.01,0,1);
            
            ctmp3 = R*c;
            ctmp4 = floor(ctmp3);
            c = maskbilevel(ctmp4,-0.01,1.01,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
            
            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
    
    Infection(eloop,:) = In;
    
end     % end eloop

figure(6)
imagesc(Infection)

Y = mean(Infection);

figure(7)
plot(Y)