Skip to content

Commit 9a8ea44

Browse files
move another labelledarray
1 parent ae9ceb3 commit 9a8ea44

File tree

2 files changed

+49
-49
lines changed

2 files changed

+49
-49
lines changed

test/labelledarrays.jl

Lines changed: 48 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -41,3 +41,51 @@ d = LVector(x = 1.0, y = 2.0, z = 3.0)
4141
@test ff.jac(d, p, ForwardDiff.Dual(0.0, 1.0)) isa Array
4242
@inferred ff.jac(d, p, ForwardDiff.Dual(0.0, 1.0))
4343
@test eltype(ff.jac(d, p, ForwardDiff.Dual(0.0, 1.0))) <: ForwardDiff.Dual
44+
45+
## https://github.com/SciML/ModelingToolkit.jl/issues/1054
46+
using LabelledArrays
47+
using ModelingToolkit
48+
49+
# ODE model: simple SIR model with seasonally forced contact rate
50+
function SIR!(du, u, p, t)
51+
52+
# Unknowns
53+
(S, I, R) = u[1:3]
54+
N = S + I + R
55+
56+
# params
57+
β = p.β
58+
η = p.η
59+
φ = p.φ
60+
ω = 1.0 / p.ω
61+
μ = p.μ
62+
σ = p.σ
63+
64+
# FOI
65+
βeff = β * (1.0 + η * cos(2.0 * π * (t - φ) / 365.0))
66+
λ = βeff * I / N
67+
68+
# change in unknowns
69+
du[1] =* N - λ * S - μ * S + ω * R)
70+
du[2] =* S - σ * I - μ * I)
71+
du[3] =* I - μ * R - ω * R)
72+
du[4] =* I) # cumulative incidence
73+
end
74+
75+
# Solver settings
76+
tmin = 0.0
77+
tmax = 10.0 * 365.0
78+
tspan = (tmin, tmax)
79+
80+
# Initiate ODE problem
81+
theta_fix = [1.0 / (80 * 365)]
82+
theta_est = [0.28, 0.07, 1.0 / 365.0, 1.0, 1.0 / 5.0]
83+
p = @LArray [theta_est; theta_fix] (, , , , , )
84+
u0 = @LArray [9998.0, 1.0, 1.0, 1.0] (:S, :I, :R, :C)
85+
86+
# Initiate ODE problem
87+
problem = ODEProblem(SIR!, u0, tspan, p)
88+
sys = complete(modelingtoolkitize(problem))
89+
90+
@test all(isequal.(parameters(sys), getproperty.(@variables(β, η, ω, φ, σ, μ), :val)))
91+
@test all(isequal.(Symbol.(unknowns(sys)), Symbol.(@variables(S(t), I(t), R(t), C(t)))))

test/modelingtoolkitize.jl

Lines changed: 1 addition & 49 deletions
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,6 @@
11
using OrdinaryDiffEq, ModelingToolkit, DataStructures, Test
22
using Optimization, RecursiveArrayTools, OptimizationOptimJL
3-
using LabelledArrays, SymbolicIndexingInterface
3+
using SymbolicIndexingInterface
44
using ModelingToolkit: t_nounits as t, D_nounits as D
55
using SciMLBase: parameterless_type
66

@@ -210,54 +210,6 @@ tspan = (0.0, 1.0)
210210
prob = ODEProblem(k, x0, tspan)
211211
sys = modelingtoolkitize(prob)
212212

213-
## https://github.com/SciML/ModelingToolkit.jl/issues/1054
214-
using LabelledArrays
215-
using ModelingToolkit
216-
217-
# ODE model: simple SIR model with seasonally forced contact rate
218-
function SIR!(du, u, p, t)
219-
220-
# Unknowns
221-
(S, I, R) = u[1:3]
222-
N = S + I + R
223-
224-
# params
225-
β = p.β
226-
η = p.η
227-
φ = p.φ
228-
ω = 1.0 / p.ω
229-
μ = p.μ
230-
σ = p.σ
231-
232-
# FOI
233-
βeff = β * (1.0 + η * cos(2.0 * π * (t - φ) / 365.0))
234-
λ = βeff * I / N
235-
236-
# change in unknowns
237-
du[1] =* N - λ * S - μ * S + ω * R)
238-
du[2] =* S - σ * I - μ * I)
239-
du[3] =* I - μ * R - ω * R)
240-
du[4] =* I) # cumulative incidence
241-
end
242-
243-
# Solver settings
244-
tmin = 0.0
245-
tmax = 10.0 * 365.0
246-
tspan = (tmin, tmax)
247-
248-
# Initiate ODE problem
249-
theta_fix = [1.0 / (80 * 365)]
250-
theta_est = [0.28, 0.07, 1.0 / 365.0, 1.0, 1.0 / 5.0]
251-
p = @LArray [theta_est; theta_fix] (, , , , , )
252-
u0 = @LArray [9998.0, 1.0, 1.0, 1.0] (:S, :I, :R, :C)
253-
254-
# Initiate ODE problem
255-
problem = ODEProblem(SIR!, u0, tspan, p)
256-
sys = complete(modelingtoolkitize(problem))
257-
258-
@test all(isequal.(parameters(sys), getproperty.(@variables(β, η, ω, φ, σ, μ), :val)))
259-
@test all(isequal.(Symbol.(unknowns(sys)), Symbol.(@variables(S(t), I(t), R(t), C(t)))))
260-
261213
# https://github.com/SciML/ModelingToolkit.jl/issues/1158
262214

263215
function ode_prob(du, u, p::NamedTuple, t)

0 commit comments

Comments
 (0)