Skip to content

Commit

Permalink
polish the plots
Browse files Browse the repository at this point in the history
  • Loading branch information
moeinkh88 committed Dec 6, 2022
1 parent 774a357 commit fcd4500
Show file tree
Hide file tree
Showing 6 changed files with 94,035 additions and 41,084 deletions.
74,312 changes: 37,174 additions & 37,138 deletions Allplots.svg
View file
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
60 changes: 33 additions & 27 deletions Allplots_Models.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3,7 +3,7 @@ using CSV
using DataFrames
using FdeSolver, Plots,SpecialFunctions, Optim, StatsBase, Random
using Interpolations, LinearAlgebra

using LaTeXStrings

# Dataset
Data=CSV.read("datoswho.csv", DataFrame)
Expand Down Expand Up @@ -85,7 +85,7 @@ T1=250
tSpan2=[T1,T2]

#solve model 1: ODE
t, x= FDEsolver(F, [4,T2], x0, ones(8), par, nc=4,h=.01)
_, 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)
Expand All @@ -106,7 +106,7 @@ 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)
_, 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)
Expand Down Expand Up @@ -157,22 +157,25 @@ t, x= FDEsolver(F, [4,450], x0, αf, par, nc=4,h=.01)
μ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)
_, x1_2= 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,:])
_, x2_2= FDEsolver(F, [T1,450], x02, α, par, nc=4,h=.01)
xf=vcat(x1_2,x2_2[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)
AppxF2=@.If2+Rf2+Lf2+Hf2+Bf2+Cf2-μ3*(Nf2-Sf2-Ef2)''

##plot dynamics
gr()
# pyplot()
plot_font="Computer Modern"
# plotfonts = Plots.font(16, "Helvetica")
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,AppxC,ylabel=L"\textrm{Cumulative\;\;cases} \; (\mathrm{10^3})",lw=2, label="Model1, RMSD=2085", c="darkorange3", linestyle=:dashdot)
plot!(t,Appx1,lw=2, label="Model2, RMSD=276.5", c="deepskyblue1",guidefont=font(12,plot_font),
yticks = ([0,5000,10000,15000],string.([0,5,10,15])),fontfamily= plot_font)
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))
Expand All @@ -181,23 +184,24 @@ scatter(data[:,1],data[:,2], label="Real data", c="khaki3",markerstrokewidth=0)
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))
plotN=plot!(t,Nf2,label="Model4",linestyle=:dot, lw=2, xlabel="Days",guidefont=font(12,plot_font),
ylabel=L"\textrm{Population} \;\; (\mathrm{10^4})", c="blue1", fontfamily="computer modern" ,
yticks = ([20000,30000,40000,50000],string.([2,3,4,5])),
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
# 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)
Expand All @@ -214,7 +218,7 @@ end
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
# 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,:]
Expand All @@ -234,18 +238,18 @@ end
RMSD1= log10(Err1) #error for integer order

# plot heatmap model 4
CoLor = cgrad([:yellow2,:darkgoldenrod1,:mediumpurple], [.45,.4], categorical = true) # define color
CoLor = cgrad([:darkblue,:blue, :lightblue,:orange], [.4,.5]) # 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",
xlabel="Individual classes", ylabel="Order of derivative",fontfamily= plot_font,
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)
pltheat=plot!(Shape(0 .+ [1,2,2,1], 0 .+ [.595,.595,.995,.995]), fillcolor=:false, linecolor=:khaki1,linestyle=:dot, legend=:false, w=3)


# 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",
xlabel="Individual classes", ylabel="Order of derivative",fontfamily= plot_font,
title = "(d) Model3 ", titleloc = :left, titlefont = font(10), colorbar=:false)


Expand All @@ -258,35 +262,37 @@ tt, x= FDEsolver(F, [4,438], x0, ones(8), par, nc=2,h=.01)
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)
Appx11=@.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
Nn=30
err=zeros(Nn)
for i=1:Nn
α=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,:])
ttt=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")
plot!(ttt[25000:100:end],Appx[25000:100:end,:],lw=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,
pltFtry=plot!(tt[25000:50:end],Appx11[25000:50:end,:],xlabel="Days",
lw=2, c="deepskyblue1",label="Integer model",legend=:false,
ylabel=L"\textrm{Cumulative\;\;cases} \; (\mathrm{10^3})",fontfamily= plot_font,
yticks = (1000 .*[11,12,13,14,15],string.([11,12,13,14,15])),
title = "(c)", titleloc = :left, titlefont = font(10))

# savefig(pltFtry,"pltFtry.svg")
Expand Down
Loading

0 comments on commit fcd4500

Please sign in to comment.