Hello,
I’m using the following SIRD eipdemiological model. This model takes into account 16 age classes. The last dD_tot compartment, an array of 16 components all equal is a trick to later calibrate with respect to cumulated deaths.
function SIRD!(du,u,p,t)
# Parameter set
β, λ_R, λ_D = p
δ₁, δ₂, δ₃, δ₄ = [0.003/100, 0.004/100, (0.015+0.030+0.064+0.213+0.718)/(5*100), (2.384+8.466+12.497+1.117)/(4*100)]
δ = vcat(repeat([δ₁],2),repeat([δ₂],2),repeat([δ₃],9),repeat([δ₄],16-2-2-9))
# Contact matrix
C = regional_all_contact_matrix
# State variables
S = @view u[16*0+1:16*1]
I = @view u[16*1+1:16*2]
R = @view u[16*2+1:16*3]
D = @view u[16*3+1:16*4]
D_tot = @view u[16*4+1:16*5]
# Differentials
dS = @view du[16*0+1:16*1]
dI = @view du[16*1+1:16*2]
dR = @view du[16*2+1:16*3]
dD = @view du[16*3+1:16*4]
dD_tot = @view du[16*4+1:16*5]
# Force of infection
Λ = β*[sum([C[i,j]*I[j]/N[j] for j in 1:size(C)[1]]) for i in 1:size(C)[2]]
# System of equations
@. dS = -Λ*S
@. dI = Λ*S - ((1-δ)*λ_R + δ*λ_D)*I
@. dR = λ_R*(1-δ)*I
@. dD = λ_D*δ*I
@. dD_tot = dD[1]+dD[2]+dD[3]+dD[4]+dD[5]+dD[6]+dD[7]+dD[8]+dD[9]+dD[10]+dD[11]+dD[12]+dD[13]+dD[14]+dD[15]+dD[16]
end;
And I would like to calibrate its parameters β, λ_R, λ_D in the range [0.0, 1.0] using real data about cumulative deaths.
I used the following code (piedmont_cumulative_deaths[1:final_time] is an array of integers, each of representing cumulative daily deaths):
initial = [0.049001 , 0.541157 , 0.825232]
lower = [0.0, 0,0, 0.0]
upper = [1.0,1.0,1.0]
problem = ODEProblem(SIRD!, ℬ, 𝒯, [ β ,λ_R , λ_D ])
cost_function = build_loss_objective(problem,Tsit5(),L2Loss(t,piedmont_cumulative_deaths[1:final_time]), maxiters=10000,verbose=true, save_idxs= [16*4+1])
result_optim = Optim.optimize(cost_function,lower,upper, initial, Optim.Fminbox(BFGS()) ,Optim.Options(show_trace = true , iterations =50 ))
But I get:
DimensionMismatch("array could not be broadcast to match destination")
Stacktrace:
[1] check_broadcast_shape at .\broadcast.jl:520 [inlined]
[2] check_broadcast_axes at C:\Users\claud\.julia\packages\DiffEqBase\V7P18\src\diffeqfastbc.jl:26 [inlined]
[3] check_broadcast_axes at .\broadcast.jl:526 [inlined]
[4] check_broadcast_axes at .\broadcast.jl:527 [inlined]
[5] instantiate at .\broadcast.jl:269 [inlined]
[6] materialize! at .\broadcast.jl:848 [inlined]
[7] materialize! at .\broadcast.jl:845 [inlined]
[8] value_gradient!!(::Optim.BarrierWrapper{OnceDifferentiable{Float64,Array{Float64,2},Array{Float64,2}},Optim.BoxBarrier{Array{Float64,1},Array{Float64,1}},Float64,Float64,Array{Float64,2}}, ::Array{Float64,2}) at C:\Users\claud\.julia\packages\Optim\auGGa\src\multivariate\solvers\constrained\fminbox.jl:65
[9] initial_state(::BFGS{LineSearches.InitialStatic{Float64},LineSearches.HagerZhang{Float64,Base.RefValue{Bool}},Nothing,Nothing,Flat}, ::Optim.Options{Float64,Nothing}, ::Optim.BarrierWrapper{OnceDifferentiable{Float64,Array{Float64,2},Array{Float64,2}},Optim.BoxBarrier{Array{Float64,1},Array{Float64,1}},Float64,Float64,Array{Float64,2}}, ::Array{Float64,2}) at C:\Users\claud\.julia\packages\Optim\auGGa\src\multivariate\solvers\first_order\bfgs.jl:85
[10] optimize(::OnceDifferentiable{Float64,Array{Float64,2},Array{Float64,2}}, ::Array{Float64,1}, ::Array{Float64,1}, ::Array{Float64,2}, ::Fminbox{BFGS{LineSearches.InitialStatic{Float64},LineSearches.HagerZhang{Float64,Base.RefValue{Bool}},Nothing,Nothing,Flat},Float64,Optim.var"#45#47"}, ::Optim.Options{Float64,Nothing}) at C:\Users\claud\.julia\packages\Optim\auGGa\src\multivariate\solvers\constrained\fminbox.jl:307
[11] #optimize#61 at C:\Users\claud\.julia\packages\Optim\auGGa\src\multivariate\solvers\constrained\fminbox.jl:254 [inlined]
[12] optimize(::Function, ::Array{Float64,1}, ::Array{Float64,1}, ::Array{Float64,2}, ::Fminbox{BFGS{LineSearches.InitialStatic{Float64},LineSearches.HagerZhang{Float64,Base.RefValue{Bool}},Nothing,Nothing,Flat},Float64,Optim.var"#45#47"}, ::Optim.Options{Float64,Nothing}) at C:\Users\claud\.julia\packages\Optim\auGGa\src\multivariate\solvers\constrained\fminbox.jl:253
[13] top-level scope at In[98]:7
[14] include_string(::Function, ::Module, ::String, ::String) at .\loading.jl:1091
[15] execute_code(::String, ::String) at C:\Users\claud\.julia\packages\IJulia\a1SNk\src\execute_request.jl:27
[16] execute_request(::ZMQ.Socket, ::IJulia.Msg) at C:\Users\claud\.julia\packages\IJulia\a1SNk\src\execute_request.jl:86
[17] #invokelatest#1 at .\essentials.jl:710 [inlined]
[18] invokelatest at .\essentials.jl:709 [inlined]
[19] eventloop(::ZMQ.Socket) at C:\Users\claud\.julia\packages\IJulia\a1SNk\src\eventloop.jl:8
[20] (::IJulia.var"#15#18")() at .\task.jl:356
Additionally, how may one calibrate with respect to multiple compartments ( say, for example, I also have another array of measurements about cumulative daily recovered, called piedmont_cumulative_recovered)?
Thank you very much
.