BoundsError: attempt to access Float64 using ModelingToolkit.jl and AutoOptimize.jl to optimize a DifferentialEquations.jl model


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]


# 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]

 [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]

 [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 ``
  [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.


Thanks for the well-written report. It looks like we introduced a bug here in our ModelingToolkit v4 release. I fixed it in and will merge and tag the patch tomorrow morning. Also, I made this model into a regression test to catch this in the future.

P.S. The speedup you get from this is before:

0.000285 seconds (4.43 k allocations: 197.484 KiB)


0.000060 seconds (73 allocations: 13.453 KiB)

mostly due to the allocation of the Λ.