1
1
# ## Prepares Tests ###
2
2
3
3
# Fetch packages
4
- using ModelingToolkit, JumpProcesses, NonlinearSolve, OrdinaryDiffEq, SteadyStateDiffEq, StochasticDiffEq, Test
4
+ using ModelingToolkit, JumpProcesses, NonlinearSolve, OrdinaryDiffEq, SteadyStateDiffEq,
5
+ StochasticDiffEq, Test
5
6
using ModelingToolkit: t_nounits as t, D_nounits as D
6
7
7
8
# Sets rnd number.
@@ -12,19 +13,19 @@ seed = rand(rng, 1:100)
12
13
# ## Basic Tests ###
13
14
14
15
# Prepares a models and initial conditions/parameters (of different forms) to be used as problem inputs.
15
- begin
16
+ begin
16
17
# Prepare system components.
17
18
@parameters kp kd k1 k2= 0.5 Z0
18
- @variables X (t) Y (t) Z (t) = Z0
19
+ @variables X (t) Y (t) Z (t)= Z0
19
20
alg_eqs = [
20
- 0 ~ kp - k1* X + k2* Y - kd* X,
21
- 0 ~ - k1* Y + k1* X - k2* Y + k2* Z,
22
- 0 ~ k1* Y - k2* Z
21
+ 0 ~ kp - k1 * X + k2 * Y - kd * X,
22
+ 0 ~ - k1 * Y + k1 * X - k2 * Y + k2 * Z,
23
+ 0 ~ k1 * Y - k2 * Z
23
24
]
24
25
diff_eqs = [
25
- D (X) ~ kp - k1* X + k2* Y - kd* X,
26
- D (Y) ~ - k1* Y + k1* X - k2* Y + k2* Z,
27
- D (Z) ~ k1* Y - k2* Z
26
+ D (X) ~ kp - k1 * X + k2 * Y - kd * X,
27
+ D (Y) ~ - k1 * Y + k1 * X - k2 * Y + k2 * Z,
28
+ D (Z) ~ k1 * Y - k2 * Z
28
29
]
29
30
noise_eqs = fill (0.01 , 3 , 6 )
30
31
jumps = [
38
39
39
40
# Create systems (without structural_simplify, since that might modify systems to affect intended tests).
40
41
osys = complete (ODESystem (diff_eqs, t; name = :osys ))
41
- ssys = complete (SDESystem (diff_eqs, noise_eqs, t, [X, Y, Z], [kp, kd, k1, k2]; name = :ssys ))
42
+ ssys = complete (SDESystem (
43
+ diff_eqs, noise_eqs, t, [X, Y, Z], [kp, kd, k1, k2]; name = :ssys ))
42
44
jsys = complete (JumpSystem (jumps, t, [X, Y, Z], [kp, kd, k1, k2]; name = :jsys ))
43
45
nsys = complete (NonlinearSystem (alg_eqs; name = :nsys ))
44
46
60
62
(osys. X => 4 , osys. Y => 5 ),
61
63
# Tuples providing default values.
62
64
(X => 4 , Y => 5 , Z => 10 ),
63
- (osys. X => 4 , osys. Y => 5 , osys. Z => 10 ),
65
+ (osys. X => 4 , osys. Y => 5 , osys. Z => 10 )
64
66
]
65
67
tspan = (0.0 , 10.0 )
66
68
p_alts = [
@@ -75,18 +77,19 @@ begin
75
77
Dict ([osys. kp => 1.0 , osys. kd => 0.1 , osys. k1 => 0.25 , osys. Z0 => 10 ]),
76
78
# Dicts providing default values.
77
79
Dict ([kp => 1.0 , kd => 0.1 , k1 => 0.25 , k2 => 0.5 , Z0 => 10 ]),
78
- Dict ([osys. kp => 1.0 , osys. kd => 0.1 , osys. k1 => 0.25 , osys. k2 => 0.5 , osys. Z0 => 10 ]),
80
+ Dict ([osys. kp => 1.0 , osys. kd => 0.1 , osys. k1 => 0.25 ,
81
+ osys. k2 => 0.5 , osys. Z0 => 10 ]),
79
82
# Tuples not providing default values.
80
83
(kp => 1.0 , kd => 0.1 , k1 => 0.25 , Z0 => 10 ),
81
84
(osys. kp => 1.0 , osys. kd => 0.1 , osys. k1 => 0.25 , osys. Z0 => 10 ),
82
85
# Tuples providing default values.
83
86
(kp => 1.0 , kd => 0.1 , k1 => 0.25 , k2 => 0.5 , Z0 => 10 ),
84
- (osys. kp => 1.0 , osys. kd => 0.1 , osys. k1 => 0.25 , osys. k2 => 0.5 , osys. Z0 => 10 ),
87
+ (osys. kp => 1.0 , osys. kd => 0.1 , osys. k1 => 0.25 , osys. k2 => 0.5 , osys. Z0 => 10 )
85
88
]
86
89
end
87
90
88
91
# Perform ODE simulations (singular and ensemble).
89
- let
92
+ let
90
93
# Creates normal and ensemble problems.
91
94
base_oprob = ODEProblem (osys, u0_alts[1 ], tspan, p_alts[1 ])
92
95
base_sol = solve (base_oprob, Tsit5 (); saveat = 1.0 )
97
100
# test failure.
98
101
for u0 in u0_alts, p in p_alts
99
102
oprob = remake (base_oprob; u0, p)
100
- # @test base_sol == solve(oprob, Tsit5(); saveat = 1.0)
103
+ # @test base_sol == solve(oprob, Tsit5(); saveat = 1.0)
101
104
eprob = remake (base_eprob; u0, p)
102
- # @test base_esol == solve(eprob, Tsit5(); trajectories = 2, saveat = 1.0)
105
+ # @test base_esol == solve(eprob, Tsit5(); trajectories = 2, saveat = 1.0)
103
106
end
104
107
end
105
108
106
109
# Perform SDE simulations (singular and ensemble).
107
- let
110
+ let
108
111
# Creates normal and ensemble problems.
109
112
base_sprob = SDEProblem (ssys, u0_alts[1 ], tspan, p_alts[1 ])
110
113
base_sol = solve (base_sprob, ImplicitEM (); seed, saveat = 1.0 )
@@ -114,15 +117,15 @@ let
114
117
# Simulates problems for all input types, checking that identical solutions are found.
115
118
@test_broken false # first remake in subsequent test yields a `ERROR: type Nothing has no field portion`.
116
119
for u0 in u0_alts, p in p_alts
117
- # sprob = remake(base_sprob; u0, p)
118
- # @test base_sol == solve(sprob, ImplicitEM(); seed, saveat = 1.0)
119
- # eprob = remake(base_eprob; u0, p)
120
- # @test base_esol == solve(eprob, ImplicitEM(); seed, trajectories = 2, saveat = 1.0)
120
+ # sprob = remake(base_sprob; u0, p)
121
+ # @test base_sol == solve(sprob, ImplicitEM(); seed, saveat = 1.0)
122
+ # eprob = remake(base_eprob; u0, p)
123
+ # @test base_esol == solve(eprob, ImplicitEM(); seed, trajectories = 2, saveat = 1.0)
121
124
end
122
125
end
123
126
124
127
# Perform jump simulations (singular and ensemble).
125
- let
128
+ let
126
129
# Creates normal and ensemble problems.
127
130
base_dprob = DiscreteProblem (jsys, u0_alts[1 ], tspan, p_alts[1 ])
128
131
base_jprob = JumpProblem (jsys, base_dprob, Direct (); rng)
@@ -133,10 +136,10 @@ let
133
136
# Simulates problems for all input types, checking that identical solutions are found.
134
137
@test_broken false # first remake in subsequent test yields a `ERROR: type Nothing has no field portion`.
135
138
for u0 in u0_alts, p in p_alts
136
- # jprob = remake(base_jprob; u0, p)
137
- # @test base_sol == solve(base_jprob, SSAStepper(); seed, saveat = 1.0)
138
- # eprob = remake(base_eprob; u0, p)
139
- # @test base_esol == solve(eprob, SSAStepper(); seed, trajectories = 2, saveat = 1.0)
139
+ # jprob = remake(base_jprob; u0, p)
140
+ # @test base_sol == solve(base_jprob, SSAStepper(); seed, saveat = 1.0)
141
+ # eprob = remake(base_eprob; u0, p)
142
+ # @test base_esol == solve(eprob, SSAStepper(); seed, trajectories = 2, saveat = 1.0)
140
143
end
141
144
end
142
145
@@ -148,12 +151,12 @@ let
148
151
# test failure.
149
152
for u0 in u0_alts, p in p_alts
150
153
nlprob = remake (base_nlprob; u0, p)
151
- # @test base_sol == solve(nlprob, NewtonRaphson())
154
+ # @test base_sol == solve(nlprob, NewtonRaphson())
152
155
end
153
156
end
154
157
155
158
# Perform steady state simulations (singular and ensemble).
156
- let
159
+ let
157
160
# Creates normal and ensemble problems.
158
161
base_ssprob = SteadyStateProblem (osys, u0_alts[1 ], p_alts[1 ])
159
162
base_sol = solve (base_ssprob, DynamicSS (Tsit5 ()))
170
173
end
171
174
end
172
175
173
-
174
176
# ## Checks Errors On Faulty Inputs ###
175
177
176
178
# Checks various erroneous problem inputs, ensuring that these throw errors.
@@ -179,12 +181,12 @@ let
179
181
@parameters k1 k2 k3
180
182
@variables X1 (t) X2 (t) X3 (t)
181
183
alg_eqs = [
182
- 0 ~ - k1* X1 + k2* X2,
183
- 0 ~ k1* X1 - k2* X2
184
+ 0 ~ - k1 * X1 + k2 * X2,
185
+ 0 ~ k1 * X1 - k2 * X2
184
186
]
185
187
diff_eqs = [
186
- D (X1) ~ - k1* X1 + k2* X2,
187
- D (X2) ~ k1* X1 - k2* X2
188
+ D (X1) ~ - k1 * X1 + k2 * X2,
189
+ D (X2) ~ k1 * X1 - k2 * X2
188
190
]
189
191
noise_eqs = fill (0.01 , 2 , 2 )
190
192
jumps = [
215
217
# Contain an additional value.
216
218
[X1 => 1 , X2 => 2 , X3 => 3 ],
217
219
Dict ([X1 => 1 , X2 => 2 , X3 => 3 ]),
218
- (X1 => 1 , X2 => 2 , X3 => 3 ),
220
+ (X1 => 1 , X2 => 2 , X3 => 3 )
219
221
]
220
222
ps_invalid = [
221
223
# Missing a value.
@@ -228,16 +230,18 @@ let
228
230
# Contain an additional value.
229
231
[k1 => 1.0 , k2 => 2.0 , k3 => 3.0 ],
230
232
Dict ([k1 => 1.0 , k2 => 2.0 , k3 => 3.0 ]),
231
- (k1 => 1.0 , k2 => 2.0 , k3 => 3.0 ),
233
+ (k1 => 1.0 , k2 => 2.0 , k3 => 3.0 )
232
234
]
233
235
234
236
# Loops through all potential parameter sets, checking their inputs yield errors.
235
237
# Broken tests are due to this issue: https://github.com/SciML/ModelingToolkit.jl/issues/2779
236
238
for ps in [ps_valid; ps_invalid], u0 in [u0_valid; u0s_invalid]
237
239
# Handles problems with/without tspan separately. Special check ensuring that valid inputs passes.
238
- for (xsys, XProblem) in zip ([osys, ssys, jsys], [ODEProblem, SDEProblem, DiscreteProblem])
240
+ for (xsys, XProblem) in zip (
241
+ [osys, ssys, jsys], [ODEProblem, SDEProblem, DiscreteProblem])
239
242
if (ps == ps_valid) && (u0 == u0_valid)
240
- XProblem (xsys, u0, (0.0 , 1.0 ), ps); @test true ;
243
+ XProblem (xsys, u0, (0.0 , 1.0 ), ps)
244
+ @test true
241
245
else
242
246
@test_broken false
243
247
continue
246
250
end
247
251
for (xsys, XProblem) in zip ([nsys, osys], [NonlinearProblem, SteadyStateProblem])
248
252
if (ps == ps_valid) && (u0 == u0_valid)
249
- XProblem (xsys, u0, ps); @test true ;
253
+ XProblem (xsys, u0, ps)
254
+ @test true
250
255
else
251
256
@test_broken false
252
257
continue
@@ -258,21 +263,21 @@ end
258
263
259
264
# ## Vector Parameter/Variable Inputs ###
260
265
261
- begin
266
+ begin
262
267
# Declares the model (with vector species/parameters, with/without default values, and observables).
263
- @variables X (t)[1 : 2 ] Y (t)[1 : 2 ] = [10.0 , 20.0 ] XY (t)[1 : 2 ]
264
- @parameters p[1 : 2 ] d[1 : 2 ] = [0.2 , 0.5 ]
268
+ @variables X (t)[1 : 2 ] Y (t)[1 : 2 ]= [10.0 , 20.0 ] XY (t)[1 : 2 ]
269
+ @parameters p[1 : 2 ] d[1 : 2 ]= [0.2 , 0.5 ]
265
270
alg_eqs = [
266
- 0 ~ p[1 ] - d[1 ]* X[1 ],
267
- 0 ~ p[2 ] - d[2 ]* X[2 ],
268
- 0 ~ p[1 ] - d[1 ]* Y[1 ],
269
- 0 ~ p[2 ] - d[2 ]* Y[2 ],
271
+ 0 ~ p[1 ] - d[1 ] * X[1 ],
272
+ 0 ~ p[2 ] - d[2 ] * X[2 ],
273
+ 0 ~ p[1 ] - d[1 ] * Y[1 ],
274
+ 0 ~ p[2 ] - d[2 ] * Y[2 ]
270
275
]
271
276
diff_eqs = [
272
- D (X[1 ]) ~ p[1 ] - d[1 ]* X[1 ],
273
- D (X[2 ]) ~ p[2 ] - d[2 ]* X[2 ],
274
- D (Y[1 ]) ~ p[1 ] - d[1 ]* Y[1 ],
275
- D (Y[2 ]) ~ p[2 ] - d[2 ]* Y[2 ],
277
+ D (X[1 ]) ~ p[1 ] - d[1 ] * X[1 ],
278
+ D (X[2 ]) ~ p[2 ] - d[2 ] * X[2 ],
279
+ D (Y[1 ]) ~ p[1 ] - d[1 ] * Y[1 ],
280
+ D (Y[2 ]) ~ p[2 ] - d[2 ] * Y[2 ]
276
281
]
277
282
noise_eqs = fill (0.01 , 4 , 8 )
278
283
jumps = [
@@ -289,8 +294,10 @@ begin
289
294
290
295
# Create systems (without structural_simplify, since that might modify systems to affect intended tests).
291
296
osys = complete (ODESystem (diff_eqs, t; observed, name = :osys ))
292
- ssys = complete (SDESystem (diff_eqs, noise_eqs, t, [X[1 ], X[2 ], Y[1 ], Y[2 ]], [p, d]; observed, name = :ssys ))
293
- jsys = complete (JumpSystem (jumps, t, [X[1 ], X[2 ], Y[1 ], Y[2 ]], [p, d]; observed, name = :jsys ))
297
+ ssys = complete (SDESystem (
298
+ diff_eqs, noise_eqs, t, [X[1 ], X[2 ], Y[1 ], Y[2 ]], [p, d]; observed, name = :ssys ))
299
+ jsys = complete (JumpSystem (
300
+ jumps, t, [X[1 ], X[2 ], Y[1 ], Y[2 ]], [p, d]; observed, name = :jsys ))
294
301
nsys = complete (NonlinearSystem (alg_eqs; observed, name = :nsys ))
295
302
296
303
# Declares various u0 versions (scalarised and vector forms).
@@ -354,7 +361,7 @@ begin
354
361
end
355
362
356
363
# Perform ODE simulations (singular and ensemble).
357
- let
364
+ let
358
365
# Creates normal and ensemble problems.
359
366
base_oprob = ODEProblem (osys, u0_alts_vec[1 ], tspan, p_alts_vec[1 ])
360
367
base_sol = solve (base_oprob, Tsit5 (); saveat = 1.0 )
372
379
end
373
380
374
381
# Perform SDE simulations (singular and ensemble).
375
- let
382
+ let
376
383
# Creates normal and ensemble problems.
377
384
base_sprob = SDEProblem (ssys, u0_alts_vec[1 ], tspan, p_alts_vec[1 ])
378
385
base_sol = solve (base_sprob, ImplicitEM (); seed, saveat = 1.0 )
391
398
392
399
# Perform jump simulations (singular and ensemble).
393
400
# Fails. At least partially related to https://github.com/SciML/ModelingToolkit.jl/issues/2804.
394
- @test_broken let
401
+ @test_broken let
395
402
# Creates normal and ensemble problems.
396
403
base_dprob = DiscreteProblem (jsys, u0_alts_vec[1 ], tspan, p_alts_vec[1 ])
397
404
base_jprob = JumpProblem (jsys, base_dprob, Direct (); rng)
422
429
423
430
# Perform steady state simulations (singular and ensemble).
424
431
# Fails. At least partially related to https://github.com/SciML/ModelingToolkit.jl/issues/2804.
425
- let
432
+ let
426
433
# Creates normal and ensemble problems.
427
434
base_ssprob = SteadyStateProblem (osys, u0_alts_vec[1 ], p_alts_vec[1 ])
428
435
base_sol = solve (base_ssprob, DynamicSS (Tsit5 ()))
0 commit comments