explicit handling of LabelledArrays (and NamedTuples?)

[fkrauer]

I have an example, where LabelledArrays are not accepted. I'd love to use modelingtoolkitize, but I prefer to use labelled arrays to keep track of my parameters. Would be great if this could work. The modelingtoolkitized solution also loses the names for the states and returns x1..xn, which is a bit of a pity.

This code:

using DifferentialEquations
using StatsPlots
using Turing
using Distributions
using Random
using StatsBase
using LabelledArrays
using ModelingToolkit
#using DynamicHMC

Random.seed!(42)

# ODE model: simple SIR model with seasonally forced contact rate
function SIR!(du,u,p,t)

    # states
    (S, I, R) = u[1:3]
    N = S + I + R
  
    # params
    β = p.β
    η = p.η
    φ = p.φ
    ω = 1.0/p.ω
    μ = p.μ
    σ = p.σ
      
    # FOI
    βeff = β * (1.0+η*cos(2.0*π*(t-φ)/365.0))
    λ = βeff*I/N

    # change in states
    du[1] = (μ*N - λ*S - μ*S + ω*R)
    du[2] = (λ*S - σ*I - μ*I)
    du[3] = (σ*I - μ*R - ω*R)
    du[4] = (σ*I) # cumulative incidence

end
  
# Solver settings
tmin = 0.0
tmax = 10.0*365.0
tspan = (tmin, tmax)
solvsettings = (abstol = 1.0e-6, 
                reltol = 1.0e-3, 
                saveat = 7.0,
                solver = AutoTsit5(Rosenbrock23()))

# Initiate ODE problem
theta_fix =  [1.0/(80*365)]
theta_est =  [0.28, 0.07, 1.0/365.0, 1.0 ,1.0/5.0]
p = @LArray [theta_est; theta_fix] (:β, :η, :ω, :φ, :σ, :μ)
u0 = @LArray [9998.0,1.0,1.0,1.0] (:S,:I,:R,:C)

# Initiate ODE problem
problem = ODEProblem(SIR!,u0,tspan,p)
problem_acc = ODEProblem(modelingtoolkitize(problem), u0, tspan, p, jac=true, sparse=true)

sol = solve(problem_acc, 
            solvsettings.solver, 
            abstol=solvsettings.abstol, 
            reltol=solvsettings.reltol, 
            isoutofdomain=(u,p,t)->any(x->x<0.0,u), 
            saveat=solvsettings.saveat)

throws an error:

ERROR: LoadError: type Array has no field β

But this works:

function SIR!(du,u,p,t)

    # states
    (S, I, R) = u[1:3]
    N = S + I + R
  
    # params
    β, η, ω, φ, σ, μ = p
      
    # FOI
    βeff = β * (1.0+η*cos(2.0*π*(t-φ)/365.0))
    λ = βeff*I/N

    # change in states
    du[1] = (μ*N - λ*S - μ*S + ω*R)
    du[2] = (λ*S - σ*I - μ*I)
    du[3] = (σ*I - μ*R - ω*R)
    du[4] = (σ*I) # cumulative incidence

end

p = [theta_est; theta_fix]

#--> plug p in code above

Originally posted by @fkrauer in #1009 (comment)

[ChrisRackauckas]

Thanks for the example. The PR is almost complete: #1055

[fkrauer]

awesome, thanks

[jbiffl]

Hi All,

I'm not sure if I should open a new issue, but I'm having the same issues but with nested labelled vectors with the latest version (and previous versions). The code below gives pretty much the same error as the original example by fkrauer, but only when using nested LVectors.

using DifferentialEquations
using LabelledArrays
using ModelingToolkit

# ODE model: nonsense based on lorenz
function ODEFun!(du,u,p,t)

    # params in two layers does not work

    a = p.firstset.a
    b = p.firstset.b
    c = p.firstset.c

    d = p.secondset.d
    e = p.secondset.e
    f = p.secondset.f

    #params in one layer does work

    #=
    a = p.a
    b = p.b
    c = p.c
    d = p.d
    e = p.e
    f = p.f
    =#

    # change in states
    du[1] = a*(u[2] - u[1])+d
    du[2] = u[1]*(b-u[3]-u[2])+e
    du[3] = u[1]*u[2] - c*u[3]+f

end

# Solver settings
tmin = 0.0
tmax = 10.0
tspan = (tmin, tmax)
solvsettings = (abstol = 1.0e-6, 
                reltol = 1.0e-6, 
                saveat = 1.0,
                solver = Rodas4())

# ODE Parameters as nested LVector does not work
p = LVector(firstset = LVector(a = 10, b = 29, c = 8/3), secondset = LVector(d = 0, e = 0, f = 0)) #value of variables has no impact on error
# ODE Parameters as single layer LVector works
#p = LVector(a = 10, b = 28, c = 8/3, d = 0, e = 0, f = 0)
u0 = [1.0 0.0 0.0]

# Initiate ODE problem
problem = ODEProblem(ODEFun!,u0,tspan,p)
problem_acc = ODEProblem(modelingtoolkitize(problem), u0, tspan, p, jac=true, sparse=true)

sol = solve(problem_acc, 
            solvsettings.solver, 
            abstol=solvsettings.abstol, 
            reltol=solvsettings.reltol, 
            isoutofdomain=(u,p,t)->any(x->x<0.0,u), 
            saveat=solvsettings.saveat)

That code throws this error:

ERROR: LoadError: type Num has no field a
Stacktrace:
 [1] getproperty(x::Num, f::Symbol)
   @ Base .\Base.jl:33
 [2] ODEFun!(du::Matrix{Num}, u::Matrix{Num}, p::LArray{Num, 1, Vector{Num}, (:firstset, :secondset)}, t::Num)

When you inspect p while it is running, you obtain

2-element LArray{Num, 1, Vector{Num}, (:firstset, :secondset)}:
:firstset => firstset
:secondset => secondset

edit: fixed typo odes ->does, fixed wrong error message pasted

[ChrisRackauckas]

New issue.