[Scilab-users] Covid19 model

classic Classic list List threaded Threaded
3 messages Options
Claus Futtrup Claus Futtrup
Reply | Threaded
Open this post in threaded view
|

[Scilab-users] Covid19 model

Hi Scilabers

On 17. April I sent an email, but I think it wasn't released to the mailing list, I only find it online here:

http://mailinglists.scilab.org/Scilab-users-Covid19-model-td4040626.html

Here's the email again (without attachment, see above link if you'd like to view the attachment)

A friend (Lars Risbo) published some MATLAB code in LinkedIn for simulating the infection with a company lockdown after some time. I figured I'd try to convert it to Scilab.

The code uses ODE45 and I'm not sure that I understand how to convert this MATLAB code to Scilab. I hope you can explain what I need to do to make the code work. Some of the original MATLAB code is found in the comments.

// covid19risbo.sce
//
// SEIRsim1
// Susceptible-Exposed-Infectious-Recovered (SEIR)

function dydt=odefun(t, y)
    if t<days
      R=R0;
    else
      R=0.6; // Change reproductive rate after lockdown
    end
    A=[0    0   -delta*R*y(1)/Npop  0 0 0;...
      0   -gam delta*R*y(1)/Npop  0 0 0;...
      0    gam -delta        0 0 0;...
      0    0   delta*(1-Fhosp)   0 0 0;...
      0    0   delta*Fhosp     0 -1/Thosp 0;...
      0    0   delta*R*y(1)/Npop  0 0 0];
    dydt=A*y;
end

Begin=datenum(2020,02,20,0,0,0); // begin date
Dlock=datenum(2020,03,12,0,0,0); // date of lockdown
days=Dlock-Begin;
R0=2.6;   // inital R value
gam= 1/3; // gamma, the inverse of average latent time
delta= 1/5; // inv time constant which infectious people either recover or enter hospital
Fhosp=0.16; // fraction of recovering people going to hospital
Thosp=14;  // average time of hospitalisation
Npop=6e6;  // total initial population of sensitive
y0=[Npop;50;50;0;0;0];// [S E I R Hosp ] intial cond.
t1 = 0:60;
y1 = ode(y0, 0, 60, odefun); // [t1,y1] = ode45(@odefun,[0 60],y0); // run 1st scernario
idx=[3:6];
scf();
a = gca();
plot(Begin+t1,y1(:,idx),'-'); // semilogy(Begin+t1,y1(:,idx),'-','LineWidth',3)
a.log_flags = "nln";
xlabel('Date')
ylabel('Number of cases')
xgrid();
// ax=gca;
// ax.YLim=[10 max(max(y2(:,idx)))];
xtitle({'Danish Corona lock down on 20.03.12, vs 14 days later','R goes from 2.6 to 0.6 at lockdown'});
legend('Infected','Recovered','Hospitalised','Total Cases');

P.S. Attached a simple graph of what the output should look like (file covid19risbo.jpg < 50 kb), at least partially - because I deleted an alternative case with a later lockdown date = dashed lines.

Best regards,

Claus


_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users
mottelet mottelet
Reply | Threaded
Open this post in threaded view
|

Re: Covid19 model

Hi Clauss,

Just modify the following lines :

t1 = 0:0.1:60;
y1 = ode(y0, 0, t1, odefun); // [t1,y1] = ode45(@odefun,[0 60],y0); // run 1st scernario
idx=[3:6];
scf();
a = gca();
plot(Begin+t1,y1(idx,:),'-'); // semilogy(Begin+t1,y1(:,idx),'-','LineWidth',3)

And you will be fine. ode() needs a vector of time values and the output is oriented transposed w.r.t Matlab output.

S.

Le 30/04/2020 à 16:40, Claus Futtrup a écrit :
function dydt=odefun(t, y)
    if t<days
      R=R0;
    else
      R=0.6; // Change reproductive rate after lockdown
    end
    A=[0    0   -delta*R*y(1)/Npop  0 0 0;...
      0   -gam delta*R*y(1)/Npop  0 0 0;...
      0    gam -delta        0 0 0;...
      0    0   delta*(1-Fhosp)   0 0 0;...
      0    0   delta*Fhosp     0 -1/Thosp 0;...
      0    0   delta*R*y(1)/Npop  0 0 0];
    dydt=A*y;
end

Begin=datenum(2020,02,20,0,0,0); // begin date
Dlock=datenum(2020,03,12,0,0,0); // date of lockdown
days=Dlock-Begin;
R0=2.6;   // inital R value
gam= 1/3; // gamma, the inverse of average latent time
delta= 1/5; // inv time constant which infectious people either recover or enter hospital
Fhosp=0.16; // fraction of recovering people going to hospital
Thosp=14;  // average time of hospitalisation
Npop=6e6;  // total initial population of sensitive
y0=[Npop;50;50;0;0;0];// [S E I R Hosp ] intial cond.
t1 = 0:60;
y1 = ode(y0, 0, 60, odefun); // [t1,y1] = ode45(@odefun,[0 60],y0); // run 1st scernario
idx=[3:6];
scf();
a = gca();
plot(Begin+t1,y1(:,idx),'-'); // semilogy(Begin+t1,y1(:,idx),'-','LineWidth',3)
a.log_flags = "nln";
xlabel('Date')
ylabel('Number of cases')
xgrid();
// ax=gca;
// ax.YLim=[10 max(max(y2(:,idx)))];
xtitle({'Danish Corona lock down on 20.03.12, vs 14 days later','R goes from 2.6 to 0.6 at lockdown'});
legend('Infected','Recovered','Hospitalised','Total Cases');
-- 
Stéphane Mottelet
Ingénieur de recherche
EA 4297 Transformations Intégrées de la Matière Renouvelable
Département Génie des Procédés Industriels
Sorbonne Universités - Université de Technologie de Compiègne
CS 60319, 60203 Compiègne cedex
Tel : +33(0)344234688
http://www.utc.fr/~mottelet

_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users
Claus Futtrup Claus Futtrup
Reply | Threaded
Open this post in threaded view
|

Re: Covid19 model

Hi Stéphane

Thank you for showing me how to work ODE.

P.S. about Corona modeling, the model by Risbo is nice in that you can implement a lockdown of some sort. It's a simple one-step lockdown, where you have to tune the R-value ... in reality each geographic area is different, related to how much of a lockdown is executed, how dense the population is, etc.

Cheers,
Claus

On 4/30/2020 4:51 PM, Stéphane Mottelet wrote:

Hi Clauss,

Just modify the following lines :

t1 = 0:0.1:60;
y1 = ode(y0, 0, t1, odefun); // [t1,y1] = ode45(@odefun,[0 60],y0); // run 1st scernario
idx=[3:6];
scf();
a = gca();
plot(Begin+t1,y1(idx,:),'-'); // semilogy(Begin+t1,y1(:,idx),'-','LineWidth',3)

And you will be fine. ode() needs a vector of time values and the output is oriented transposed w.r.t Matlab output.

S.

Le 30/04/2020 à 16:40, Claus Futtrup a écrit :
function dydt=odefun(t, y)
    if t<days
      R=R0;
    else
      R=0.6; // Change reproductive rate after lockdown
    end
    A=[0    0   -delta*R*y(1)/Npop  0 0 0;...
      0   -gam delta*R*y(1)/Npop  0 0 0;...
      0    gam -delta        0 0 0;...
      0    0   delta*(1-Fhosp)   0 0 0;...
      0    0   delta*Fhosp     0 -1/Thosp 0;...
      0    0   delta*R*y(1)/Npop  0 0 0];
    dydt=A*y;
end

Begin=datenum(2020,02,20,0,0,0); // begin date
Dlock=datenum(2020,03,12,0,0,0); // date of lockdown
days=Dlock-Begin;
R0=2.6;   // inital R value
gam= 1/3; // gamma, the inverse of average latent time
delta= 1/5; // inv time constant which infectious people either recover or enter hospital
Fhosp=0.16; // fraction of recovering people going to hospital
Thosp=14;  // average time of hospitalisation
Npop=6e6;  // total initial population of sensitive
y0=[Npop;50;50;0;0;0];// [S E I R Hosp ] intial cond.
t1 = 0:60;
y1 = ode(y0, 0, 60, odefun); // [t1,y1] = ode45(@odefun,[0 60],y0); // run 1st scernario
idx=[3:6];
scf();
a = gca();
plot(Begin+t1,y1(:,idx),'-'); // semilogy(Begin+t1,y1(:,idx),'-','LineWidth',3)
a.log_flags = "nln";
xlabel('Date')
ylabel('Number of cases')
xgrid();
// ax=gca;
// ax.YLim=[10 max(max(y2(:,idx)))];
xtitle({'Danish Corona lock down on 20.03.12, vs 14 days later','R goes from 2.6 to 0.6 at lockdown'});
legend('Infected','Recovered','Hospitalised','Total Cases');
-- 
Stéphane Mottelet
Ingénieur de recherche
EA 4297 Transformations Intégrées de la Matière Renouvelable
Département Génie des Procédés Industriels
Sorbonne Universités - Université de Technologie de Compiègne
CS 60319, 60203 Compiègne cedex
Tel : +33(0)344234688
http://www.utc.fr/~mottelet

_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users



_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users