From 774a357a8c53f17acd04abe204930ea3ce5e0dae Mon Sep 17 00:00:00 2001 From: moeinkh88 Date: Wed, 23 Nov 2022 18:30:34 +0200 Subject: [PATCH] clear the plots --- Allplots.svg | 37342 +++++++++++++++++++++++++++++++ Allplots_Models.jl | 298 + Optim-FDE-Data-Order.jl | 239 + PLOT-Code.jl | 141 - PLOT-Code1.jl | 139 - PLOT-Code2.jl | 120 - PLOT-Code3.jl | 145 - TestNewModel.jl | 117 + Turring-ODE--ConstantN-Data.jl | 133 + Turring-ODE-Data.jl | 266 + pltFtry.svg | 3153 +++ pltheat.svg | 385 + pltheat1.svg | 381 + 13 files changed, 42314 insertions(+), 545 deletions(-) create mode 100644 Allplots.svg create mode 100644 Allplots_Models.jl create mode 100644 Optim-FDE-Data-Order.jl delete mode 100644 PLOT-Code.jl delete mode 100644 PLOT-Code1.jl delete mode 100644 PLOT-Code2.jl delete mode 100644 PLOT-Code3.jl create mode 100644 TestNewModel.jl create mode 100644 Turring-ODE--ConstantN-Data.jl create mode 100644 Turring-ODE-Data.jl create mode 100644 pltFtry.svg create mode 100644 pltheat.svg create mode 100644 pltheat1.svg diff --git a/Allplots.svg b/Allplots.svg new file mode 100644 index 0000000..f7dda3b --- /dev/null +++ b/Allplots.svg @@ -0,0 +1,37342 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/Allplots_Models.jl b/Allplots_Models.jl new file mode 100644 index 0000000..933362a --- /dev/null +++ b/Allplots_Models.jl @@ -0,0 +1,298 @@ +# plots ODE + FDE, variable and constant N, heatmaps +using CSV +using DataFrames +using FdeSolver, Plots,SpecialFunctions, Optim, StatsBase, Random +using Interpolations, LinearAlgebra + + +# Dataset +Data=CSV.read("datoswho.csv", DataFrame) +data=(Matrix(Float64.(Data))) +Data2=LinearInterpolation(data[:,1], data[:,2]) #Interpolate the data + +#initial conditons and parameters + +x0=[18000,0,15,0,0,0,0,0]# initial conditons S0,E0,I0,R0,L0,H0,B0,C0 + + α1=35.37e-3 # Density independent part of the birth rate for individuals. + α2=8.334570273681649e-8# Density dependent part of the birth rate for individuals. + σ=1/11.4 # Per capita rate at which exposed individuals become infectious. + γ1=0.1 # Per capita rate of progression of individuals from the infectious class to the asymptomatic class. + γ2=1/5 # Per capita rate of progression of individuals from the hospitalized class to the asymptomatic class. + γ3=1/30 # Per capita recovery rate of individuals from the asymptomatic class to the complete recovered class. + ϵ=1/9.6 # Fatality rate. + δ1=1/2 # Per capita rate of progression of individuals from the dead class to the buried class. + δ2=1/4.6 # Per capita rate of progression of individuals from the hospitalized class to the buried class. + τ=1/5 # Per capita rate of progression of individuals from the infectious class to the hospitalized class. + βi=0.14 # Contact rate of infective individuals and susceptible. + βd=0.4 # Contact rate of infective individuals and dead. + βh=0.29 # Contact rate of infective individuals and hospitalized. + βr=0.2174851498937417 # Contact rate of infective individuals and asymptomatic. + μ=14e-3 + μ1=10.17e-3 # Density independent part of the death rate for individuals. + μ2=4.859965522407347e-7 # Density dependent part of the death rate for individuals. + ξ=14e-3 # Incineration rate + + par=[α1,σ,γ1,γ2,γ3,ϵ,δ1,δ2,τ,βi,βd,βh,βr,μ1,ξ,α2,μ2] + +#Define the equation with variable N +function F(t, x, par) + + α1,σ,γ1,γ2,γ3,ϵ,δ1,δ2,τ,βi,βd,βh,βr,μ1,ξ,α2,μ2=par + S, E, I1, R, L, H, B, C=x + N=sum(x) + + dS=(α1-α2*N)*N - βi/N*S*I1 - βh/N*S*H - βd/N*S*L - βr/N*S*R - (μ1+μ2*N)*S + dE=βi/N*S*I1 + βh/N*S*H + βd/N*S*L + βr/N*S*R - σ*E - (μ1+μ2*N)*E + dI=σ*E - (γ1 + ϵ + τ + μ1 + μ2*N)*I1 + dR=γ1*I1 + γ2*H - (γ3 + μ1 + μ2*N)*R + dL=ϵ*I1 - (δ1+ξ)*L + dH=τ*I1 - (γ2 + δ2 + μ1 + μ2*N)*H + dB=δ1*L + δ2*H - ξ*B + dC=γ3*R - (μ1 + μ2*N)*C + + return [dS, dE, dI, dR, dL, dH, dB, dC] + +end + +#eq with constant N +parC=[α1,σ,γ1,γ2,γ3,ϵ,δ1,δ2,τ,βi,βd,βh,βr,μ1,ξ,sum(x0)] +function Fc(t, x, par) + + α1,σ,γ1,γ2,γ3,ϵ,δ1,δ2,τ,βi,βd,βh,βr,μ1,ξ,N=parC + S, E, I, R, L, H, B, C=x + βr=0.7265002737432911 + μ=0.06918229886616623 + + dS=μ*N - βi/N*S*I - βh/N*S*H - βd/N*S*L - βr/N*S*R - μ*S + dE=βi/N*S*I + βh/N*S*H + βd/N*S*L + βr/N*S*R - σ*E - μ*E + dI=σ*E - (γ1 + ϵ + τ + μ)*I + dR=γ1*I + γ2*H - (γ3 + μ)*R + dL=ϵ*I - (δ1+ξ)*L + dH=τ*I - (γ2 + δ2 + μ)*H + dB=δ1*L + δ2*H - ξ*B + dC=γ3*R - μ*C + + return [dS, dE, dI, dR, dL, dH, dB, dC] + +end + +## +T1=250 + T2=438 + tSpan=[4,T2] + tSpan1=[4,T1] + tSpan2=[T1,T2] + +#solve model 1: ODE +t, x= FDEsolver(F, [4,T2], x0, ones(8), par, nc=4,h=.01) + S=x[:,1]; E=x[:,2]; I=x[:,3]; R=x[:,4]; + L=x[:,5]; H=x[:,6]; B=x[:,7]; C=x[:,8]; + N=sum(x,dims=2) + μ3=mean([μ1 .+ μ2 .* N, α1 .- α2 .*N]) + Appx1=@.I+R+L+H+B+C-μ3*(N-S-E) + Err1=rmsd(Appx1[1:10:end,:], Data2(4:.1:T2)) # Root-mean-square error +print("RMSD(model1)=$Err1") + +#solve Eq constant N: model 2 +_, x = FDEsolver(Fc, tSpan, x0, ones(8), parC,h=.01,nc=4) + Sc=x[:,1]; Ec=x[:,2]; Ic=x[:,3]; Rc=x[:,4]; + Lc=x[:,5]; Hc=x[:,6]; Bc=x[:,7]; Cc=x[:,8]; + Nc=sum(x0) + AppxC=@.Ic+Rc+Lc+Hc+Bc+Cc-μ*(Nc-Sc-Ec) + ErrC=rmsd(AppxC[1:10:end,:], Data2(4:.1:T2)) +print("RMSD(model2)=$ErrC") + +# model3: all time fractional +αf=[0.9414876354377308, 0.9999999999997804, 0.999999999999543, 0.9999999999998147, 0.9999999999806818, 0.9999999999981127, 0.9430213976522912, 0.7932329945305198] + +t, x= FDEsolver(F, tSpan, x0, αf, par, nc=4,h=.01) + Sf=x[:,1]; Ef=x[:,2]; If=x[:,3]; Rf=x[:,4]; + Lf=x[:,5]; Hf=x[:,6]; Bf=x[:,7]; Cf=x[:,8]; + Nf=sum(x,dims=2) + μ3=mean([μ1 .+ μ2 .* Nf, α1 .- α2 .*Nf]) + AppxF=@.If+Rf+Lf+Hf+Bf+Cf-μ3*(Nf-Sf-Ef) + ErrF=rmsd(AppxF[1:10:end,:], Data2(4:.1:T2)) # +print("RMSD(model3)=$ErrF") + +# model4 using 2 sections: first integer until 250 then fractional +#optimized order of derivatives +# α=[0.9999826367080396, 0.583200182055744, 0.9999997090382831, 0.9999552779111426, 0.9999931233485699, 0.999931061595331, 0.9999999960137106, 0.9999997372162407] +α=ones(8);α[2]=0.583200182055744 +t1, x1= FDEsolver(F, tSpan1, x0, ones(8), par, nc=4,h=.01) + x02=x1[end,:] + t2, x2= FDEsolver(F, tSpan2, x02, α, par, nc=4,h=.01) + xf=vcat(x1,x2[2:end,:]) + t=vcat(t1,t2[2:end,:]) + Sf2=xf[:,1]; Ef2=xf[:,2]; If2=xf[:,3]; Rf2=xf[:,4]; + Lf2=xf[:,5]; Hf2=xf[:,6]; Bf2=xf[:,7]; Cf2=xf[:,8]; + Nf2=sum(xf,dims=2) + μ3=mean([μ1 .+ μ2 .* Nf2, α1 .- α2 .*Nf2]) + AppxF2=@.If2+Rf2+Lf2+Hf2+Bf2+Cf2-μ3*(Nf2-Sf2-Ef2) + ErrF2=rmsd(AppxF2[1:10:end,:], Data2(4:.1:T2)) # root-mean-square error +print("RMSD(model4)=$ErrF2") +# Norm/abs +# Mf8=@.abs(Appx[1:10:end,:] - Data2(4:.1:T2)) +# M1=@.abs(Appx1[1:10:end,:] - Data2(4:.1:T2)) + + +#solutions +_, x = FDEsolver(Fc, [4,450], x0, ones(8), parC,h=.01,nc=4) + Sc=x[:,1]; Ec=x[:,2]; Ic=x[:,3]; Rc=x[:,4]; + Lc=x[:,5]; Hc=x[:,6]; Bc=x[:,7]; Cc=x[:,8]; + Nc=sum(x0) + AppxC=@.Ic+Rc+Lc+Hc+Bc+Cc-μ*(Nc-Sc-Ec) + +t, x= FDEsolver(F, [4,450], x0, ones(8), par, nc=4,h=.01) + S=x[:,1]; E=x[:,2]; I=x[:,3]; R=x[:,4]; + L=x[:,5]; H=x[:,6]; B=x[:,7]; C=x[:,8]; + N=sum(x,dims=2) + μ3=mean([μ1 .+ μ2 .* N, α1 .- α2 .*N]) + Appx1=@.I+R+L+H+B+C-μ3*(N-S-E) + +t, x= FDEsolver(F, [4,450], x0, αf, par, nc=4,h=.01) + Sf=x[:,1]; Ef=x[:,2]; If=x[:,3]; Rf=x[:,4]; + Lf=x[:,5]; Hf=x[:,6]; Bf=x[:,7]; Cf=x[:,8]; + Nf=sum(x,dims=2) + μ3=mean([μ1 .+ μ2 .* Nf, α1 .- α2 .*Nf]) + AppxF=@.If+Rf+Lf+Hf+Bf+Cf-μ3*(Nf-Sf-Ef) + +t1, x1= FDEsolver(F, tSpan1, x0, ones(8), par, nc=4,h=.01) + x02=x1[end,:] + t2, x2= FDEsolver(F, [T1,450], x02, α, par, nc=4,h=.01) + xf=vcat(x1,x2[2:end,:]) + t=vcat(t1,t2[2:end,:]) + Sf2=xf[:,1]; Ef2=xf[:,2]; If2=xf[:,3]; Rf2=xf[:,4]; + Lf2=xf[:,5]; Hf2=xf[:,6]; Bf2=xf[:,7]; Cf2=xf[:,8]; + Nf2=sum(xf,dims=2) + μ3=mean([μ1 .+ μ2 .* Nf2, α1 .- α2 .*Nf2]) + AppxF2=@.If2+Rf2+Lf2+Hf2+Bf2+Cf2-μ3*(Nf2-Sf2-Ef2) + +##plot dynamics +gr() +scatter(data[:,1],data[:,2], label="Real data", c="khaki3",markerstrokewidth=0) + plot!(t,AppxC,ylabel="Cumulative confirmed cases",lw=2, label="Model1, RMSD=2085", c="darkorange3", linestyle=:dashdot) + plot!(t,Appx1,lw=2, label="Model2, RMSD=276.5", c="deepskyblue1") + plot!(t,AppxF,lw=2, label="Model3, RMSD=245.9", c="deeppink", linestyle=:dash) + plotModels=plot!(t,AppxF2,xlabel="Days", legendposition=:bottomright, lw=2, label="Model4, RMSD=211.0", linestyle=:dot, c="blue1", + title = "(a)", titleloc = :left, titlefont = font(10)) + +#plot population +plot(t,Nc*ones(length(t)),c="darkorange3",lw=2,legendposition=:right,label="Model1", linestyle=:dashdot) + plot!(t,N,lw=2,legendposition=:right,label="Model2", c="deepskyblue1") + plot!(t,Nf,lw=2,legendposition=:right,label="Model3",linestyle=:dash,c="deeppink") + plotN=plot!(t,Nf2,label="Model4",linestyle=:dot, lw=2, xlabel="Days", ylabel="Variable population",c="blue1", + title = "(b)", titleloc = :left, titlefont = font(10)) + + + +##plot heatmaps + +#optimized order of derivatives +tSpan1=[4,438] + +Err=zeros(8,40) +Errr=zeros(8,40) +#FDE: model 3 +for i=1:8 + print(i)# to check in which order is counting + for j=1:40 + print(j)# to check in which itteration is counting + α=ones(8) + α[i]=1-j*.01 + t, x= FDEsolver(F, [4,438], x0, α, par, nc=2,h=.01) + S1=x[:,1]; E1=x[:,2]; I1=x[:,3]; R1=x[:,4]; + L1=x[:,5]; H1=x[:,6]; B1=x[:,7]; C1=x[:,8]; + N1=sum(x,dims=2) + μ3=mean([μ1 .+ μ2 .* N1, α1 .- α2 .*N1]) + Appx1=@.I1+R1+L1+H1+B1+C1-μ3*(N1-S1-E1) + Err[i,j]=rmsd(Appx1[1:10:end,:], Data2(4:.1:438)) # Normalized root-mean-square error + end +end + +#ODE+FDE: model 4 +for i=1:8 + print(i)# to check in which order is counting + for j=1:40 + print(j)# to check in which itteration is counting + α=ones(8) + t1, x1= FDEsolver(F, [4,250], x0, α, par, nc=2,h=.01) + x02=x1[end,:] + α[i]=1-j*.01 + t2, x2= FDEsolver(F, [250,438], x02, α, par, nc=2,h=.01) + x=vcat(x1,x2[2:end,:]) + t=vcat(t1,t2[2:end,:]) + S1=x[:,1]; E1=x[:,2]; I1=x[:,3]; R1=x[:,4]; + L1=x[:,5]; H1=x[:,6]; B1=x[:,7]; C1=x[:,8]; + N1=sum(x,dims=2) + μ3=mean([μ1 .+ μ2 .* N1, α1 .- α2 .*N1]) + Appx1=@.I1+R1+L1+H1+B1+C1-μ3*(N1-S1-E1) + Errr[i,j]=rmsd(Appx1[1:10:end,:], Data2(4:.1:438)) # Normalized root-mean-square error + end +end + +RMSD1= log10(Err1) #error for integer order + +# plot heatmap model 4 +CoLor = cgrad([:yellow2,:darkgoldenrod1,:mediumpurple], [.45,.4], categorical = true) # define color +heatmap(["S","E","I","R","L","H","B","C"],0.6:.01:.99,log10.(Errr[:,end:-1:1]'), + color=CoLor,colorbar_title="log10(RMSD)",clim=(RMSD1-.18,RMSD1+.21), + xlabel="Individual classes", ylabel="Order of derivative", + title = "(e) Model4 ", titleloc = :left, titlefont = font(10)) + pltheat=plot!(Shape(0 .+ [1,2,2,1], 0 .+ [.595,.595,.995,.995]), fillcolor=:false, legend=:false) + + +# plot heatmap model 4 +pltheat1=heatmap(["S","E","I","R","L","H","B","C"],0.6:.01:.99,log10.(Err[:,end:-1:1]'), + color=CoLor,colorbar_title="log10(RMSD)",clim=(RMSD1-.18,RMSD1+.21), + xlabel="Individual classes", ylabel="Order of derivative", + title = "(d) Model3 ", titleloc = :left, titlefont = font(10), colorbar=:false) + + +# savefig(pltheat,"pltheat.svg") +# savefig(pltheat1,"pltheat1.svg") + + +tt, x= FDEsolver(F, [4,438], x0, ones(8), par, nc=2,h=.01) + S1=x[:,1]; E1=x[:,2]; I1=x[:,3]; R1=x[:,4]; + L1=x[:,5]; H1=x[:,6]; B1=x[:,7]; C1=x[:,8]; + N1=sum(x,dims=2) + μ3=mean([μ1 .+ μ2 .* N1, α1 .- α2 .*N1]) + Appx1=@.I1+R1+L1+H1+B1+C1-μ3*(N1-S1-E1) + # rmsd(Appx1[1:10:end,:], Data2(4:.1:438)) + + +#plot +plot() +N=30 +err=zeros(N) +for i=1:N + α=ones(8) + t1, x1= FDEsolver(F, [4,250], x0, α, par, nc=2,h=.01) + x02=x1[end,:] + α[2]=1-i*.02 + t2, x2= FDEsolver(F, [250,438], x02, α, par, nc=2,h=.01) + x=vcat(x1,x2[2:end,:]) + t=vcat(t1,t2[2:end,:]) + S1=x[:,1]; E1=x[:,2]; I1=x[:,3]; R1=x[:,4]; + L1=x[:,5]; H1=x[:,6]; B1=x[:,7]; C1=x[:,8]; + N1=sum(x,dims=2) + μ3=mean([μ1 .+ μ2 .* N1, α1 .- α2 .*N1]) + Appx=@.I1+R1+L1+H1+B1+C1-μ3*(N1-S1-E1) + # err[i]=copy(rmsd(Appx[1:10:end,:], Data2(4:.1:438))) + plot!(t[25000:100:end],Appx[25000:100:end,:],lw=2; alpha=0.2, color="blue1") + +end + +scatter!(data[74:end,1],data[74:end,2], label="Real data", c="khaki3",markerstrokewidth=0) +pltFtry=plot!(tt[25000:50:end],Appx1[25000:50:end,:],xlabel="Days",ylabel="Cumulative confirmed cases", + lw=2, c="chocolate1",label="Integer model",legend=:false, + title = "(c)", titleloc = :left, titlefont = font(10)) + +# savefig(pltFtry,"pltFtry.svg") + + +l = @layout [a{0.5h}; [grid(1,2)]; [grid(1,1) b{0.6w}]] +Allplots=plot(plotModels, plotN, pltFtry, pltheat1, pltheat, layout= l, size = (900, 900)) + +savefig(Allplots,"Allplots.svg") diff --git a/Optim-FDE-Data-Order.jl b/Optim-FDE-Data-Order.jl new file mode 100644 index 0000000..cfcd5e9 --- /dev/null +++ b/Optim-FDE-Data-Order.jl @@ -0,0 +1,239 @@ +#for finding optimized orders + +using CSV +using DataFrames +using FdeSolver, Plots,SpecialFunctions, Optim, StatsBase, Random +using Interpolations +# Dataset +Data=CSV.read("datoswho.csv", DataFrame) + +data=(Matrix(Float64.(Data))) +#initial conditons and parameters + +x0=[18000,0,15,0,0,0,0,0]# initial conditons S0,E0,I0,R0,L0,H0,B0,C0 + +α1=35.37e-3 # Density independent part of the birth rate for individuals. +# α2=2.6297211237376283e-7# Density dependent part of the birth rate for individuals. +σ=1/11.4 # Per capita rate at which exposed individuals become infectious. +γ1=0.1 # Per capita rate of progression of individuals from the infectious class to the asymptomatic class. +γ2=1/5 # Per capita rate of progression of individuals from the hospitalized class to the asymptomatic class. +γ3=1/30 # Per capita recovery rate of individuals from the asymptomatic class to the complete recovered class. +ϵ=1/9.6 # Fatality rate. +δ1=1/2 # Per capita rate of progression of individuals from the dead class to the buried class. +δ2=1/4.6 # Per capita rate of progression of individuals from the hospitalized class to the buried class. +τ=1/5 # Per capita rate of progression of individuals from the infectious class to the hospitalized class. +βi=0.14 # Contact rate of infective individuals and susceptible. +βd=0.4 # Contact rate of infective individuals and dead. +βh=0.29 # Contact rate of infective individuals and hospitalized. +# βr=0.185 # Contact rate of infective individuals and asymptomatic. +μ=14e-3 +μ1=10.17e-3 # Density independent part of the death rate for individuals. +# μ2=3.750689119502623e-7 # Density dependent part of the death rate for individuals. +ξ=14e-3 # Incineration rate + + +α2=8.334570273681649e-8 +μ2=4.859965522407347e-7 +βr=0.2174851498937417 + +par=[α1,σ,γ1,γ2,γ3,ϵ,δ1,δ2,τ,βi,βd,βh,βr,μ1,ξ,α2,μ2] +# Par=vcat(par,α2,μ2) +# α=ones(8) # order of derivatives +T=215 +tSpan=[4,T] + +#Define the equation +function F(t, x, par) + + α1,σ,γ1,γ2,γ3,ϵ,δ1,δ2,τ,βi,βd,βh,βr,μ1,ξ,α2,μ2=par + S, E, I1, R, L, H, B, C=x + N=sum(x) + + dS=(α1-α2*N)*N - βi/N*S*I1 - βh/N*S*H - βd/N*S*L - βr/N*S*R - (μ1+μ2*N)*S + dE=βi/N*S*I1 + βh/N*S*H + βd/N*S*L + βr/N*S*R - σ*E - (μ1+μ2*N)*E + dI=σ*E - (γ1 + ϵ + τ + μ1 + μ2*N)*I1 + dR=γ1*I1 + γ2*H - (γ3 + μ1 + μ2*N)*R + dL=ϵ*I1 - (δ1+ξ)*L + dH=τ*I1 - (γ2 + δ2 + μ1 + μ2*N)*H + dB=δ1*L + δ2*H - ξ*B + dC=γ3*R - (μ1 + μ2*N)*C + + return [dS, dE, dI, dR, dL, dH, dB, dC] + +end + +## optimazation of μ2 and α2 for integer order model +Data2=LinearInterpolation(data[:,1], data[:,2]) + +function loss_1(p)# loss function + # α2, μ2=p[1:2] + # Par=vcat(par,α2,μ2) + # α=p[3:10] + α=p + if size(x0,2) != Int64(ceil(maximum(α))) # to prevent any errors regarding orders higher than 1 + indx=findall(x-> x>1, α) + α[indx]=ones(length(indx)) + end + _, x = FDEsolver(F, tSpan, x0, α, par,h=.01,nc=4) + S=x[:,1]; E=x[:,2]; I1=x[:,3]; R=x[:,4]; + L=x[:,5]; H=x[:,6]; B=x[:,7]; C=x[:,8]; + N=sum(x,dims=2) + μ3=mean([μ1 .+ μ2 .* N, α1 .- α2 .*N]) + Appx=@.I1+R+L+H+B+C-μ3*(N-S-E) + rmsd(Appx[1:10:end,:], Data2(4:.1:T); normalize=:true) # Normalized root-mean-square error +end + +p_lo_1=[.5, .5, .5, .5, .5, .5, .5, .5] #lower bound +p_up_1=[1, 1, 1, 1, 1, 1, 1, 1] # upper bound +p_vec_1=[0.939531849631367, 0.9998860578731849, 0.9997741401361016, 0.9999076787478497, 0.991141852453178, 0.9990686884413122, 0.9432066029272619, 0.7936761619660204] +Res1=optimize(loss_1,p_lo_1,p_up_1,p_vec_1,Fminbox(LBFGS()),# Broyden–Fletcher–Goldfarb–Shanno algorithm +# Res1=optimize(loss_1,p_lo_1,p_up_1,p_vec_1,SAMIN(rt=.99), # Simulated Annealing algorithm (sometimes it has better perfomance than (L-)BFGS) + Optim.Options(outer_iterations = 10, + iterations=300, + show_trace=true, + show_every=1) + ) + +α=vcat(Optim.minimizer(Res1)) +α=[1,.9958980293056472,0.8623938517464153,1,1,1,1,1] + +T2=250 +tSpan2=[200,T2] +_, xx = FDEsolver(F, [4,200], x0, α, par,h=.01,nc=4) +x02=xx[end,:] + + +function loss_2(p)# loss function + # α2, μ2=p[1:2] + # Par=vcat(par,α2,μ2) + # α=p[3:10] + α=p + if size(x0,2) != Int64(ceil(maximum(α))) # to prevent any errors regarding orders higher than 1 + indx=findall(x-> x>1, α) + α[indx]=ones(length(indx)) + end + _, x = FDEsolver(F, tSpan2, x02, α, par,h=.01,nc=4) + S=x[:,1]; E=x[:,2]; I1=x[:,3]; R=x[:,4]; + L=x[:,5]; H=x[:,6]; B=x[:,7]; C=x[:,8]; + N=sum(x,dims=2) + μ3=mean([μ1 .+ μ2 .* N, α1 .- α2 .*N]) + Appx=@.I1+R+L+H+B+C-μ3*(N-S-E) + rmsd(Appx[1:10:end,:], Data2(200:.1:T2); normalize=:true) # Normalized root-mean-square error +end + +p_lo_1=[.5, .5, .5, .5, .5, .5, .5, .5] #lower bound +p_up_1=[1, 1, 1, 1, 1, 1, 1, 1] # upper bound +p_vec_1=[0.939531849631367, 0.9998860578731849, 0.9997741401361016, 0.9999076787478497, 0.991141852453178, 0.9990686884413122, 0.9432066029272619, 0.7936761619660204] +Res2=optimize(loss_2,p_lo_1,p_up_1,p_vec_1,Fminbox(LBFGS()),# Broyden–Fletcher–Goldfarb–Shanno algorithm +# Res1=optimize(loss_1,p_lo_1,p_up_1,p_vec_1,SAMIN(rt=.99), # Simulated Annealing algorithm (sometimes it has better perfomance than (L-)BFGS) + Optim.Options(outer_iterations = 5, + iterations=30, + show_trace=true, + show_every=1) + ) + +αα=vcat(Optim.minimizer(Res2)) +αα=[0.9624227617441116, 0.9999999994058398, 0.9999999981457667, 0.9999999999999999, 0.9999954941897224, 0.9999999850132952, 0.9999999999465237, 0.9235549807693278] + + +T3=438 +tSpan3=[T2,T3] +# _, xxx = FDEsolver(F, tSpan2, x02, αα, par,h=.01,nc=4) +# x03=xxx[end,:] +_, xxx = FDEsolver(F, [0,250], x0, ones(8), par,h=.01,nc=4) +x03=xxx[end,:] + +function loss_3(p)# loss function + # α2, μ2=p[1:2] + # Par=vcat(par,α2,μ2) + # α=p[3:10] + α=p + if size(x0,2) != Int64(ceil(maximum(α))) # to prevent any errors regarding orders higher than 1 + indx=findall(x-> x>1, α) + α[indx]=ones(length(indx)) + end + _, x = FDEsolver(F, tSpan3, x03, α, par,h=.01,nc=4) + S=x[:,1]; E=x[:,2]; I1=x[:,3]; R=x[:,4]; + L=x[:,5]; H=x[:,6]; B=x[:,7]; C=x[:,8]; + N=sum(x,dims=2) + μ3=mean([μ1 .+ μ2 .* N, α1 .- α2 .*N]) + Appx=@.I1+R+L+H+B+C-μ3*(N-S-E) + rmsd(Appx[1:10:end,:], Data2(T2:.1:T3); normalize=:true) # Normalized root-mean-square error +end + +p_lo_1=[.5, .5, .5, .5, .5, .5, .5, .5] #lower bound +p_up_1=[1, 1, 1, 1, 1, 1, 1, 1] # upper bound +p_vec_1=[0.9999826367080396, 0.583200182055744, 0.9999997090382831, 0.9999552779111426, 0.9999931233485699, 0.999931061595331, 0.9999999960137106, 0.9999997372162407] +Res3=optimize(loss_3,p_lo_1,p_up_1,p_vec_1,Fminbox(LBFGS()),# Broyden–Fletcher–Goldfarb–Shanno algorithm +# Res1=optimize(loss_1,p_lo_1,p_up_1,p_vec_1,SAMIN(rt=.99), # Simulated Annealing algorithm (sometimes it has better perfomance than (L-)BFGS) + Optim.Options(outer_iterations = 5, + iterations=30, + show_trace=true, + show_every=1) + ) +ααα=vcat(Optim.minimizer(Res3)) +ααα=[0.9999826367080396, 0.583200182055744, 0.9999997090382831, 0.9999552779111426, 0.9999931233485699, 0.999931061595331, 0.9999999960137106, 0.9999997372162407] +ααα=ones(8);ααα[2]=0.583200182055744 + +#plotting + +#3 section of fractional orders +t1, x1= FDEsolver(F, [4,200], x0, α, par, nc=4,h=.01) +x02=x1[end,:] +t2, x2= FDEsolver(F, [200,250], x02, αα, par, nc=4,h=.01) +x03=x2[end,:] +t3, x3= FDEsolver(F, [T2,T3], x03, ααα, par, nc=4,h=.01) +x=vcat(x1,x2[2:end,:],x3[2:end,:]) +t=vcat(t1,t2[2:end,:],t3[2:end,:]) +S=x[:,1]; E=x[:,2]; I1=x[:,3]; R=x[:,4]; +L=x[:,5]; H=x[:,6]; B=x[:,7]; C=x[:,8]; +N=sum(x,dims=2) +μ3=mean([μ1 .+ μ2 .* N, α1 .- α2 .*N]) +Appx=@.I1+R+L+H+B+C-μ3*(N-S-E) +Err=rmsd(Appx[1:10:end,:], Data2(4:.1:T3); normalize=:true) # Normalized root-mean-square error +plot(t,Appx,xlabel="time") +scatter!(data[:,1],data[:,2]) + +# 2 sections: first integer until 250 then fractional +t1, x1= FDEsolver(F, [4,250], x0, ones(8), par, nc=4,h=.01) +x02=x1[end,:] +t2, x2= FDEsolver(F, tSpan3, x02, ααα, par, nc=4,h=.01) +xf=vcat(x1,x2[2:end,:]) +t=vcat(t1,t2[2:end,:]) +S=xf[:,1]; E=xf[:,2]; I=xf[:,3]; R=xf[:,4]; +L=xf[:,5]; H=xf[:,6]; B=xf[:,7]; C=xf[:,8]; +N=sum(xf,dims=2) +μ3=mean([μ1 .+ μ2 .* N, α1 .- α2 .*N]) +Appx=@.I+R+L+H+B+C-μ3*(N-S-E) +Err=rmsd(Appx[1:10:end,:], Data2(4:.1:T3); normalize=:false) # Normalized root-mean-square error +scatter(data[:,1],data[:,2]) +plot!(t,Appx,xlabel="time",lw=2) + +#ODE +t, x= FDEsolver(F, [4,T3], x0, ones(8), par, nc=4,h=.01) +S1=x[:,1]; E1=x[:,2]; I1=x[:,3]; R1=x[:,4]; +L1=x[:,5]; H1=x[:,6]; B1=x[:,7]; C1=x[:,8]; +N1=sum(x,dims=2) +μ3=mean([μ1 .+ μ2 .* N1, α1 .- α2 .*N1]) +Appx1=@.I1+R1+L1+H1+B1+C1-μ3*(N1-S1-E1) +Err1=rmsd(Appx1[1:10:end,:], Data2(4:.1:T3); normalize=:false) # Normalized root-mean-square error +plot!(t,Appx1,xlabel="time", legendposition=:right, lw=2) +# scatter!(data[:,1],data[:,2]) + +#plot population +plot(t,N) +plot!(t,N1,legendposition=:right) + +plot(t,S) +plot!(t,S1,legendposition=:right) + +plot(t,xf) +plot!(t,x, legend=:false) + +# Norm plot +using LinearAlgebra +Mf8=@.abs(Appx[1:10:end,:] - Data2(4:.1:T)) +M1=@.abs(Appx1[1:10:end,:] - Data2(4:.1:T)) + +plot(t[1:10:end,:],Mf8) +plot!(t[1:10:end,:],M1) diff --git a/PLOT-Code.jl b/PLOT-Code.jl deleted file mode 100644 index 66117bb..0000000 --- a/PLOT-Code.jl +++ /dev/null @@ -1,141 +0,0 @@ -using FdeSolver, Plots, Statistics -using CSV -using DataFrames,StatsBase -using Interpolations,LinearAlgebra - -#initial conditons and parameters - -x0=[18000,0,15,0,0,0,0,0]# initial conditons S0,E0,I0,R0,L0,H0,B0,C0 - -α1=35.37e-3 # Density independent part of the birth rate for individuals. -σ=1/11.4 # Per capita rate at which exposed individuals become infectious. -γ1=0.1 # Per capita rate of progression of individuals from the infectious class to the asymptomatic class. -γ2=1/5 # Per capita rate of progression of individuals from the hospitalized class to the asymptomatic class. -γ3=1/30 # Per capita recovery rate of individuals from the asymptomatic class to the complete recovered class. -ϵ=1/9.6 # Fatality rate. -δ1=1/2 # Per capita rate of progression of individuals from the dead class to the buried class. -δ2=1/4.6 # Per capita rate of progression of individuals from the hospitalized class to the buried class. -τ=1/5 # Per capita rate of progression of individuals from the infectious class to the hospitalized class. -βi=0.14 # Contact rate of infective individuals and susceptible. -βd=0.4 # Contact rate of infective individuals and dead. -βh=0.29 # Contact rate of infective individuals and hospitalized. -μ1=10.17e-3 # Density independent part of the death rate for individuals. -ξ=14e-3 # Incineration rate -μ=14e-3 -#fitted 4 parameters by Turing (HMC) -α2=4.596016101361542e-7 -μ2=4.286684252177388e-7 -βr=0.11565848502778753 -γ3=0.00015462168178164486 - -par=[α1,α2,σ,γ1,γ2,γ3,ϵ,δ1,δ2,τ,βi,βd,βh,βr,μ1,μ2,ξ] -# fitted orders by Optim (LBFGS) -# α=[0.9999999999547146, 0.9999999869785005, 0.9494471201561786, 0.9999999999534154, 0.9504717277935494, 0.9999999852655446, 0.9999999749309598, 0.89616744140632] -α=[1, 1, 0.9494471201561786, 1, 0.9504717277935494, 1, 1, 0.89616744140632] -# α=[1,1,.9,.9,.9,.9,.9,.9] # order of derivatives -T=438 -tSpan=[4,T] - -#Define the equation - -function F(t, x, par) - - α1,α2,σ,γ1,γ2,γ3,ϵ,δ1,δ2,τ,βi,βd,βh,βr,μ1,μ2,ξ=par - S, E, I, R, L, H, B, C=x - N=sum(x) - - dS=(α1-α2*N)*N - βi/N*S*I - βh/N*S*H - βd/N*S*L - βr/N*S*R - (μ1+μ2*N)*S - dE=βi/N*S*I + βh/N*S*H + βd/N*S*L + βr/N*S*R - σ*E - (μ1+μ2*N)*E - dI=σ*E - (γ1 + ϵ + τ + μ1 + μ2*N)*I - dR=γ1*I + γ2*H - (γ3 + μ1 + μ2*N)*R - dL=ϵ*I - (δ1+ξ)*L - dH=τ*I - (γ2 + δ2 + μ1 + μ2*N)*H - dB=δ1*L + δ2*H - ξ*B - dC=γ3*R - (μ1 + μ2*N)*C - - return [dS, dE, dI, dR, dL, dH, dB, dC] - -end -#eq with constant N -parC=[α1,σ,γ1,γ2,γ3,ϵ,δ1,δ2,τ,βi,βd,βh,βr,μ1,ξ,sum(x0)] -function Fc(t, x, par) - - α1,σ,γ1,γ2,γ3,ϵ,δ1,δ2,τ,βi,βd,βh,βr,μ1,ξ,N=par - S, E, I, R, L, H, B, C=x - βr=0.185 - γ3=1/30 - μ=14e-3 - - dS=μ*N - βi/N*S*I - βh/N*S*H - βd/N*S*L - βr/N*S*R - μ*S - dE=βi/N*S*I + βh/N*S*H + βd/N*S*L + βr/N*S*R - σ*E - μ*E - dI=σ*E - (γ1 + ϵ + τ + μ)*I - dR=γ1*I + γ2*H - (γ3 + μ)*R - dL=ϵ*I - (δ1+ξ)*L - dH=τ*I - (γ2 + δ2 + μ)*H - dB=δ1*L + δ2*H - ξ*B - dC=γ3*R - μ*C - - return [dS, dE, dI, dR, dL, dH, dB, dC] - -end -#solving - -#solve Eq constant N -_, x = FDEsolver(Fc, tSpan, x0, ones(8), parC,h=.01,nc=4) -S=x[:,1]; E=x[:,2]; I1=x[:,3]; R=x[:,4]; -L=x[:,5]; H=x[:,6]; B=x[:,7]; C=x[:,8]; -N1C=sum(x,dims=2) -AppxC=@.I1+R+L+H+B+C-μ*(N1C-S-E) - - -#model with optimized orders -t, x= FDEsolver(F, tSpan, x0, α, par, nc=4,h=.01) -S=x[:,1]; E=x[:,2]; I1=x[:,3]; R=x[:,4]; -L=x[:,5]; H=x[:,6]; B=x[:,7]; C=x[:,8]; -N=sum(x,dims=2) -μ3=mean([μ1 .+ μ2 .* N, α1 .- α2 .*N]) -app=@.I1+R+L+H+B+C-μ3*(N-S-E) - -#lets test with integer order -_, x = FDEsolver(F, tSpan, x0, ones(8), par,h=.01,nc=4) -S=x[:,1]; E=x[:,2]; I1=x[:,3]; R=x[:,4]; -L=x[:,5]; H=x[:,6]; B=x[:,7]; C=x[:,8]; -N1=sum(x,dims=2) -μ3=mean([μ1 .+ μ2 .* N1, α1 .- α2 .*N1]) -Appx1=@.I1+R+L+H+B+C-μ3*(N1-S-E) - -# Dataset -Data=CSV.read("datoswho.csv", DataFrame) -data=(Matrix(Float64.(Data))) -Data2=LinearInterpolation(data[:,1], data[:,2]) - - -Err=rmsd(app[1:10:end,:], Data2(4:.1:T); normalize=:true) # Normalized root-mean-square error -Err1=rmsd(Appx1[1:10:end,:], Data2(4:.1:T); normalize=:true) # Normalized root-mean-square error -Err0=rmsd(AppxC[1:10:end,:], Data2(4:.1:T); normalize=:true) # Normalized root-mean-square error -#plotting -plt=plot(t,x,xlabel="time",lw = 2, label=["S" "E" "I" "R" "L" "H" "B" "C"], - thickness_scaling =1.2) -savefig(plt,"plt.png") - -scatter(data[:,1],data[:,2],label="Real data",c="black") -plot!(t,Appx1, lw=3,label="Model 1", xlabel="Time (days)", ylabel="Cumulative confirmed cases") -plot!(t,app,lw=3,label="Model f8",linestyle=:dash,xlabel="time",legendposition=:bottomright) -pltData=plot!(t,AppxC, lw=3,label="Model 0") -savefig(pltData,"pltData.png") - -#plot population -plot(t,N,lw=3,label="Model 1",xlabel="Days",ylabel="Population") -pltN=plot!(t,N1,lw=3,linestyle=:dash,label="Model f8",legendposition=:right) -savefig(pltN,"pltN.png") - -# Norm plot -Mf8=(app[1:10:end,:] - Data2(4:.1:T)) -M1=(Appx1[1:10:end,:] - Data2(4:.1:T)) - -plot(t[1:10:end,:],M1,label="Model 1", xlabel="Days",ylabel="Error (appX - exactX)") -pltErr=plot!(t[1:10:end,:],Mf8,label="Model f8") -savefig(pltErr,"pltErr.png") - -norm(Mf8,2) -norm(M1,2) diff --git a/PLOT-Code1.jl b/PLOT-Code1.jl deleted file mode 100644 index 7fcc844..0000000 --- a/PLOT-Code1.jl +++ /dev/null @@ -1,139 +0,0 @@ -using FdeSolver, Plots, Statistics -using CSV -using DataFrames,StatsBase -using Interpolations,LinearAlgebra - -#initial conditons and parameters - -x0=[18000,0,15,0,0,0,0,0]# initial conditons S0,E0,I0,R0,L0,H0,B0,C0 - -α1=35.37e-3 # Density independent part of the birth rate for individuals. -σ=1/11.4 # Per capita rate at which exposed individuals become infectious. -γ1=0.1 # Per capita rate of progression of individuals from the infectious class to the asymptomatic class. -γ2=1/5 # Per capita rate of progression of individuals from the hospitalized class to the asymptomatic class. -γ3=1/30 # Per capita recovery rate of individuals from the asymptomatic class to the complete recovered class. -ϵ=1/9.6 # Fatality rate. -δ1=1/2 # Per capita rate of progression of individuals from the dead class to the buried class. -δ2=1/4.6 # Per capita rate of progression of individuals from the hospitalized class to the buried class. -τ=1/5 # Per capita rate of progression of individuals from the infectious class to the hospitalized class. -βi=0.14 # Contact rate of infective individuals and susceptible. -βd=0.4 # Contact rate of infective individuals and dead. -βh=0.29 # Contact rate of infective individuals and hospitalized. -μ1=10.17e-3 # Density independent part of the death rate for individuals. -ξ=14e-3 # Incineration rate -μ=14e-3 -#fitted 3 parameters by Turing (HMC) -α2=8.334570273681649e-8 -μ2=4.859965522407347e-7 -βr=0.2174851498937417 - -par=[α1,α2,σ,γ1,γ2,γ3,ϵ,δ1,δ2,τ,βi,βd,βh,βr,μ1,μ2,ξ] -#fitted order based on 3 optimized parameters -# α=[0.9414876125410915, 0.999999999635505, 0.999999999265876, 0.9999999997019355, 0.9999999682296614, 0.999999996978037, 0.9430213826071122, 0.793233020972691] -α=[0.9414876354377308, 0.9999999999997804, 0.999999999999543, 0.9999999999998147, 0.9999999999806818, 0.9999999999981127, 0.9430213976522912, 0.7932329945305198] -T=438 -tSpan=[4,T] - -#Define the equation - -function F(t, x, par) - - α1,α2,σ,γ1,γ2,γ3,ϵ,δ1,δ2,τ,βi,βd,βh,βr,μ1,μ2,ξ=par - S, E, I, R, L, H, B, C=x - N=sum(x) - - dS=(α1-α2*N)*N - βi/N*S*I - βh/N*S*H - βd/N*S*L - βr/N*S*R - (μ1+μ2*N)*S - dE=βi/N*S*I + βh/N*S*H + βd/N*S*L + βr/N*S*R - σ*E - (μ1+μ2*N)*E - dI=σ*E - (γ1 + ϵ + τ + μ1 + μ2*N)*I - dR=γ1*I + γ2*H - (γ3 + μ1 + μ2*N)*R - dL=ϵ*I - (δ1+ξ)*L - dH=τ*I - (γ2 + δ2 + μ1 + μ2*N)*H - dB=δ1*L + δ2*H - ξ*B - dC=γ3*R - (μ1 + μ2*N)*C - - return [dS, dE, dI, dR, dL, dH, dB, dC] - -end -#eq with constant N -parC=[α1,σ,γ1,γ2,γ3,ϵ,δ1,δ2,τ,βi,βd,βh,βr,μ1,ξ,sum(x0)] -function Fc(t, x, par) - - α1,σ,γ1,γ2,γ3,ϵ,δ1,δ2,τ,βi,βd,βh,βr,μ1,ξ,N=par - S, E, I, R, L, H, B, C=x - βr=0.185 - γ3=1/30 - μ=14e-3 - - dS=μ*N - βi/N*S*I - βh/N*S*H - βd/N*S*L - βr/N*S*R - μ*S - dE=βi/N*S*I + βh/N*S*H + βd/N*S*L + βr/N*S*R - σ*E - μ*E - dI=σ*E - (γ1 + ϵ + τ + μ)*I - dR=γ1*I + γ2*H - (γ3 + μ)*R - dL=ϵ*I - (δ1+ξ)*L - dH=τ*I - (γ2 + δ2 + μ)*H - dB=δ1*L + δ2*H - ξ*B - dC=γ3*R - μ*C - - return [dS, dE, dI, dR, dL, dH, dB, dC] - -end -#solving - -#solve Eq constant N -_, x = FDEsolver(Fc, tSpan, x0, ones(8), parC,h=.01,nc=4) -S=x[:,1]; E=x[:,2]; I1=x[:,3]; R=x[:,4]; -L=x[:,5]; H=x[:,6]; B=x[:,7]; C=x[:,8]; -N1C=sum(x,dims=2) -AppxC=@.I1+R+L+H+B+C-μ*(N1C-S-E) - - -#model with optimized orders -t, x= FDEsolver(F, tSpan, x0, α, par, nc=4,h=.01) -S=x[:,1]; E=x[:,2]; I1=x[:,3]; R=x[:,4]; -L=x[:,5]; H=x[:,6]; B=x[:,7]; C=x[:,8]; -N=sum(x,dims=2) -μ3=mean([μ1 .+ μ2 .* N, α1 .- α2 .*N]) -app=@.I1+R+L+H+B+C-μ3*(N-S-E) - -#lets test with integer order -_, x = FDEsolver(F, tSpan, x0, ones(8), par,h=.01,nc=4) -S=x[:,1]; E=x[:,2]; I1=x[:,3]; R=x[:,4]; -L=x[:,5]; H=x[:,6]; B=x[:,7]; C=x[:,8]; -N1=sum(x,dims=2) -μ3=mean([μ1 .+ μ2 .* N1, α1 .- α2 .*N1]) -Appx1=@.I1+R+L+H+B+C-μ3*(N1-S-E) - -# Dataset -Data=CSV.read("datoswho.csv", DataFrame) -data=(Matrix(Float64.(Data))) -Data2=LinearInterpolation(data[:,1], data[:,2]) - - -Err=rmsd(app[1:10:end,:], Data2(4:.1:T); normalize=:true) # Normalized root-mean-square error -Err1=rmsd(Appx1[1:10:end,:], Data2(4:.1:T); normalize=:true) # Normalized root-mean-square error -Err0=rmsd(AppxC[1:10:end,:], Data2(4:.1:T); normalize=:true) # Normalized root-mean-square error -#plotting -plt=plot(t,x,xlabel="time",lw = 2, label=["S" "E" "I" "R" "L" "H" "B" "C"], - thickness_scaling =1.2) -savefig(plt,"plt.png") - -scatter(data[:,1],data[:,2],label="Real data",c="black") -plot!(t,Appx1, lw=3,label="Model 1", xlabel="Time (days)", ylabel="Cumulative confirmed cases") -plot!(t,app,lw=3,label="Model f8",linestyle=:dash,xlabel="time",legendposition=:bottomright) -pltData=plot!(t,AppxC, lw=3,label="Model 0") -savefig(pltData,"pltData.png") - -#plot population -plot(t,N,lw=3,label="Model 1",xlabel="Days",ylabel="Population") -pltN=plot!(t,N1,lw=3,linestyle=:dash,label="Model f8",legendposition=:right) -savefig(pltN,"pltN.png") - -# Norm plot -Mf8=(app[1:10:end,:] - Data2(4:.1:T)) -M1=(Appx1[1:10:end,:] - Data2(4:.1:T)) - -plot(t[1:10:end,:],M1,label="Model 1", xlabel="Days",ylabel="Error (appX - exactX)") -pltErr=plot!(t[1:10:end,:],Mf8,label="Model f8") -savefig(pltErr,"pltErr.png") - -norm(Mf8,2) -norm(M1,2) diff --git a/PLOT-Code2.jl b/PLOT-Code2.jl deleted file mode 100644 index f258fd5..0000000 --- a/PLOT-Code2.jl +++ /dev/null @@ -1,120 +0,0 @@ -using CSV -using DataFrames -using FdeSolver, Plots,SpecialFunctions, Optim, StatsBase, Random -using Interpolations, LinearAlgebra -# Dataset -Data=CSV.read("datoswho.csv", DataFrame) -data=(Matrix(Float64.(Data))) -Data2=LinearInterpolation(data[:,1], data[:,2]) #Interpolate the data -#initial conditons and parameters - -x0=[18000,0,15,0,0,0,0,0]# initial conditons S0,E0,I0,R0,L0,H0,B0,C0 - -α1=35.37e-3 # Density independent part of the birth rate for individuals. -α2=8.334570273681649e-8# Density dependent part of the birth rate for individuals. -σ=1/11.4 # Per capita rate at which exposed individuals become infectious. -γ1=0.1 # Per capita rate of progression of individuals from the infectious class to the asymptomatic class. -γ2=1/5 # Per capita rate of progression of individuals from the hospitalized class to the asymptomatic class. -γ3=1/30 # Per capita recovery rate of individuals from the asymptomatic class to the complete recovered class. -ϵ=1/9.6 # Fatality rate. -δ1=1/2 # Per capita rate of progression of individuals from the dead class to the buried class. -δ2=1/4.6 # Per capita rate of progression of individuals from the hospitalized class to the buried class. -τ=1/5 # Per capita rate of progression of individuals from the infectious class to the hospitalized class. -βi=0.14 # Contact rate of infective individuals and susceptible. -βd=0.4 # Contact rate of infective individuals and dead. -βh=0.29 # Contact rate of infective individuals and hospitalized. -βr=0.2174851498937417 # Contact rate of infective individuals and asymptomatic. -μ=14e-3 -μ1=10.17e-3 # Density independent part of the death rate for individuals. -μ2=4.859965522407347e-7 # Density dependent part of the death rate for individuals. -ξ=14e-3 # Incineration rate - -par=[α1,σ,γ1,γ2,γ3,ϵ,δ1,δ2,τ,βi,βd,βh,βr,μ1,ξ,α2,μ2] - -#Define the equation -function F(t, x, par) - - α1,σ,γ1,γ2,γ3,ϵ,δ1,δ2,τ,βi,βd,βh,βr,μ1,ξ,α2,μ2=par - S, E, I1, R, L, H, B, C=x - N=sum(x) - - dS=(α1-α2*N)*N - βi/N*S*I1 - βh/N*S*H - βd/N*S*L - βr/N*S*R - (μ1+μ2*N)*S - dE=βi/N*S*I1 + βh/N*S*H + βd/N*S*L + βr/N*S*R - σ*E - (μ1+μ2*N)*E - dI=σ*E - (γ1 + ϵ + τ + μ1 + μ2*N)*I1 - dR=γ1*I1 + γ2*H - (γ3 + μ1 + μ2*N)*R - dL=ϵ*I1 - (δ1+ξ)*L - dH=τ*I1 - (γ2 + δ2 + μ1 + μ2*N)*H - dB=δ1*L + δ2*H - ξ*B - dC=γ3*R - (μ1 + μ2*N)*C - - return [dS, dE, dI, dR, dL, dH, dB, dC] - -end - -#optimized order of derivatives -# α=[0.9999826367080396, 0.583200182055744, 0.9999997090382831, 0.9999552779111426, 0.9999931233485699, 0.999931061595331, 0.9999999960137106, 0.9999997372162407] -α=ones(8);α[2]=0.583200182055744 -T1=250 -T2=438 -tSpan1=[4,T1] -tSpan2=[T1,T2] - -#ODE -t, x= FDEsolver(F, [4,T2], x0, ones(8), par, nc=4,h=.01) -S1=x[:,1]; E1=x[:,2]; I1=x[:,3]; R1=x[:,4]; -L1=x[:,5]; H1=x[:,6]; B1=x[:,7]; C1=x[:,8]; -N1=sum(x,dims=2) -μ3=mean([μ1 .+ μ2 .* N1, α1 .- α2 .*N1]) -Appx1=@.I1+R1+L1+H1+B1+C1-μ3*(N1-S1-E1) -Err1=rmsd(Appx1[1:10:end,:], Data2(4:.1:T2)) # Normalized root-mean-square error - -# 2 sections: first integer until 250 then fractional -t1, x1= FDEsolver(F, tSpan1, x0, ones(8), par, nc=4,h=.01) -x02=x1[end,:] -t2, x2= FDEsolver(F, tSpan2, x02, α, par, nc=4,h=.01) -xf=vcat(x1,x2[2:end,:]) -t=vcat(t1,t2[2:end,:]) -S=xf[:,1]; E=xf[:,2]; I=xf[:,3]; R=xf[:,4]; -L=xf[:,5]; H=xf[:,6]; B=xf[:,7]; C=xf[:,8]; -N=sum(xf,dims=2) -μ3=mean([μ1 .+ μ2 .* N, α1 .- α2 .*N]) -Appx=@.I+R+L+H+B+C-μ3*(N-S-E) -Err=rmsd(Appx[1:10:end,:], Data2(4:.1:T2)) # Normalized root-mean-square error - -# Norm/abs -# Mf8=@.abs(Appx[1:10:end,:] - Data2(4:.1:T2)) -# M1=@.abs(Appx1[1:10:end,:] - Data2(4:.1:T2)) - - -#plotting -t, x= FDEsolver(F, [4,450], x0, ones(8), par, nc=4,h=.01) -S1=x[:,1]; E1=x[:,2]; I1=x[:,3]; R1=x[:,4]; -L1=x[:,5]; H1=x[:,6]; B1=x[:,7]; C1=x[:,8]; -N1=sum(x,dims=2) -μ3=mean([μ1 .+ μ2 .* N1, α1 .- α2 .*N1]) -Appx1=@.I1+R1+L1+H1+B1+C1-μ3*(N1-S1-E1) - -t1, x1= FDEsolver(F, tSpan1, x0, ones(8), par, nc=4,h=.01) -x02=x1[end,:] -t2, x2= FDEsolver(F, [T1,450], x02, α, par, nc=4,h=.01) -xf=vcat(x1,x2[2:end,:]) -t=vcat(t1,t2[2:end,:]) -S=xf[:,1]; E=xf[:,2]; I=xf[:,3]; R=xf[:,4]; -L=xf[:,5]; H=xf[:,6]; B=xf[:,7]; C=xf[:,8]; -N=sum(xf,dims=2) -μ3=mean([μ1 .+ μ2 .* N, α1 .- α2 .*N]) -Appx=@.I+R+L+H+B+C-μ3*(N-S-E) - -scatter(data[:,1],data[:,2], label="Real data", c="khaki3",markerstrokewidth=0) -plot!(t,Appx1,ylabel="Cumulative confirmed cases",lw=2, label="Integer model, RMSD=276.5", c="darkorange3") -plot!(t,Appx,xlabel="Days", legendposition=:bottomright, lw=2, label="Fractional model, RMSD=211.0" - , linestyle=:dot, c="blue1") - -#plot population -plot(t,N1,c="darkorange3",lw=4,legendposition=:right,label="Integer model") -plot!(t,N, c="blue1", label="Fractional model",linestyle=:dot, lw=4, xlabel="Days", ylabel="Variable population") - -# -# plot(t[1:10:end,:],Mf8) -# plot!(t[1:10:end,:],M1) -# diff --git a/PLOT-Code3.jl b/PLOT-Code3.jl deleted file mode 100644 index 21a77f7..0000000 --- a/PLOT-Code3.jl +++ /dev/null @@ -1,145 +0,0 @@ -using CSV -using DataFrames -using FdeSolver, Plots,SpecialFunctions, Optim, StatsBase, Random -using Interpolations, LinearAlgebra -# Dataset -Data=CSV.read("datoswho.csv", DataFrame) -data=(Matrix(Float64.(Data))) -Data2=LinearInterpolation(data[:,1], data[:,2]) #Interpolate the data -#initial conditons and parameters - -x0=[18000,0,15,0,0,0,0,0]# initial conditons S0,E0,I0,R0,L0,H0,B0,C0 - -α1=35.37e-3 # Density independent part of the birth rate for individuals. -α2=8.334570273681649e-8# Density dependent part of the birth rate for individuals. -σ=1/11.4 # Per capita rate at which exposed individuals become infectious. -γ1=0.1 # Per capita rate of progression of individuals from the infectious class to the asymptomatic class. -γ2=1/5 # Per capita rate of progression of individuals from the hospitalized class to the asymptomatic class. -γ3=1/30 # Per capita recovery rate of individuals from the asymptomatic class to the complete recovered class. -ϵ=1/9.6 # Fatality rate. -δ1=1/2 # Per capita rate of progression of individuals from the dead class to the buried class. -δ2=1/4.6 # Per capita rate of progression of individuals from the hospitalized class to the buried class. -τ=1/5 # Per capita rate of progression of individuals from the infectious class to the hospitalized class. -βi=0.14 # Contact rate of infective individuals and susceptible. -βd=0.4 # Contact rate of infective individuals and dead. -βh=0.29 # Contact rate of infective individuals and hospitalized. -βr=0.2174851498937417 # Contact rate of infective individuals and asymptomatic. -μ=14e-3 -μ1=10.17e-3 # Density independent part of the death rate for individuals. -μ2=4.859965522407347e-7 # Density dependent part of the death rate for individuals. -ξ=14e-3 # Incineration rate - -par=[α1,σ,γ1,γ2,γ3,ϵ,δ1,δ2,τ,βi,βd,βh,βr,μ1,ξ,α2,μ2] - -#Define the equation -function F(t, x, par) - - α1,σ,γ1,γ2,γ3,ϵ,δ1,δ2,τ,βi,βd,βh,βr,μ1,ξ,α2,μ2=par - S, E, I1, R, L, H, B, C=x - N=sum(x) - - dS=(α1-α2*N)*N - βi/N*S*I1 - βh/N*S*H - βd/N*S*L - βr/N*S*R - (μ1+μ2*N)*S - dE=βi/N*S*I1 + βh/N*S*H + βd/N*S*L + βr/N*S*R - σ*E - (μ1+μ2*N)*E - dI=σ*E - (γ1 + ϵ + τ + μ1 + μ2*N)*I1 - dR=γ1*I1 + γ2*H - (γ3 + μ1 + μ2*N)*R - dL=ϵ*I1 - (δ1+ξ)*L - dH=τ*I1 - (γ2 + δ2 + μ1 + μ2*N)*H - dB=δ1*L + δ2*H - ξ*B - dC=γ3*R - (μ1 + μ2*N)*C - - return [dS, dE, dI, dR, dL, dH, dB, dC] - -end - -#optimized order of derivatives -tSpan1=[4,438] - -Err1=zeros(8,40) -#FDE -for i=1:8 - print(i) - for j=1:40 - print(j) - α=ones(8) - α[i]=1-j*.01 - t, x= FDEsolver(F, [4,438], x0, α, par, nc=2,h=.01) - S1=x[:,1]; E1=x[:,2]; I1=x[:,3]; R1=x[:,4]; - L1=x[:,5]; H1=x[:,6]; B1=x[:,7]; C1=x[:,8]; - N1=sum(x,dims=2) - μ3=mean([μ1 .+ μ2 .* N1, α1 .- α2 .*N1]) - Appx1=@.I1+R1+L1+H1+B1+C1-μ3*(N1-S1-E1) - Err1[i,j]=rmsd(Appx1[1:10:end,:], Data2(4:.1:438)) # Normalized root-mean-square error - end -end - -#ODE+FDE -for i=1:8 - print(i) - for j=1:40 - print(j) - α=ones(8) - t1, x1= FDEsolver(F, [4,250], x0, α, par, nc=2,h=.01) - x02=x1[end,:] - α[i]=1-j*.01 - t2, x2= FDEsolver(F, [250,438], x02, α, par, nc=2,h=.01) - x=vcat(x1,x2[2:end,:]) - t=vcat(t1,t2[2:end,:]) - S1=x[:,1]; E1=x[:,2]; I1=x[:,3]; R1=x[:,4]; - L1=x[:,5]; H1=x[:,6]; B1=x[:,7]; C1=x[:,8]; - N1=sum(x,dims=2) - μ3=mean([μ1 .+ μ2 .* N1, α1 .- α2 .*N1]) - Appx1=@.I1+R1+L1+H1+B1+C1-μ3*(N1-S1-E1) - Err1[i,j]=rmsd(Appx1[1:10:end,:], Data2(4:.1:438)) # Normalized root-mean-square error - end -end - -# Err=Err1 .- 276.5100875995236 - -backend(:plotly) -pltheat=heatmap(["S","E","I","R","L","H","B","C"],0.6:.01:.99,log10.(Err1[:,end:-1:1]'), - c=:balance,colorbar_title="log10(RMSD)", - xlabel="Individual classes", ylabel="Order of derivative") - -savefig(pltheat,"pltheat.svg") - -t, x= FDEsolver(F, [4,438], x0, ones(8), par, nc=2,h=.01) -S1=x[:,1]; E1=x[:,2]; I1=x[:,3]; R1=x[:,4]; -L1=x[:,5]; H1=x[:,6]; B1=x[:,7]; C1=x[:,8]; -N1=sum(x,dims=2) -μ3=mean([μ1 .+ μ2 .* N1, α1 .- α2 .*N1]) -Appx1=@.I1+R1+L1+H1+B1+C1-μ3*(N1-S1-E1) -rmsd(Appx1[1:10:end,:], Data2(4:.1:438)) - - -#plot - - -# plot(; legend=:false) -scatter(data[74:end,1],data[74:end,2], label="Real data", c="khaki3",markerstrokewidth=0) -N=30 -err=zeros(N) -for i=1:N - α=ones(8) - t1, x1= FDEsolver(F, [4,250], x0, α, par, nc=2,h=.01) - x02=x1[end,:] - α[2]=1-i*.02 - t2, x2= FDEsolver(F, [250,438], x02, α, par, nc=2,h=.01) - x=vcat(x1,x2[2:end,:]) - t=vcat(t1,t2[2:end,:]) - S1=x[:,1]; E1=x[:,2]; I1=x[:,3]; R1=x[:,4]; - L1=x[:,5]; H1=x[:,6]; B1=x[:,7]; C1=x[:,8]; - N1=sum(x,dims=2) - μ3=mean([μ1 .+ μ2 .* N1, α1 .- α2 .*N1]) - Appx=@.I1+R1+L1+H1+B1+C1-μ3*(N1-S1-E1) - # err[i]=copy(rmsd(Appx[1:10:end,:], Data2(4:.1:438))) - plot!(t[25000:20:end],Appx[25000:20:end,:],lw=2; alpha=0.2, color="blue1") - -end - -# plot!(t[25000:end],Appx1[25000:end,:],ylabel="Cumulative confirmed cases", - # lw=2, c="darkorange3",label="Integer model",legendposition=:bottomright) -pltFtry=plot!(t[25000:20:end],Appx1[25000:20:end,:],xlabel="Days",ylabel="Cumulative confirmed cases", - lw=2, c="darkorange3",label="Integer model",legend=:false) - - -savefig(pltFtry,"pltFtry.svg") diff --git a/TestNewModel.jl b/TestNewModel.jl new file mode 100644 index 0000000..a1d4530 --- /dev/null +++ b/TestNewModel.jl @@ -0,0 +1,117 @@ +# EBOLA VIRUS DISEASE DYNAMICS WITH SOME +# PREVENTIVE MEASURES: A CASE STUDY +# OF THE 2018–2020 KIVU OUTBREAK + +using FdeSolver, Plots, Statistics +using CSV +using DataFrames +using Interpolations,LinearAlgebra + +#initial conditons and parameters + +x0=[11000000,300,300,300,400,200,5] # initial conditons S0,J0,I0,H0,D0,R0,V0 + +PI=400 +p=0.5708 +τ=0.005 +θ=0.9 +β=0.16 +ν1=1.5267 +ν2=1.5 +ϵ=0.3875 +γj=0.3796 +σ1=0.3074 +σ2=0.7086 +ηj=0.2311 +μ=10.13/1000 +ηi=0.25 +δi=0.3079 +δj=0.0152 +δh=0.1808 +γh=0.4594 +γi=0.0153 +α=0.5723 +b=1/2.01 +u=0.4019 +q1=0.05 +q2=0.05 +q3=0.1173 +q4=0.06 +βv=5.07e-8 + +σh = ϵ*(1 - σ1) +θ1 = θ*τ +σd = ν2*(1 - σ2) +θ2 = μ + θ*τ +φ1 = (α + δj + γj + ηj) +φ3 = (γh + δh) +φ2 = (δi + ηi) + +par=[PI,p,τ,θ,β,ν1,ν2,ϵ,γj,σ1,σ2,ηj,μ,ηi,δi,δj,δh,γh,γi,α,b,u,q1,q2,q3,q4,βv] + +tSpan=[0,50] + +#Define the equation + +function F(t, x, par) + + PI,p,τ,θ,β,ν1,ν2,ϵ,γj,σ1,σ2,ηj,μ,ηi,δi,δj,δh,γh,γi,α,b,u,q1,q2,q3,q4,βv=par + S, J, I, H, D, R, V=copy(x) + N=sum(x[1:6]) + + λ=(β*(J + ϵ*(1-σ1)*H + ν1*I + ν2*(1-σ2)*D))/N+βv*V + + dS = PI - λ*S - θ2*S + dJ= p*λ*S - φ1*J + dI=(1-p)*λ*S + α*J - φ2*I + dH= ηj*J + ηi*I - φ3*H + dD=δj*J + δi*I + δh*H - b*D + dR=θ1*S + γj*J + γi*I + γh*H - μ*R + dV=q1*J + q2*I + q3*H + q4*D - u*V + + return [dS, dJ, dI, dH, dD, dR, dV] + +end + +#solving +t, x= FDEsolver(F, tSpan, x0, ones(7), par, nc=4,h=.01) +plot(t,sum(x[:,2:6],dims=2)) +S=x[:,1]; E=x[:,2]; I1=x[:,3]; R=x[:,4] +L=x[:,5]; H=x[:,6]; B=x[:,7]; C=x[:,8] +N=sum(x,dims=2) +μ3=mean([μ1 .+ μ2 .* N, α1 .- α2 .*N]) +app=@.I1+R+L+H+B+C-μ3*(N-S-E) +#lets test with integer order +_, x = FDEsolver(F, tSpan, x0, ones(8), par,h=.01,nc=4) +S=x[:,1]; E=x[:,2]; I1=x[:,3]; R=x[:,4]; +L=x[:,5]; H=x[:,6]; B=x[:,7]; C=x[:,8]; +N1=sum(x,dims=2) +μ3=mean([μ1 .+ μ2 .* N1, α1 .- α2 .*N1]) +Appx1=@.I1+R+L+H+B+C-μ3*(N1-S-E) + +# Dataset +Data=CSV.read("datoswho.csv", DataFrame) +data=(Matrix(Float64.(Data))) +Data2=LinearInterpolation(data[:,1], data[:,2]) + +#plotting +plt=plot(t,x,xlabel="time",lw = 2, label=["S" "E" "I" "R" "L" "H" "B" "C"], + thickness_scaling =1.2) +# savefig(plt,"plt.png") + +Err=rmsd(app[1:10:end,:], Data2(4:.1:T); normalize=:true) # Normalized root-mean-square error +Err1=rmsd(Appx1[1:10:end,:], Data2(4:.1:T); normalize=:true) # Normalized root-mean-square error +scatter(data[:,1],data[:,2],label="realdata",c="black") +plot!(t,app, lw=3,label="model1") +plot!(t,Appx1,lw=3,label="model2",linestyle=:dash,xlabel="time",legendposition=:right) +#plot population +plot(t,N,lw=3) +plot!(t,N1,lw=3,linestyle=:dash) + + +# Norm plot +Mf8=@.abs(Appx[1:10:end,:] - Data2(4:.1:T)) +M1=@.abs(Appx1[1:10:end,:] - Data2(4:.1:T)) + +plot(t[1:10:end,:],Mf8) +plot!(t[1:10:end,:],M1) diff --git a/Turring-ODE--ConstantN-Data.jl b/Turring-ODE--ConstantN-Data.jl new file mode 100644 index 0000000..8997364 --- /dev/null +++ b/Turring-ODE--ConstantN-Data.jl @@ -0,0 +1,133 @@ +# this code is for for optimizing μ and βr + +using Plots,SpecialFunctions, StatsBase, Random +using DifferentialEquations, Turing, LinearAlgebra +# Load StatsPlots for visualizations and diagnostics. +using StatsPlots +#initial conditons and parameters +using CSV +using DataFrames +using FdeSolver, Plots,SpecialFunctions, Optim, StatsBase, Random +using Interpolations +# Dataset +Data=CSV.read("datoswho.csv", DataFrame) + +data=(Matrix(Float64.(Data))) +#initial conditons and parameters + + +x0=[18000,0,15,0,0,0,0,0]# initial conditons S0,E0,I0,R0,L0,H0,B0,C0 +N=sum(x0) + +σ=1/11.4 # Per capita rate at which exposed individuals become infectious. +γ1=0.1 # Per capita rate of progression of individuals from the infectious class to the asymptomatic class. +γ2=1/5 # Per capita rate of progression of individuals from the hospitalized class to the asymptomatic class. +γ3=1/30 # Per capita recovery rate of individuals from the asymptomatic class to the complete recovered class. +ϵ=1/9.6 # Fatality rate. +δ1=1/2 # Per capita rate of progression of individuals from the dead class to the buried class. +δ2=1/4.6 # Per capita rate of progression of individuals from the hospitalized class to the buried class. +τ=1/5 # Per capita rate of progression of individuals from the infectious class to the hospitalized class. +βi=0.14 # Contact rate of infective individuals and susceptible. +βd=0.4 # Contact rate of infective individuals and dead. +βh=0.29 # Contact rate of infective individuals and hospitalized. +βr=0.305 # Contact rate of infective individuals and asymptomatic. +μ=14e-3 +ξ=14e-3 # Incineration rate + +par1=[σ,γ1,γ2,γ3,ϵ,δ1,δ2,τ,βi,βd,βh,ξ,N,βr,μ] # for model with constant N +par=[σ,γ1,γ2,γ3,ϵ,δ1,δ2,τ,βi,βd,βh,ξ,N] # for fitting + +tSpan=(4,438) + +#Define the equation + +function F1(dx, x, par, t) # model with constant N + + σ,γ1,γ2,γ3,ϵ,δ1,δ2,τ,βi,βd,βh,ξ,N,βr,μ=par + S, E, I, R, L, H, B, C=x + + dx[1]=μ*N - βi/N*S*I - βh/N*S*H - βd/N*S*L - βr/N*S*R - μ*S + dx[2]=βi/N*S*I + βh/N*S*H + βd/N*S*L + βr/N*S*R - σ*E - μ*E + dx[3]=σ*E - (γ1 + ϵ + τ + μ)*I + dx[4]=γ1*I + γ2*H - (γ3 + μ)*R + dx[5]=ϵ*I - (δ1+ξ)*L + dx[6]=τ*I - (γ2 + δ2 + μ)*H + dx[7]=δ1*L + δ2*H - ξ*B + dx[8]=γ3*R - μ*C + + return nothing + +end + +## optimazation of μ2 and α2 for integer order model +prob1 = ODEProblem(F1, x0, tSpan, par1) + +sol = solve(prob1, alg_hints=[:stiff]; saveat=0.1) +RealData=LinearInterpolation(data[:,1], data[:,2]) + + +@model function fitlv(data, prob2) + # Prior distributions. + σ ~ InverseGamma(2, 3) + βr ~ truncated(Normal(0.1, .9); lower=0.1, upper=.9) + μ ~ truncated(Normal(0, .1); lower=0, upper=.1) + + # Simulate model. + p=vcat(par,βr,μ) + xx = solve(prob2, alg_hints=[:stiff],reltol=1e-8,abstol=1e-8; p=p, saveat=1) + x=hcat(xx.u...) + S=x[1,:]; E=x[2,:]; I1=x[3,:]; R=x[4,:]; + L=x[5,:]; H=x[6,:]; B=x[7,:]; C=x[8,:]; + Appx=@.I1+R+L+H+B+C-μ*(N-S-E) + # Observations. + for i in 1:length(S) + data[i] ~ Normal(Appx[i], σ^2) + end + + return nothing +end + +model = fitlv(RealData(4:438), prob1) + +# Sample 3 independent chains with forward-mode automatic differentiation (the default). +chain = sample(model, NUTS(0.65), MCMCSerial(), 2000, 3; progress=false) + + +#plotting +pltChainConstN=plot(chain) +savefig(pltChainConstN,"pltChainConstN.svg") + +# # Save a chain. +# write("chain-file2.jls", chain) +# +# # Read a chain. +# chn2 = read("chain-file2.jls", Chains) +# julia> mean(chain[:βr]) +βr=0.7265002737432911 +# julia> mean(chain[:μ]) +μ=0.06918229886616623 + +plot(; legend=false) +posterior_samples = sample(chain[[:βr, :μ]], 2000; replace=false) +for p1 in eachrow(Array(posterior_samples)) + p=vcat(par,p1[1],p1[2]) + sol_p = solve(prob1, Tsit5(); p=p, saveat=0.1) + x=reduce(vcat,sol_p.u') + S=x[:,1]; E=x[:,2]; I1=x[:,3]; R=x[:,4]; + L=x[:,5]; H=x[:,6]; B=x[:,7]; C=x[:,8]; + Appx=@.I1+R+L+H+B+C-μ*(N-S-E) + plot!(sol_p.t,Appx; alpha=0.1, color="#BBBBBB") +end +# Plot simulation and noisy observations. +x=reduce(vcat,sol.u') +S=x[:,1]; E=x[:,2]; I1=x[:,3]; R=x[:,4]; +L=x[:,5]; H=x[:,6]; B=x[:,7]; C=x[:,8]; +Ex=@.I1+R+L+H+B+C-μ*(N-S-E) +pltConstN=scatter!(data[:,1],data[:,2]; color=[1 2], linewidth=1) + +savefig(pltConstN,"pltConstN.svg") + + +# \optimzed values' +βr=0.7265002737432911 +μ=0.06918229886616623 diff --git a/Turring-ODE-Data.jl b/Turring-ODE-Data.jl new file mode 100644 index 0000000..9fcedac --- /dev/null +++ b/Turring-ODE-Data.jl @@ -0,0 +1,266 @@ +# this code is for finding initial guess for the acceptable interval of α2 and μ2 + +using Plots,SpecialFunctions, StatsBase, Random +using DifferentialEquations, Turing, LinearAlgebra +# Load StatsPlots for visualizations and diagnostics. +using StatsPlots +#initial conditons and parameters +using CSV +using DataFrames +using FdeSolver, Plots,SpecialFunctions, Optim, StatsBase, Random +using Interpolations +# Dataset +Data=CSV.read("datoswho.csv", DataFrame) + +data=(Matrix(Float64.(Data))) +#initial conditons and parameters + + +x0=[18000,0,15,0,0,0,0,0]# initial conditons S0,E0,I0,R0,L0,H0,B0,C0 +N=sum(x0) + +α1=35.37e-3 # Density independent part of the birth rate for individuals. +α2=.1*α1 # Density dependent part of the birth rate for individuals. +σ=1/11.4 # Per capita rate at which exposed individuals become infectious. +γ1=0.1 # Per capita rate of progression of individuals from the infectious class to the asymptomatic class. +γ2=1/5 # Per capita rate of progression of individuals from the hospitalized class to the asymptomatic class. +γ3=1/30 # Per capita recovery rate of individuals from the asymptomatic class to the complete recovered class. +ϵ=1/9.6 # Fatality rate. +δ1=1/2 # Per capita rate of progression of individuals from the dead class to the buried class. +δ2=1/4.6 # Per capita rate of progression of individuals from the hospitalized class to the buried class. +τ=1/5 # Per capita rate of progression of individuals from the infectious class to the hospitalized class. +βi=0.14 # Contact rate of infective individuals and susceptible. +βd=0.4 # Contact rate of infective individuals and dead. +βh=0.29 # Contact rate of infective individuals and hospitalized. +βr=0.305 # Contact rate of infective individuals and asymptomatic. +# βr=0.185 +μ=14e-3 +μ1=10.17e-3 # Density independent part of the death rate for individuals. +μ2=.1*μ1 # Density dependent part of the death rate for individuals. +ξ=14e-3 # Incineration rate + +par1=[α1,σ,γ1,γ2,γ3,ϵ,δ1,δ2,τ,βi,βd,βh,βr,μ1,ξ,N] # for model with constant N +par2=[α1,σ,γ1,γ2,γ3,ϵ,δ1,δ2,τ,βi,βd,βh,βr,μ1,ξ,α2,μ2] # for model with variable N +# par3=[α1,σ,γ1,γ2,ϵ,δ1,δ2,τ,βi,βd,βh,μ1,ξ] # I need it for fitting +par3=[α1,σ,γ1,γ2,γ3,ϵ,δ1,δ2,τ,βi,βd,βh,μ1,ξ] # I need it for fitting + +tSpan=(4,438) + +#Define the equation + +function F1(dx, x, par, t) # model with constant N + + α1,σ,γ1,γ2,ϵ,δ1,δ2,τ,βi,βd,βh,βr,μ1,ξ,N=par + S, E, I, R, L, H, B, C=x + + dx[1]=μ*N - βi/N*S*I - βh/N*S*H - βd/N*S*L - βr/N*S*R - μ*S + dx[2]=βi/N*S*I + βh/N*S*H + βd/N*S*L + βr/N*S*R - σ*E - μ*E + dx[3]=σ*E - (γ1 + ϵ + τ + μ)*I + dx[4]=γ1*I + γ2*H - (γ3 + μ)*R + dx[5]=ϵ*I - (δ1+ξ)*L + dx[6]=τ*I - (γ2 + δ2 + μ)*H + dx[7]=δ1*L + δ2*H - ξ*B + dx[8]=γ3*R - μ*C + + return nothing + +end + +function F2(dx, x, par, t) # model with variable N + + α1,σ,γ1,γ2,γ3,ϵ,δ1,δ2,τ,βi,βd,βh,μ1,ξ,βr,α2,μ2=par + S, E, I, R, L, H, B, C=x + N=sum(x) + + dx[1]=(α1-α2*N)*N - βi/N*S*I - βh/N*S*H - βd/N*S*L - βr/N*S*R - (μ1+μ2*N)*S + dx[2]=βi/N*S*I + βh/N*S*H + βd/N*S*L + βr/N*S*R - σ*E - (μ1+μ2*N)*E + dx[3]=σ*E - (γ1 + ϵ + τ + μ1 + μ2*N)*I + dx[4]=γ1*I + γ2*H - (γ3 + μ1 + μ2*N)*R + dx[5]=ϵ*I - (δ1+ξ)*L + dx[6]=τ*I - (γ2 + δ2 + μ1 + μ2*N)*H + dx[7]=δ1*L + δ2*H - ξ*B + dx[8]=γ3*R - (μ1 + μ2*N)*C + + return nothing + +end + +## optimazation of μ2 and α2 for integer order model +prob1 = ODEProblem(F1, x0, tSpan, par1) +prob2 = ODEProblem(F2, x0, tSpan, par2) + +sol = solve(prob1, alg_hints=[:stiff]; saveat=0.1) +RealData=LinearInterpolation(data[:,1], data[:,2]) + + +@model function fitlv(data, prob2) + # Prior distributions. + σ ~ InverseGamma(2, 3) + α2 ~ truncated(Normal(0, .0001); lower=0, upper=.0001) + μ2 ~ truncated(Normal(0, .0001); lower=0, upper=.0001) + βr ~ truncated(Normal(0.1, .9); lower=0.1, upper=.9) + # γ3 ~ truncated(Normal(0.001, .9); lower=0.001, upper=.9) + + # Simulate model. + # p=vcat(par3,γ3,βr,α2,μ2) + p=vcat(par3,βr,α2,μ2) + xx = solve(prob2, alg_hints=[:stiff],reltol=1e-8,abstol=1e-8; p=p, saveat=1) + x=hcat(xx.u...) + S=x[1,:]; E=x[2,:]; I1=x[3,:]; R=x[4,:]; + L=x[5,:]; H=x[6,:]; B=x[7,:]; C=x[8,:]; + N=sum(x,dims=1)' + μ3=mean([μ1 .+ μ2 .* N, α1 .- α2 .*N]) + Appx=@.I1+R+L+H+B+C-μ3[:,1]*(N[:,1]-S-E) + # Observations. + for i in 1:length(S) + data[i] ~ Normal(Appx[i], σ^2) + end + + return nothing +end + +model = fitlv(RealData(4:438), prob2) + +# Sample 3 independent chains with forward-mode automatic differentiation (the default). +chain = sample(model, NUTS(0.65), MCMCSerial(), 2000, 3; progress=false) + + +#plotting +pltChain4=plot(chain) +savefig(pltChain4,"pltChain4.svg") + +# Save a chain. +write("chain-file2.jls", chain) + +# Read a chain. +chn2 = read("chain-file2.jls", Chains) + +plot(; legend=false) +posterior_samples = sample(chain[[:α2, :μ2]], 2000; replace=false) +for p1 in eachrow(Array(posterior_samples)) + p=vcat(par3,p1[1],p1[2]) + sol_p = solve(prob2, Tsit5(); p=p, saveat=0.1) + x=reduce(vcat,sol_p.u') + S=x[:,1]; E=x[:,2]; I1=x[:,3]; R=x[:,4]; + L=x[:,5]; H=x[:,6]; B=x[:,7]; C=x[:,8]; + N1=sum(x,dims=2) + μ3=mean([μ1 .+ p1[2] .* N1, α1 .- p1[1] .*N1]) + Appx=@.I1+R+L+H+B+C-μ3*(N1-S-E) + plot!(sol_p.t,Appx; alpha=0.1, color="#BBBBBB") +end +# Plot simulation and noisy observations. +x=reduce(vcat,sol.u') +S=x[:,1]; E=x[:,2]; I1=x[:,3]; R=x[:,4]; +L=x[:,5]; H=x[:,6]; B=x[:,7]; C=x[:,8]; +Ex=@.I1+R+L+H+B+C-μ*(N-S-E) +plt305=scatter!(data[:,1],data[:,2]; color=[1 2], linewidth=1) + +savefig(plt305,"pltFit305.svg") + +#optimized values + +# mean(chain[:α2]) +α2=8.334570273681649e-8 + +# mean(chain[:μ2]) +μ2=4.859965522407347e-7 + +# mean(chain[:βr]) +βr=0.2174851498937417 + +#let's plot the results of two models +x=reduce(vcat,sol.u') +S=x[:,1]; E=x[:,2]; I1=x[:,3]; R=x[:,4]; +L=x[:,5]; H=x[:,6]; B=x[:,7]; C=x[:,8]; +Ex=@.I1+R+L+H+B+C-μ*(N-S-E) +plot(sol.t, Ex; lw=1, label="Model1") + +# p=vcat(par3,γ3,βr,α2,μ2) +p=vcat(par3,βr,α2,μ2) +# p=vcat(par3,p1[1],p1[2]) +sol_p = solve(prob2, alg_hints=[:stiff]; p=p, saveat=0.1) +x=reduce(vcat,sol_p.u') +S=x[:,1]; E=x[:,2]; I1=x[:,3]; R=x[:,4]; +L=x[:,5]; H=x[:,6]; B=x[:,7]; C=x[:,8]; +N1=sum(x,dims=2) +μ3=mean([μ1 .+ μ2 .* N1, α1 .- α2 .*N1]) +Appx=@.I1+R+L+H+B+C-μ3*(N1-S-E) +plot!(sol_p.t,Appx, label="Model2") +scatter!(data[:,1],data[:,2],legendposition=:right) + +plot(sol) +plot(N1) + +using CSV +using DataFrames +using Interpolations +# Dataset +Data=CSV.read("plot-data.csv", DataFrame) + +data=(Matrix(Data)) +Data2=LinearInterpolation(data[:,1], data[:,2]) +plt=scatter!(Data2(1:85),label="Real data",xlabel="days", ylabel="Cumulative confirmed cases", + legendposition=:bottomright) + +savefig(plt,"plt.png") +pltN=plot(sol_p.t,N1,ylabel="Population N",xlabel="days", + legend=:false) + +savefig(pltN,"pltN.png") +#################for βr=305 + +function F11(dx, x, par, t) # model with constant N + + α1,σ,γ1,γ2,γ3,ϵ,δ1,δ2,τ,βi,βd,βh,βr,μ1,ξ,N=par + # βr=0.185 # why is it different value in the papers + S, E, I, R, L, H, B, C=x + + dx[1]=μ*N - βi/N*S*I - βh/N*S*H - βd/N*S*L - βr/N*S*R - μ*S + dx[2]=βi/N*S*I + βh/N*S*H + βd/N*S*L + βr/N*S*R - σ*E - μ*E + dx[3]=σ*E - (γ1 + ϵ + τ + μ)*I + dx[4]=γ1*I + γ2*H - (γ3 + μ)*R + dx[5]=ϵ*I - (δ1+ξ)*L + dx[6]=τ*I - (γ2 + δ2 + μ)*H + dx[7]=δ1*L + δ2*H - ξ*B + dx[8]=γ3*R - μ*C + + return nothing + +end +prob1 = ODEProblem(F11, x0, tSpan, par1) + +sol = solve(prob1, alg_hints=[:stiff]; saveat=0.1) +x=reduce(vcat,sol.u') +S=x[:,1]; E=x[:,2]; I1=x[:,3]; R=x[:,4]; +L=x[:,5]; H=x[:,6]; B=x[:,7]; C=x[:,8]; +Ex=@.I1+R+L+H+B+C-μ*(N-S-E) +plot(sol.t, Ex; label="Model1", lw=3) + +p1=[1.1862350479039923e-6,2.1259887426666564e-7] +p=vcat(par3,p1[1],p1[2]) +sol_p = solve(prob2, alg_hints=[:stiff]; p=p, saveat=0.1) +x=reduce(vcat,sol_p.u') +S=x[:,1]; E=x[:,2]; I1=x[:,3]; R=x[:,4]; +L=x[:,5]; H=x[:,6]; B=x[:,7]; C=x[:,8]; +N1=sum(x,dims=2) +μ3=mean([μ1 .+ p1[2] .* N1, α1 .- p1[1] .*N1]) +Appx=@.I1+R+L+H+B+C-μ3*(N1-S-E) +plot!(sol_p.t,Appx, label="Model2",lw=3, linestyle=:dash) + + +using CSV +using DataFrames +using Interpolations +# Dataset +Data=CSV.read("plot-data.csv", DataFrame) + +data=(Matrix(Data)) +Data2=LinearInterpolation(data[:,1], data[:,2]) +plt=scatter!(Data2(1:85),label="Real data",xlabel="days", ylabel="Cumulative confirmed cases", + legendposition=:left) + +savefig(plt,"plt2.png") +pltN2=plot(sol_p.t,N1,ylabel="Population N",xlabel="days", + legend=:false) + +savefig(pltN2,"pltN2.png") diff --git a/pltFtry.svg b/pltFtry.svg new file mode 100644 index 0000000..8653358 --- /dev/null +++ b/pltFtry.svg @@ -0,0 +1,3153 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/pltheat.svg b/pltheat.svg new file mode 100644 index 0000000..d578829 --- /dev/null +++ b/pltheat.svg @@ -0,0 +1,385 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/pltheat1.svg b/pltheat1.svg new file mode 100644 index 0000000..adbd3aa --- /dev/null +++ b/pltheat1.svg @@ -0,0 +1,381 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +