Friday, June 25, 2010

Solve SIR model -C++/MatLab/Python

1) C++
------------------------------------------------------
/* Solve simple SIR Model
Author: Jiansen Lu
Date: June 25, 2010
Description:
dS/dt = - beta*S*I
dI/dt = beta*S - gamma*I
dR/dt = gamma*I
where S: Susceptibel
I: Infectious
R: Recovered
beta: transmission rate
gamma: recover rate:
set initial condition S(0)=1e-6, I(0)=1.0-S(0)
beta = 1.5, gamma = 0.2
Algorithm: using Runge-Kutta method
dy/dt = f(t,y), y(t_0) = y_0
t_(n+1) = t_n+h
y_(n+1) = y_n + (k_1+2*k_2+2*_k3+k_4)/6
where: k_1=f(t_n,y_n), k_2=f(t_n+h/2,y_n+h*k_1/2)
k_3=f(t_n+h/2,y_n+h*k_2/2), k_4=f(t_n+h,y_n+h*k_3)
*/

#include<iostream>

#include<cmath>

using namespace std;

class SIR{
private:

double t,S,I,R,Pop[3];
double dPop[3],step;

double beta, gamma;
double tmax;
public:


SIR();
SIR(double beta0, double gamma0, double step0, double S00, \
          double I00, double tmax0);
~SIR();

void Diff(double Pop[3]);
void Runge_Kutta();

void Solve_Eq();
};
// Initialise the equations and Runge-Kutta integration
SIR::SIR(double beta0, double gamma0, double step0,double S00, \
       double I00, double tmax0)
{

beta = beta0;
gamma =gamma0;
step = step0;

S =S00;
I = I00;
R = 1 - S - I;

tmax = tmax0;
}
SIR::~SIR(){
cout <<"delete SIR"<<endl;
}

void SIR::Diff(double Pop[3])
{

// The differential equations


dPop[0] = - beta*Pop[0]*Pop[1];              // dS/dt

  dPop[1] = beta*Pop[0]*Pop[1] - gamma*Pop[1];   // dI/dt

  dPop[2] = gamma*Pop[1];                    // dR/dt

}

void SIR::Runge_Kutta(){
int i;
double dPop1[3], dPop2[3], dPop3[3], dPop4[3];

double tmpPop[3], initialPop[3];

/* Integrates the equations one step, using Runge-Kutta 4
Note: we work with arrays rather than variables to make the
coding easier */

initialPop[0]=S; initialPop[1]=I; initialPop[2]=R;

Diff(initialPop);
for(i=0;i<3;i++)
{

dPop1[i]=dPop[i];
tmpPop[i]=initialPop[i]+step*dPop1[i]/2;
}

Diff(tmpPop);
for(i=0;i<3;i++)
{


dPop2[i]=dPop[i];
tmpPop[i]=initialPop[i]+step*dPop2[i]/2;
}

Diff(tmpPop);
for(i=0;i<3;i++)
{

dPop3[i]=dPop[i];
tmpPop[i]=initialPop[i]+step*dPop3[i];
}

Diff(tmpPop);

for(i=0;i<3;i++)
{

dPop4[i]=dPop[i];
tmpPop[i]=initialPop[i]+(dPop1[i]/6 + dPop2[i]/3 + dPop3[i]/3 + dPop4[i]/6)*step;
}


S=tmpPop[0]; I=tmpPop[1]; R=tmpPop[2];


}

void SIR::Solve_Eq(){
t=0;
cout <<"t    S    I       R"<<endl;

do
{
Runge_Kutta();
t+=step;
cout<<t<<"   "<<S<<"   "<<I<<"   "<<R<<"   "<<endl;
}

while(t<tmax);


}
int main(int argc, char** argv)
{

double beta0 = 1.5;
double gamma0 =0.2;

double I00 = 1e-6;
double S00 =1-I00;

double tmax0 = 50;
/* Find a suitable time-scale for outputs */

double step0=0.01/((beta0+gamma0)*S00);

SIR mySIR(beta0, gamma0,step0,S00,  I00,  tmax0);

mySIR.Solve_Eq();
return(0);

}

2) Matlab
in SIR.m
--------------------------------------------------------
function [t,S,I,R] =SIR()
%Solve SIR equation in Matlab

beta=1.5; gamma=0.2; I0=1e-6;
S0=1-I0; tmax=50;
S=S0; I=I0; R=1-S-I;

% The main iteration
[t, pop]=ode45(@Diff_SIR,[0 tmax],[S I R],[],[beta gamma]);

S=pop(:,1); I=pop(:,2); R=pop(:,3);

% plots the graphs with scaled colours
subplot(2,1,1)
h=plot(t,S,'-g',t,R,'-k');
xlabel 'Time';
ylabel 'Susceptibles and Recovereds'

subplot(2,1,2)
h=plot(t,I,'-r');
xlabel 'Time';
ylabel 'Infectious'

% Calculates the differential rates used in the integration.
function dPop=Diff_SIR(t,pop, parameter)
beta=parameter(1); gamma=parameter(2);
S=pop(1); I=pop(2); R=pop(3);
dPop=zeros(3,1);
dPop(1)= -beta*S*I;
dPop(2)= beta*S*I - gamma*I;
dPop(3)= gamma*I;
-------------------------------------------------
3. Python
in SIR.py (chmod  +x SIR.py, ./SIR.py)
------------------------------------------------
#!/usr/bin/env python



import scipy.integrate as spi
import numpy as np
import pylab as pl

beta=1.5
gamma=0.2
TS=1.0
tmax =50
I0=1e-6
S0=1-I0
INPUT = (S0, I0, 0.0)

def diff_eqs(INP,t):
'''The main set of equations'''
Y=np.zeros((3))
V = INP
Y[0] = - beta * V[0] * V[1]
Y[1] = beta * V[0] * V[1] - gamma * V[1]
Y[2] = gamma * V[1]
return Y   # For odeint

t_start = 0.0; t_end = tmax; t_inc = TS
t_range = np.arange(t_start, t_end+t_inc, t_inc)
RES = spi.odeint(diff_eqs,INPUT,t_range)
print RES #Ploting pl.subplot(211) pl.plot(RES[:,0], '-g', label='Susceptibles') pl.plot(RES[:,2], '-k', label='Recovereds') pl.legend(loc=0) pl.title('SIR.py') pl.xlabel('Time') pl.ylabel('Susceptibles and Recovereds') pl.subplot(212) pl.plot(RES[:,1], '-r', label='Infectious') pl.xlabel('Time') pl.ylabel('Infectious') pl.show()

3 comments:

  1. good day. Might I explain in detail the code that you generated in matlab for the SIR model or function that meets this term "pop" and "DPOP"

    ReplyDelete
  2. Hello There,


    Grazie! Grazie! Grazie! Your blog is indeed quite interesting around Jiansen Lu's Computing Blog with you on lot of points!

    A subset of a natural language contains these sentences:
    1) Triangle ABC.
    2) Line Segment AB.
    3) Angle A.
    In this set, names of triangles are formed by putting three letters together,
    all drawn from the alphabet {A,B,C,D,E}. Line segment names are defined by
    putting two letters together from the previous alphabet. And angle names are
    (suprise!) given by any letter in that alphabet.


    Great effort, I wish I saw it earlier. Would have saved my day :)


    Thank you,
    Uday

    ReplyDelete
  3. Python Training, SVR Technologies offers Online & Corporate Classes with Free Live Demo, We cover practical sessions and all the modules in Python.
    python training course

    ReplyDelete