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

**URL:** https://discourse.julialang.org/t/boundserror-attempt-to-access-float64-using-modelingtoolkit-jl-and-autooptimize-jl-to-optimize-a-differentialequations-jl-model/51000
**Category:** Optimization (Mathematical)
**Tags:** question, diffeq, optimization, modelingtoolkit
**Created:** [November 30, 2020, 4:17pm UTC](https://discourse.julialang.org/t/boundserror-attempt-to-access-float64-using-modelingtoolkit-jl-and-autooptimize-jl-to-optimize-a-differentialequations-jl-model/51000 "2020-11-30T16:17:12Z")
**Posts on this page:** 2
**Page:** 1

<div class="post-metadata">

### Author: ![claudio20497](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/claudio20497/32/36602_2.png) [@claudio20497](https://discourse.julialang.org/u/claudio20497)
#### Post date: [November 30, 2020, 4:17pm UTC](https://discourse.julialang.org/t/boundserror-attempt-to-access-float64-using-modelingtoolkit-jl-and-autooptimize-jl-to-optimize-a-differentialequations-jl-model/51000/1 "2020-11-30T16:17:12Z")

</div>

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:

```julia
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:

```julia
problem = ODEProblem(SIRD_ac!, ℬ, 𝒯, 𝒫)
sys = modelingtoolkitize(problem)

```

But when I try to rebuild the `ODEProblem` I get:

```julia
fast_problem = ODEProblem(sys,ℬ, 𝒯, 𝒫 )

```

```julia
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`:

```julia
_prob, _alg = auto_optimize(problem)

```

```julia
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:

```julia
(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.

---

<div class="post-metadata">

### Author: ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [December 1, 2020, 4:35am UTC](https://discourse.julialang.org/t/boundserror-attempt-to-access-float64-using-modelingtoolkit-jl-and-autooptimize-jl-to-optimize-a-differentialequations-jl-model/51000/2 "2020-12-01T04:35:19Z")

</div>

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](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:

```julia
0.000285 seconds (4.43 k allocations: 197.484 KiB)

```

after:

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

```

mostly due to the allocation of the Λ.