update covid19 file

CRA_Updated_Version/CRA_Toolbox_bash/Examples/ITALIA/init_italia.m

+
+region_name='Italy_pred';
+
+%define initial conditions
+N = 60e6;         %population Italy
+
+%define time
+%time = 0:1:110;
+
+E00 = (1/N)*10^5;
+S0 = 10^5-E00;
+x0 = [S0 E00 0 0 0 0 0 0 0];
+        
+%inizio epidemia: 27 Gennaio
+%dopo 24 giorni (20 Febbraio): paziente 1
+%dopo 25 giorni (21 Febbraio): zona rossa in Lombardia e Veneto
+%dopo 28 giorni (24 Febbraio): chiusura scuole in Emilia-Romagna, Friuli-Venezia Giulia,
+%Liguria, Lombardia, parte delle Marche, Piemonte, Veneto
+%dopo 25 giorni (21 Febbraio) DPI
+%dopo 38 giorni (5 Marzo): chiusura scuole nel resto d'Italia
+%dopo 40 giorni (7 Marzo): fuga verso sud
+%dopo 41 giorni (8 Marzo): lockdown nord Italia + social distancing +
+%eventi vietati 
+%dopo 43 giorni (10 Marzo): total lockdown
+
+%RILASCIO DEL LOCKDOWN:
+%dopo 98 giorni (4 Maggio 2020)
+%dopo 112 giorni (18 Maggio 2020)
+
+Tlock=[25 28 38 43 98 112];
+        
+%s0 = [0.8 0.75 0.7 0.2];  %reduction of 30%, 70%, 80%
+%s1=[0.7 0.65 0.6 0.15];  %reduction of 40%, 75%, 85% 
+%s2=0.4;  %reduction of 60%
\ No newline at end of file

CRA_Updated_Version/CRA_Toolbox_bash/Examples/compute_intervention.m

+
+
+function [s0,s1,s2]=compute_intervention(t,Tlock,s00,s11,s22,region_name)
+
+
+    if(strcmp(region_name,'Italy'))
+       s0 = (t<Tlock(1)) + s00(1) * (t>=Tlock(1)) * (t<Tlock(2)) + s00(2)*(t>=Tlock(2)) * (t<Tlock(3)) + s00(3)*(t>=Tlock(3))*(t<Tlock(4)) + s00(4)*(t>=Tlock(4));%*(t<Tlock(5)) + s00(5)*(t>=Tlock(5))*(t<Tlock(6)) + s00(6)*(t>=Tlock(6));
+       s1 = (t<Tlock(3)) + s11(1) * (t>=Tlock(3));
+       s2 = (t<Tlock(1)) + s22 * (t>=Tlock(1));        
+    end
+
+    if(strcmp(region_name,'Italy_pred'))
+       s0 = (t<Tlock(1)) + s00(1) * (t>=Tlock(1)) * (t<Tlock(2)) + s00(2)*(t>=Tlock(2)) * (t<Tlock(3)) + s00(3)*(t>=Tlock(3))*(t<Tlock(4)) + s00(4)*(t>=Tlock(4))*(t<Tlock(5)) + s00(5)*(t>=Tlock(5))*(t<Tlock(6)) + s00(6)*(t>=Tlock(6));
+       s1 = (t<Tlock(3)) + s11(1) * (t>=Tlock(3));
+       s2 = (t<Tlock(1)) + s22 * (t>=Tlock(1));        
+    end
+    
+    if(strcmp(region_name,'Umbria'))
+       s0 = (t<Tlock(1)) + s00(1) * (t>=Tlock(1)) * (t<Tlock(2)) + s00(2)*(t>=Tlock(2)) * (t<Tlock(3)) + s00(3)*(t>=Tlock(3));%*(t<Tlock(5)) + s00(5)*(t>=Tlock(5))*(t<Tlock(6)) + s00(6)*(t>=Tlock(6));
+       s1 = (t<Tlock(3)) + s11(1) * (t>=Tlock(3));
+       s2 = (t<Tlock(1)) + s22 * (t>=Tlock(1));        
+    end
+    
+    if(strcmp(region_name,'Campania'))
+       s0 = (t<Tlock(2)) + s00(1) * (t>=Tlock(2)) * (t<Tlock(3)) + s00(2)*(t>=Tlock(3)); %*(t<Tlock(5)) + s00(5)*(t>=Tlock(5))*(t<Tlock(6)) + s00(6)*(t>=Tlock(6));
+       s1 = (t<Tlock(3)) + s11(1) * (t>=Tlock(3));
+       s2 = (t<Tlock(1)) + s22 * (t>=Tlock(1));        
+    end
+    
+        
+    if(strcmp(region_name,'Lazio'))
+       s0 = (t<Tlock(2)) + s00(1) * (t>=Tlock(2)) * (t<Tlock(3)) + s00(2)*(t>=Tlock(3));%*(t<Tlock(5)) + s00(5)*(t>=Tlock(5))*(t<Tlock(6)) + s00(6)*(t>=Tlock(6));
+       s1 = (t<Tlock(3)) + s11(1) * (t>=Tlock(3));
+       s2 = (t<Tlock(1)) + s22 * (t>=Tlock(1));        
+    end
+    
+    
+    if(strcmp(region_name,'Lombardia'))
+       s0 = (t<Tlock(1)) + s00(1) * (t>=Tlock(1)) * (t<Tlock(2)) + s00(2)*(t>=Tlock(2)) * (t<Tlock(3)) + s00(3)*(t>=Tlock(3));%*(t<Tlock(5)) + s00(5)*(t>=Tlock(5))*(t<Tlock(6)) + s00(6)*(t>=Tlock(6));
+       s1 = (t<Tlock(1)) + s11 * (t>=Tlock(3));
+       s2 = (t<Tlock(1)) + s22 * (t>=Tlock(1));        
+    end
+    
+   if(strcmp(region_name,'Veneto'))
+       s0 = (t<Tlock(1)) + s00(1) * (t>=Tlock(1)) * (t<Tlock(2)) + s00(2)*(t>=Tlock(2)) * (t<Tlock(3)) + s00(3)*(t>=Tlock(3));%*(t<Tlock(5)) + s00(5)*(t>=Tlock(5))*(t<Tlock(6)) + s00(6)*(t>=Tlock(6));
+       s1 = (t<Tlock(3)) + s11 * (t>=Tlock(3));
+       s2 = (t<Tlock(1)) + s22 * (t>=Tlock(1));        
+   end
+    
+   
+   if(strcmp(region_name,'EmiliaRomagna'))
+       s0 = (t<Tlock(2)) + s00(1) * (t>=Tlock(2)) * (t<Tlock(3)) + s00(2)*(t>=Tlock(3));%*(t<Tlock(5)) + s00(5)*(t>=Tlock(5))*(t<Tlock(6)) + s00(6)*(t>=Tlock(6));
+       s1 = (t<Tlock(3)) + s11 * (t>=Tlock(3));
+       s2 = (t<Tlock(1)) + s22 * (t>=Tlock(1));        
+   end
+    
+   
+   if(strcmp(region_name,'Calabria'))
+       s0 = (t<Tlock(1)) + s00(1) * (t>=Tlock(1)) * (t<Tlock(2)) + s00(2)*(t>=Tlock(2));%*(t<Tlock(5)) + s00(5)*(t>=Tlock(5))*(t<Tlock(6)) + s00(6)*(t>=Tlock(6));
+       s1 = (t<Tlock(2)) + s11(1) * (t>=Tlock(2));
+       s2 = (t<Tlock(1)) + s22 * (t>=Tlock(1));        
+   end
+   
+      
+   if(strcmp(region_name,'Abruzzo'))
+       s0 = (t<Tlock(2)) + s00(1) * (t>=Tlock(2)) * (t<Tlock(3)) + s00(2)*(t>=Tlock(3));%*(t<Tlock(5)) + s00(5)*(t>=Tlock(5))*(t<Tlock(6)) + s00(6)*(t>=Tlock(6));
+       s1 = (t<Tlock(3)) + s11(1) * (t>=Tlock(3));
+       s2 = (t<Tlock(1)) + s22 * (t>=Tlock(1));        
+   end
+   
+   if(strcmp(region_name,'Basilicata'))
+       s0 = (t<Tlock(1)) + s00(1) * (t>=Tlock(1)) * (t<Tlock(2)) + s00(2)*(t>=Tlock(2));%*(t<Tlock(5)) + s00(5)*(t>=Tlock(5))*(t<Tlock(6)) + s00(6)*(t>=Tlock(6));
+       s1 = (t<Tlock(2)) + s11(1) * (t>=Tlock(2));
+       s2 = (t<Tlock(1)) + s22 * (t>=Tlock(1));        
+   end
+   
+   if(strcmp(region_name,'Marche'))
+       s0 = (t<Tlock(1)) + s00(1) * (t>=Tlock(1)) * (t<Tlock(2)) + s00(2)*(t>=Tlock(2));%*(t<Tlock(5)) + s00(5)*(t>=Tlock(5))*(t<Tlock(6)) + s00(6)*(t>=Tlock(6));
+       s1 = (t<Tlock(2)) + s11(1) * (t>=Tlock(2));
+       s2 = (t<Tlock(1)) + s22 * (t>=Tlock(1));        
+   end
+   
+   if(strcmp(region_name,'Piemonte'))
+       s0 = (t<Tlock(1)) + s00(1) * (t>=Tlock(1)) * (t<Tlock(2)) + s00(2)*(t>=Tlock(2));%*(t<Tlock(5)) + s00(5)*(t>=Tlock(5))*(t<Tlock(6)) + s00(6)*(t>=Tlock(6));
+       s1 = (t<Tlock(2)) + s11(1) * (t>=Tlock(2));
+       s2 = (t<Tlock(1)) + s22 * (t>=Tlock(1));        
+   end
+       
+   if(strcmp(region_name,'Toscana'))
+       s0 = (t<Tlock(2)) + s00(1) * (t>=Tlock(2)) * (t<Tlock(3)) + s00(2)*(t>=Tlock(3));%*(t<Tlock(5)) + s00(5)*(t>=Tlock(5))*(t<Tlock(6)) + s00(6)*(t>=Tlock(6));
+       s1 = (t<Tlock(3)) + s11(1) * (t>=Tlock(3));
+       s2 = (t<Tlock(1)) + s22 * (t>=Tlock(1));        
+   end
+   
+   if(strcmp(region_name,'Puglia'))
+       s0 = (t<Tlock(2)) + s00(1) * (t>=Tlock(2)) * (t<Tlock(3)) + s00(2)*(t>=Tlock(3));%*(t<Tlock(5)) + s00(5)*(t>=Tlock(5))*(t<Tlock(6)) + s00(6)*(t>=Tlock(6));
+       s1 = (t<Tlock(3)) + s11(1) * (t>=Tlock(3));
+       s2 = (t<Tlock(1)) + s22 * (t>=Tlock(1));        
+   end
+    
+   if(strcmp(region_name,'Liguria'))
+       s0 = (t<Tlock(1)) + s00(1) * (t>=Tlock(1)) * (t<Tlock(2)) + s00(2)*(t>=Tlock(2));%*(t<Tlock(5)) + s00(5)*(t>=Tlock(5))*(t<Tlock(6)) + s00(6)*(t>=Tlock(6));
+       s1 = (t<Tlock(2)) + s11(1) * (t>=Tlock(2));
+       s2 = (t<Tlock(1)) + s22 * (t>=Tlock(1));        
+   end
+   
+   if(strcmp(region_name,'Sicilia'))
+       s0 = (t<Tlock(2)) + s00(1) * (t>=Tlock(2)) * (t<Tlock(3)) + s00(2)*(t>=Tlock(3));%*(t<Tlock(5)) + s00(5)*(t>=Tlock(5))*(t<Tlock(6)) + s00(6)*(t>=Tlock(6));
+       s1 = (t<Tlock(3)) + s11(1) * (t>=Tlock(3));
+       s2 = (t<Tlock(1)) + s22 * (t>=Tlock(1));        
+   end
+   
+   if(strcmp(region_name,'Sardegna'))
+       s0 = (t<Tlock(2)) + s00(1) * (t>=Tlock(2)) * (t<Tlock(3)) + s00(2)*(t>=Tlock(3));%*(t<Tlock(5)) + s00(5)*(t>=Tlock(5))*(t<Tlock(6)) + s00(6)*(t>=Tlock(6));
+       s1 = (t<Tlock(3)) + s11(1) * (t>=Tlock(3));
+       s2 = (t<Tlock(1)) + s22 * (t>=Tlock(1));        
+   end
+   
+   if(strcmp(region_name,'Friuli'))
+       s0 = (t<Tlock(1)) + s00(1) * (t>=Tlock(1)) * (t<Tlock(3)) + s00(2)*(t>=Tlock(3));%*(t<Tlock(5)) + s00(5)*(t>=Tlock(5))*(t<Tlock(6)) + s00(6)*(t>=Tlock(6));
+       s1 = (t<Tlock(3)) + s11(1) * (t>=Tlock(3));
+       s2 = (t<Tlock(2)) + s22 * (t>=Tlock(2));        
+   end
+   
+   if(strcmp(region_name,'Aosta'))
+       s0 = (t<Tlock(1)) + s00(1) * (t>=Tlock(1)) * (t<Tlock(3)) + s00(2)*(t>=Tlock(3));%*(t<Tlock(5)) + s00(5)*(t>=Tlock(5))*(t<Tlock(6)) + s00(6)*(t>=Tlock(6));
+       s1 = (t<Tlock(3)) + s11(1) * (t>=Tlock(3));
+       s2 = (t<Tlock(2)) + s22 * (t>=Tlock(2));        
+   end
+   
+   if(strcmp(region_name,'Molise'))
+       s0 = (t<Tlock(2)) + s00(1) * (t>=Tlock(2)) * (t<Tlock(3)) + s00(2)*(t>=Tlock(3));%*(t<Tlock(5)) + s00(5)*(t>=Tlock(5))*(t<Tlock(6)) + s00(6)*(t>=Tlock(6));
+       s1 = (t<Tlock(3)) + s11(1) * (t>=Tlock(3));
+       s2 = (t<Tlock(1)) + s22 * (t>=Tlock(1));        
+   end
+   
+      
+   if(strcmp(region_name,'Trentino'))
+       s0 = (t<Tlock(2)) + s00(1) * (t>=Tlock(2)) * (t<Tlock(3)) + s00(2)*(t>=Tlock(3));%*(t<Tlock(5)) + s00(5)*(t>=Tlock(5))*(t<Tlock(6)) + s00(6)*(t>=Tlock(6));
+       s1 = (t<Tlock(3)) + s11(1) * (t>=Tlock(3));
+       s2 = (t<Tlock(1)) + s22 * (t>=Tlock(1));        
+   end
+end

CRA_Updated_Version/CRA_Toolbox_bash/Examples/ode_covid19_v3.m → CRA_Updated_Version/CRA_Toolbox_bash/Examples/ode_covid19_v3_NO.m

@@ -53,7 +53,7 @@ u=(1/TimeICUDeath)*(CFR/FracCritical);
 g3=(1/TimeICUDeath)-u;
 
 
-s0 = (t<Tlock(2)) + s00(1) * (t>=Tlock(2)) * (t<Tlock(3)) + s00(2)*(t>=Tlock(3)) * (t<Tlock(4)) + s00(3)*(t>=Tlock(4));
+s0 = (t<Tlock(2)) + s00(1) * (t>=Tlock(2)) * (t<Tlock(3)) + s00(2)*(t>=Tlock(3)) * (t<Tlock(4)) + s00(3)*(t>=Tlock(4))* (t<Tlock(5)) + s00(4)*(t>=Tlock(5))* (t<Tlock(6)) + s00(5)* (t>=Tlock(6));
 s1 = (t<Tlock(2)) + s11(1) * (t>=Tlock(2)) * (t<Tlock(3)) + s11(2)*(t>=Tlock(3))* (t<Tlock(4)) + s11(3)*(t>=Tlock(4));
 s2 = (t<Tlock(1)) + s22 * (t>=Tlock(1));
 %s3 = (t<Tlock(1)) + s33 * (t>=Tlock(1));

CRA_Updated_Version/CRA_Toolbox_bash/Examples/ode_covid19_v3_OK.m

+%ODE model to describe the spread and clinical progression of COVID-19
+
+function dx=ode_covid19_v3_OK(t,x,p,Tlock,region_name) %s00,s11,s22,region_name)%,Tlock2,s02,s12)
+
+be = p(1);
+b0=p(2);
+b1=p(3);
+b2=p(4);
+b3=p(5);
+        
+Norm=10^5;
+
+be=be/(Norm);
+b0=b0/(Norm);
+b1=b1/(Norm);
+b2=b2/(Norm);
+b3=b3/(Norm);
+        
+FracSevere = p(6);
+FracCritical = p(7);
+FracMild = 1-FracSevere-FracCritical;
+FracAsym = p(8);
+               
+IncubPeriod = p(9);
+DurMildInf = p(10);
+
+DurAsym = p(11);
+   
+DurHosp = p(12);
+TimeICUDeath = p(13);
+ProbDeath=p(14);
+
+PresymPeriod=p(15)*IncubPeriod;
+           
+CFR=ProbDeath*FracCritical/100; 
+            
+a1=1/PresymPeriod; %presymptomatic period of transmission
+a0=1/(IncubPeriod-PresymPeriod); %true latent period  
+
+f=FracAsym;
+    
+g0=1/DurAsym;
+
+g1=(1/DurMildInf)*FracMild;
+p1=(1/DurMildInf)-g1;
+    
+p2=(1/DurHosp)*(FracCritical/(FracSevere+FracCritical));
+g2=(1/DurHosp)-p2;
+    
+u=(1/TimeICUDeath)*(CFR/FracCritical);
+    
+g3=(1/TimeICUDeath)-u;
+
+s00=[p(16) p(17) p(18) p(19) p(20) p(21)];
+s11=p(22);
+s22=p(23);
+
+[s0,s1,s2]=compute_intervention(t,Tlock,s00,s11,s22,region_name);
+
+
+%susceptibles individuals (not infected)
+S = x(1);
+%exposed individuals, infected but no symptoms and no transmission
+E0 = x(2);
+%exposed individuals, infected but no symptoms and can transmit
+E1 = x(3);
+%infected individuals but asymptomatic infection
+I0 = x(4);
+%infected individuals, mild infection
+I1 = x(5);
+%infected individuals, severe infection
+I2 = x(6);
+%infected individuals, critical infection
+I3 = x(7);
+%recovered
+R = x(8);
+%dead
+D = x(9);
+
+
+%dS = - (be *(1/(1+s0^n))* E1 + b0*(1/(1+s0^n))*I0 + b1*(1/(1+s1^n))*I1 + b2*(1/(1+s2^n))*I2 + b3*(1/(1+s2^n))*I3)*S; % + k_e0s * E0;
+%dE0 = (be *(1/(1+s0^n))* E1 + b0*(1/(1+s0^n))*I0 + b1*(1/(1+s1^n))*I1 + b2*(1/(1+s2^n))*I2 + b3*(1/(1+s2^n))*I3)*S -a0*E0; % - k_e0s * E0; 
+dS = - ((be * s0)* E1 + (b0*s0)*I0 + (b1*s1)*I1 + (b2*s2)*I2 + (b3*s2)*I3)*S; % + k_e0s * E0;
+dE0 = ((be *s0)* E1 + (b0*s0)*I0 + (b1*s1)*I1 + (b2*s2)*I2 + (b3*s2)*I3)*S -a0*E0; % - k_e0s * E0; 
+dE1 = a0 * E0 - a1*E1; % -g4*E1;
+dI0 = f * a1 * E1 - g0*I0;
+dI1 = (1-f) * a1 * E1 - g1*I1 - p1*I1;
+dI2 = p1 * I1 - g2*I2 - p2*I2;  % -u2*I2;
+dI3 = p2 * I2 - g3 *I3 - u * I3 ;
+dR = g0*I0 + g1*I1 + g2*I2 + g3*I3; % +g4 * E1;
+dD = u * I3; % + u2 * I2;
+
+
+
+dx=[dS;dE0;dE1;dI0;dI1;dI2;dI3;dR;dD];

CRA_Updated_Version/CRA_Toolbox_bash/compute_MIRI.m

@@ -88,7 +88,7 @@ try
         %values of the evaluation function are computed and stored in the
         %TimeBehavior object Results
         Results.currentEvalFunc=current_func;
-        Results.computeEvalFunc();  
+        Results.computeEvalFunc();
         
         samples=Results.evalFuncValues;
         
@@ -124,6 +124,7 @@ try
         Results.currentTailMethod=current_tm;
         
         [XiMax,XiMin]=Results.compute_tail(AllPerturbations{k,1},tail_size);
+        [XiMaxD,XiMinD]=Results.compute_tail(AllDerivedParam{k,1},tail_size);
         %[XiMax,XiMin]=handles.Results.compute_tail(handles1.perturbation,handles.tail_size);
         
         NcloudMax=size(XiMax,1);
@@ -139,7 +140,7 @@ try
         
         if(isfield(model,'derived_parameters'))
            obj_MIRI_derived=MIRI();
-           [all_xbin_max_d,all_xbin_min_d,all_ks_max_d,all_ks_min_d]=obj_MIRI_derived.estimate_parampdfs(derived_parameters, AllDerivedParam{k,1}, XiMax, XiMin);
+           [all_xbin_max_d,all_xbin_min_d,all_ks_max_d,all_ks_min_d]=obj_MIRI_derived.estimate_parampdfs(derived_parameters, AllDerivedParam{k,1}, XiMaxD, XiMinD);
            %[all_xbin_max,all_xbin_min,all_ks_max,all_ks_min]=obj_MIRI.estimate_parampdfs(handles1.nominal_p, handles1.perturbation, Ncloud, handles.XiMax, handles.XiMin);
            [MIRI_param_d,CloudCondMin_d, CloudCondMax_d]=obj_MIRI_derived.MIRI_computation(all_ks_max_d,all_ks_min_d,all_xbin_max_d,all_xbin_min_d);