Using DifferentialEquations.jl Fminbox with multiple parameters

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

2 Likes

Most of the variables do not seem to be defined in your code?

You’re right. Here’s the code.

Packages:

# Basic Utilities
using StaticArrays, RecursiveArrayTools, ParameterizedFunctions
# Numerical Computation
using DifferentialEquations, LinearAlgebra, DiffEqBayes, OrdinaryDiffEq, DiffEqParamEstim, Optim
# Statistics 
using Random, Distributions
# Data Wrangling 
using Queryverse, DataFrames, HTTP, CSV

Parameters:


final_time = 20 # day up to which calibration will be performed

# Model parameters

β = 0.05 # infection rate
λ_R = 0.55 # inverse of transition time from  infected to recovered
λ_D = 0.83 # inverse of transition time from  infected to dead
i₀ = 0.075 # fraction of initial infected people in every age class
𝒫 = vcat([β, δ, λ_R, λ_D]...)

# regional contact matrix and regional population
## regional contact matrix
all_contact_path = joinpath(@__DIR__, "All.csv");
regional_all_contact_matrix = SMatrix{16,16}(convert(Matrix,select!(DataFrame(load(all_contact_path)), Not(1)))); # 16x16 contact matrix

## regional population stratified by age
regional_population_path      = joinpath(@__DIR__, "AgeStratifiedRegionalPopulation.csv")
N= DataFrame(load(regional_population_path))[1] # array of 16 elements, each of which representing the absolute amount of population in the corresponding age class.

# Initial conditions 
I₀ = repeat([i₀],16)
S₀ = N.-I₀
R₀ = [0.0 for n in 1:length(N)]
D₀ = [0.0 for n in 1:length(N)]
D_tot₀ = [0.0 for n in 1:length(N)]
ℬ = vcat([S₀, I₀, R₀, D₀, D_tot₀]...) 

# Time 
𝒯 = (1.0,final_time); 

# data used for calibration

## get regional level epidemic data
regional_epi_data = DataFrame(CSV.File(HTTP.get("https://raw.githubusercontent.com/pcm-dpc/COVID-19/master/dati-regioni/dpc-covid19-ita-regioni.csv").body))
## extract data of the region of interest (Piedmont)
piedmont_epi_data = filter(row -> row[:codice_regione] == 1, regional_epi_data)
## extract calibration data (dauly cumulative deaths)
piedmont_cumulative_deaths = piedmont_epi_data.deceduti

We use a SIRD model. This model takes into account the 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
    # initializa this parameter (death probability stratified by age, taken from literature)
    δ₁, δ₂, δ₃, δ₄ = [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;

Then i try to do bounded calibration:

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 the error:

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

I don’t understamd the reasons behind this error.

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

1 Like

I don’t have the CSV to run this… can you try simplifying it? If you simplify it, you’ll know exactly what line causes it. That will either solve the problem or it would help you ask the right question.

Thanks for your time. Sorry for not posting a self-contained question from the beginning.
I simplified those csv and hard-coded them below. So now there are two differences:

  1. Now there are 4 age classes insteas of 16 ( these is the result of the semplification)

  2. The error changed sligthly: it seems i cannot reproduce the same for some reason. Anyway the problematic lines are always the same.

Let me repeat the purpose of this question for clarity:I have a SIRD model which has three parameters that I would like to calibrate constraining the values they may assume using Fminbox. Unfortunately, I cannot run the code that should do it.

Packages:

# Basic Utilities
using StaticArrays, RecursiveArrayTools, ParameterizedFunctions
# Numerical Computation
using DifferentialEquations, LinearAlgebra, DiffEqBayes, OrdinaryDiffEq, DiffEqParamEstim, Optim
# Statistics 
using Random, Distributions
# Data Wrangling 
using Queryverse, DataFrames, HTTP, CSV

Parameters:


final_time = 20 # day up to which calibration will be performed

# Model parameters

β = 0.05 # infection rate
λ_R = 0.55 # inverse of transition time from  infected to recovered
λ_D = 0.83 # inverse of transition time from  infected to dead
i₀ = 0.075 # fraction of initial infected people in every age class
𝒫 = vcat([β, δ, λ_R, λ_D]...)

# regional contact matrix and regional population
## regional contact matrix

regional_all_contact_matrix = [3.45536   0.485314  0.506389  0.123002 ; 0.597721  2.11738   0.911374  0.323385 ; 0.906231  1.35041   1.60756   0.67411 ; 0.237902  0.432631  0.726488  0.979258] # 4x4 contact matrix

## regional population stratified by age

N= [723208 , 874150, 1330993, 1411928] # array of 4 elements, each of which representing the absolute amount of population in the corresponding age class.

# Initial conditions 
I₀ = repeat([i₀],4)
S₀ = N.-I₀
R₀ = [0.0 for n in 1:length(N)]
D₀ = [0.0 for n in 1:length(N)]
D_tot₀ = [0.0 for n in 1:length(N)]
ℬ = vcat([S₀, I₀, R₀, D₀, D_tot₀]...) 

# Time 
𝒯 = (1.0,final_time); 

# data used for calibration
piedmont_cumulative_deaths = [0,2,4,5,5,13,17,21,26,46,59,81,111,133,154,175,209,238,283,315,374]

We use a SIRD model. This model takes into account the 4 age classes. The last dD_tot compartment, an array of 4 components all equal is just a trick to later calibrate with respect to cumulated deaths.

function SIRD!(du,u,p,t)  
    # Parameters to be calibrated
    β, λ_R, λ_D = p

    # initializa this parameter (death probability stratified by age, taken from literature)
    δ₁, δ₂, δ₃, δ₄ = [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([δ₁],1),repeat([δ₂],1),repeat([δ₃],1),repeat([δ₄],4-1-1-1))

    
    C = regional_all_contact_matrix 
    
    # State variables
    S = @view u[4*0+1:4*1]
    I = @view u[4*1+1:4*2]
    R = @view u[4*2+1:4*3]
    D = @view u[4*3+1:4*4]
    D_tot = @view u[4*4+1:4*5]

    # Differentials
    dS = @view du[4*0+1:4*1]
    dI = @view du[4*1+1:4*2]
    dR = @view du[4*2+1:4*3]
    dD = @view du[4*3+1:4*4]
    dD_tot = @view du[4*4+1:4*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]

end;

Then I try to do bounded calibration:

initial = [0.049001 , 0.541157 , 0.825232] #initial parameters values
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), maxiters=10000,verbose=true, save_idxs= [4*4+1])
result_optim = Optim.optimize(cost_function,lower,upper, initial, Optim.Fminbox(BFGS()) ,Optim.Options(show_trace = true , iterations =50 ))

But now I get the error:

BoundsError: attempt to access 3-element Array{Float64,1} at index [4]

Stacktrace:
 [1] getindex at .\array.jl:809 [inlined]
 [2] optimize(::OnceDifferentiable{Float64,Array{Float64,1},Array{Float64,1}}, ::Array{Float64,1}, ::Array{Float64,1}, ::Array{Float64,1}, ::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:282
 [3] #optimize#61 at C:\Users\claud\.julia\packages\Optim\auGGa\src\multivariate\solvers\constrained\fminbox.jl:254 [inlined]
 [4] optimize(::Function, ::Array{Float64,1}, ::Array{Float64,1}, ::Array{Float64,1}, ::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
 [5] top-level scope at In[56]:6
 [6] include_string(::Function, ::Module, ::String, ::String) at .\loading.jl:1091
 [7] execute_code(::String, ::String) at C:\Users\claud\.julia\packages\IJulia\a1SNk\src\execute_request.jl:27
 [8] execute_request(::ZMQ.Socket, ::IJulia.Msg) at C:\Users\claud\.julia\packages\IJulia\a1SNk\src\execute_request.jl:86
 [9] #invokelatest#1 at .\essentials.jl:710 [inlined]
 [10] invokelatest at .\essentials.jl:709 [inlined]
 [11] eventloop(::ZMQ.Socket) at C:\Users\claud\.julia\packages\IJulia\a1SNk\src\eventloop.jl:8
 [12] (::IJulia.var"#15#18")() at .\task.jl:356

Expected behaviour: It calibrates the three parameters [ β ,λ_R , λ_D ], specified when creating the ODEProblem in the last code cell.

Observation: The array that causes the error seems to be initial, which contains the initial estimations of the three parameters i would like to calibrate. Adding elements to it makes no sense to me.

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 ) using the methods from Optim?

Thank you very much, please don’t hesitate to ask for other information or code in case I missed it.

2 Likes

I still can’t run the code because t and piedmont_cumulative_recovered were not included. Can you just give it as sample or fake data?

You probably want to use DiffEqFlux.jl instead if you have a complicated loss function in mind.

I don’t have real piedmont_cumulative_recovered, anyway I created it as fake data below (also inserted t)

# Basic Utilities
using StaticArrays, RecursiveArrayTools, ParameterizedFunctions
# Numerical Computation
using DifferentialEquations, LinearAlgebra, DiffEqBayes, OrdinaryDiffEq, DiffEqParamEstim, Optim
# Statistics 
using Random, Distributions
# Data Wrangling 
using Queryverse, DataFrames, HTTP, CSV


final_time = 20 # day up to which calibration will be performed

# Model parameters

β = 0.05 # infection rate
λ_R = 0.55 # inverse of transition time from  infected to recovered
λ_D = 0.83 # inverse of transition time from  infected to dead
i₀ = 0.075 # fraction of initial infected people in every age class
𝒫 = vcat([β, λ_R, λ_D]...)

# regional contact matrix and regional population
## regional contact matrix

regional_all_contact_matrix = [3.45536   0.485314  0.506389  0.123002 ; 0.597721  2.11738   0.911374  0.323385 ; 0.906231  1.35041   1.60756   0.67411 ; 0.237902  0.432631  0.726488  0.979258] # 4x4 contact matrix

## regional population stratified by age

N= [723208 , 874150, 1330993, 1411928] # array of 4 elements, each of which representing the absolute amount of population in the corresponding age class.

# Initial conditions 
I₀ = repeat([i₀],4)
S₀ = N.-I₀
R₀ = [0.0 for n in 1:length(N)]
D₀ = [0.0 for n in 1:length(N)]
D_tot₀ = [0.0 for n in 1:length(N)]
ℬ = vcat([S₀, I₀, R₀, D₀, D_tot₀]...) 

# Time 
𝒯 = (1.0,final_time); 

# data used for calibration
piedmont_cumulative_deaths = [0,2,4,5,5,13,17,21,26,46,59,81,111,133,154,175,209,238,283,315]
piedmont_cumulative_recovered= [R[i] +  0.01*randn(1)[1] for i in  1:length(R)]


function SIRD!(du,u,p,t)  
    # Parameters to be calibrated
    β, λ_R, λ_D = p

    # initializa this parameter (death probability stratified by age, taken from literature)
    δ₁, δ₂, δ₃, δ₄ = [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([δ₁],1),repeat([δ₂],1),repeat([δ₃],1),repeat([δ₄],4-1-1-1))

    
    C = regional_all_contact_matrix 
    
    # State variables
    S = @view u[4*0+1:4*1]
    I = @view u[4*1+1:4*2]
    R = @view u[4*2+1:4*3]
    D = @view u[4*3+1:4*4]
    D_tot = @view u[4*4+1:4*5]

    # Differentials
    dS = @view du[4*0+1:4*1]
    dI = @view du[4*1+1:4*2]
    dR = @view du[4*2+1:4*3]
    dD = @view du[4*3+1:4*4]
    dD_tot = @view du[4*4+1:4*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]

end;

# generate fake piedmont_cumulative_recovered
problem = ODEProblem(SIRD!, ℬ, 𝒯, 𝒫)
solution_all = @time solve(problem, saveat=1:final_time);
piedmont_cumulative_recovered = [sum(solution_all.u[t][4*2+1:4*3]) +  0.01*randn(1)[1] for t in 1:length(solution_all.t)];


t = collect(1.:final_time);  
initial = [0.049001 , 0.541157 , 0.825232] #initial parameters values
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), maxiters=10000,verbose=true, save_idxs= [4*4+1])
result_optim = Optim.optimize(cost_function,lower,upper, initial, Optim.Fminbox(BFGS()) ,Optim.Options(show_trace = true , iterations =50 ))

You had a typo here. It should be lower = [0.0, 0.0, 0.0]. The error message is admittedly bad but if you clicked on the stack trace it highlights:

    for i in eachindex(l)
        thisx = x[i]

so it loops over the indices of l and errors on indexing x, so the size of the lower bound must be larger than the upper bound, and then I looked at the lower bound code :slight_smile:.

@pkofod is there a place in Optim to put some global error messages to catch this?

1 Like