From dfdc6c9966656aa6cb684e8f528fcaa42314543b Mon Sep 17 00:00:00 2001
From: David Nolte <nolte@purdue.edu>
Date: Thu, 4 Feb 2021 06:32:47 -0500
Subject: [PATCH] Create SIRS.m

---
 SIRS.m | 230 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 230 insertions(+)
 create mode 100644 SIRS.m

diff --git a/SIRS.m b/SIRS.m
new file mode 100644
index 0000000..e69930b
--- /dev/null
+++ b/SIRS.m
@@ -0,0 +1,230 @@
+% 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
+
+
+