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

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)]) for i in 1:size(C)]

# System of equations
@. dS = -Λ*S
@. dI = Λ*S - ((1-δ)*λ_R + δ*λ_D)*I
@. dR = λ_R*(1-δ)*I
@. dD = λ_D*δ*I
@. dD_tot = dD+dD+dD+dD

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 

Stacktrace:
 indexed_iterate(::Float64, ::Int64, ::Nothing) at .\tuple.jl:90
 (::ModelingToolkit.var"#208#210"{ODESystem})(::Float64) at .\none:0
 iterate at .\generator.jl:47 [inlined]
 collect at .\array.jl:686 [inlined]
 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
 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
 top-level scope at In:1
``````

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 

Stacktrace:
 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
 auto_optimize at C:\Users\claud\.julia\packages\AutoOptimize\V2uMw\src\AutoOptimize.jl:37 [inlined] (repeats 2 times)
 top-level scope at In:1
``````

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.

2 Likes

Thanks for the well-written report. It looks like we introduced a bug here in our ModelingToolkit v4 release. I fixed it in https://github.com/SciML/ModelingToolkit.jl/pull/677 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)
``````

after:

``````0.000060 seconds (73 allocations: 13.453 KiB)
``````

mostly due to the allocation of the Λ.

3 Likes