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