Hello,
I’m trying to optimize the code and finding the best solve algorithm for an epidemiological model. Anyway, I get an error that due to my inexperience I cannot solve.
So below is a simplified verison of the model:
using DifferentialEquations, Distributions, DiffEqBayes, AutoOptimize
# Model parameters
β = 0.01# infection rate
λ_R = 0.05 # 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
final_time = 20
𝒯 = (1.0,final_time);
function SIRD_ac!(du,u,p,t)
# Parameters to be calibrated
β, λ_R, λ_D = p
# initialize 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;
# create problem and check it works
problem = ODEProblem(SIRD_ac!, ℬ, 𝒯, 𝒫)
solution = solve(problem, saveat = 1:final_time);
Then this works:
problem = ODEProblem(SIRD_ac!, ℬ, 𝒯, 𝒫)
sys = modelingtoolkitize(problem)
But when I try to rebuild the ODEProblem I get:
fast_problem = ODEProblem(sys,ℬ, 𝒯, 𝒫 )
BoundsError: attempt to access Float64
at index [2]
Stacktrace:
[1] indexed_iterate(::Float64, ::Int64, ::Nothing) at .\tuple.jl:90
[2] (::ModelingToolkit.var"#208#210"{ODESystem})(::Float64) at .\none:0
[3] iterate at .\generator.jl:47 [inlined]
[4] collect at .\array.jl:686 [inlined]
[5] ODEProblem{true,tType,isinplace,P,F,K,PT} where PT where K where F where P where isinplace where tType(::ODESystem, ::Array{Float64,1}, ::Tuple{Float64,Int64}, ::Array{Float64,1}; version::Nothing, tgrad::Bool, jac::Bool, checkbounds::Bool, sparse::Bool, simplify::Bool, linenumbers::Bool, parallel::ModelingToolkit.SerialForm, eval_expression::Bool, kwargs::Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}) at C:\Users\claud\.julia\packages\ModelingToolkit\CPB4Z\src\systems\diffeqs\abstractodesystem.jl:270
[6] ODEProblem(::ODESystem, ::Array{Float64,1}, ::Vararg{Any,N} where N; kwargs::Base.Iterators.Pairs{Symbol,Bool,Tuple{Symbol},NamedTuple{(:jac,),Tuple{Bool}}}) at C:\Users\claud\.julia\packages\ModelingToolkit\CPB4Z\src\systems\diffeqs\abstractodesystem.jl:241
[7] top-level scope at In[26]:1
[8] include_string(::Function, ::Module, ::String, ::String) at .\loading.jl:1091
This seems to be related to what I get when trying to use auto_optimize:
_prob, _alg = auto_optimize(problem)
Try ModelingToolkitization
┌ Warning: ModelingToolkitization Approach Failed
└ @ AutoOptimize C:\Users\claud\.julia\packages\AutoOptimize\V2uMw\src\AutoOptimize.jl:65
BoundsError(723207.925, 2)
BoundsError: attempt to access Float64
at index [2]
Stacktrace:
[1] auto_optimize(::ODEProblem{Array{Float64,1},Tuple{Float64,Float64},true,Array{Float64,1},ODEFunction{true,typeof(SIRD_ac!),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}},DiffEqBase.StandardODEProblem}, ::Nothing; verbose::Bool, stiff::Bool, mtkify::Bool, sparsify::Bool, gpuify::Bool, static::Bool, gpup::Nothing) at C:\Users\claud\.julia\packages\AutoOptimize\V2uMw\src\AutoOptimize.jl:67
[2] auto_optimize at C:\Users\claud\.julia\packages\AutoOptimize\V2uMw\src\AutoOptimize.jl:37 [inlined] (repeats 2 times)
[3] top-level scope at In[31]:1
[4] include_string(::Function, ::Module, ::String, ::String) at .\loading.jl:1091
So what am I missing?
Also, would ModelingToolkit and/or AutoOptimize perform such optimization and/or best solve algorithm selection also for SDEProblem and EnsembleProblem built upon these?
I’m using:
(computationalEpi) pkg> st
Status `E:\IlMIoDrive\magistrale\2anno\primo_periodo\computationalEpi\Project.toml`
[ff3c4d4f] AutoOptimize v0.1.0 `https://github.com/SciML/AutoOptimize.jl#master`
[8f4d0f93] Conda v1.5.0
[071ae1c0] DiffEqGPU v1.8.0
[0c46a032] DifferentialEquations v6.15.0
[7073ff75] IJulia v1.23.0
[961ee093] ModelingToolkit v4.0.6
[d330b81b] PyPlot v2.9.0
Thank you very much.