1
+ # ! format: off
2
+
3
+ # ## Prepares Tests ###
4
+
5
+ # Fetch packages
6
+ using Catalyst, JumpProcesses, NonlinearSolve, OrdinaryDiffEq, SteadyStateDiffEq, StochasticDiffEq, Test
7
+
8
+ # Sets rnd number.
9
+ using StableRNGs
10
+ rng = StableRNG (12345 )
11
+ seed = rand (rng, 1 : 100 )
12
+
13
+
14
+ # ## Basic Tests ###
15
+
16
+ # Prepares a models and initial conditions/parameters (of different forms) to be used as problem inputs.
17
+ begin
18
+ model = @reaction_network begin
19
+ @species Z (t) = Z0
20
+ @parameters k2= 0.5 Z0
21
+ (kp,kd), 0 <--> X
22
+ (k1,k2), X <--> Y
23
+ (k1,k2), Y <--> Z
24
+ end
25
+ @unpack X, Y, Z, kp, kd, k1, k2, Z0 = model
26
+
27
+ u0_alts = [
28
+ # Vectors not providing default values.
29
+ [X => 4 , Y => 5 ],
30
+ [model. X => 4 , model. Y => 5 ],
31
+ [:X => 4 , :Y => 5 ],
32
+ # Vectors providing default values.
33
+ [X => 4 , Y => 5 , Z => 10 ],
34
+ [model. X => 4 , model. Y => 5 , model. Z => 10 ],
35
+ [:X => 4 , :Y => 5 , :Z => 10 ],
36
+ # Dicts not providing default values.
37
+ Dict ([X => 4 , Y => 5 ]),
38
+ Dict ([model. X => 4 , model. Y => 5 ]),
39
+ Dict ([:X => 4 , :Y => 5 ]),
40
+ # Dicts providing default values.
41
+ Dict ([X => 4 , Y => 5 , Z => 10 ]),
42
+ Dict ([model. X => 4 , model. Y => 5 , model. Z => 10 ]),
43
+ Dict ([:X => 4 , :Y => 5 , :Z => 10 ]),
44
+ # Tuples not providing default values.
45
+ (X => 4 , Y => 5 ),
46
+ (model. X => 4 , model. Y => 5 ),
47
+ (:X => 4 , :Y => 5 ),
48
+ # Tuples providing default values.
49
+ (X => 4 , Y => 5 , Z => 10 ),
50
+ (model. X => 4 , model. Y => 5 , model. Z => 10 ),
51
+ (:X => 4 , :Y => 5 , :Z => 10 )
52
+ ]
53
+ tspan = (0.0 , 10.0 )
54
+ p_alts = [
55
+ # Vectors not providing default values.
56
+ [kp => 1.0 , kd => 0.1 , k1 => 0.25 , Z0 => 10 ],
57
+ [model. kp => 1.0 , model. kd => 0.1 , model. k1 => 0.25 , model. Z0 => 10 ],
58
+ [:kp => 1.0 , :kd => 0.1 , :k1 => 0.25 , :Z0 => 10 ],
59
+ # Vectors providing default values.
60
+ [kp => 1.0 , kd => 0.1 , k1 => 0.25 , k2 => 0.5 , Z0 => 10 ],
61
+ [model. kp => 1.0 , model. kd => 0.1 , model. k1 => 0.25 , model. k2 => 0.5 , model. Z0 => 10 ],
62
+ [:kp => 1.0 , :kd => 0.1 , :k1 => 0.25 , :k2 => 0.5 , :Z0 => 10 ],
63
+ # Dicts not providing default values.
64
+ Dict ([kp => 1.0 , kd => 0.1 , k1 => 0.25 , Z0 => 10 ]),
65
+ Dict ([model. kp => 1.0 , model. kd => 0.1 , model. k1 => 0.25 , model. Z0 => 10 ]),
66
+ Dict ([:kp => 1.0 , :kd => 0.1 , :k1 => 0.25 , :Z0 => 10 ]),
67
+ # Dicts providing default values.
68
+ Dict ([kp => 1.0 , kd => 0.1 , k1 => 0.25 , k2 => 0.5 , Z0 => 10 ]),
69
+ Dict ([model. kp => 1.0 , model. kd => 0.1 , model. k1 => 0.25 , model. k2 => 0.5 , model. Z0 => 10 ]),
70
+ Dict ([:kp => 1.0 , :kd => 0.1 , :k1 => 0.25 , :k2 => 0.5 , :Z0 => 10 ]),
71
+ # Tuples not providing default values.
72
+ (kp => 1.0 , kd => 0.1 , k1 => 0.25 , Z0 => 10 ),
73
+ (model. kp => 1.0 , model. kd => 0.1 , model. k1 => 0.25 , model. Z0 => 10 ),
74
+ (:kp => 1.0 , :kd => 0.1 , :k1 => 0.25 , :Z0 => 10 ),
75
+ # Tuples providing default values.
76
+ (kp => 1.0 , kd => 0.1 , k1 => 0.25 , k2 => 0.5 , Z0 => 10 ),
77
+ (model. kp => 1.0 , model. kd => 0.1 , model. k1 => 0.25 , model. k2 => 0.5 , model. Z0 => 10 ),
78
+ (:kp => 1.0 , :kd => 0.1 , :k1 => 0.25 , :k2 => 0.5 , :Z0 => 10 ),
79
+ ]
80
+ end
81
+
82
+ # Perform ODE simulations (singular and ensemble).
83
+ let
84
+ # Creates normal and ensemble problems.
85
+ base_oprob = ODEProblem (model, u0_alts[1 ], tspan, p_alts[1 ])
86
+ base_sol = solve (base_oprob, Tsit5 (); saveat = 1.0 )
87
+ base_eprob = EnsembleProblem (base_oprob)
88
+ base_esol = solve (base_eprob, Tsit5 (); trajectories = 2 , saveat = 1.0 )
89
+
90
+ # Simulates problems for all input types, checking that identical solutions are found.
91
+ for u0 in u0_alts, p in p_alts
92
+ oprob = remake (base_oprob; u0, p)
93
+ @test base_sol == solve (oprob, Tsit5 (); saveat = 1.0 )
94
+ eprob = remake (base_eprob; u0, p)
95
+ @test base_esol == solve (eprob, Tsit5 (); trajectories = 2 , saveat = 1.0 )
96
+ end
97
+ end
98
+
99
+ # Perform SDE simulations (singular and ensemble).
100
+ let
101
+ # Creates normal and ensemble problems.
102
+ base_sprob = SDEProblem (model, u0_alts[1 ], tspan, p_alts[1 ])
103
+ base_sol = solve (base_sprob, ImplicitEM (); seed, saveat = 1.0 )
104
+ base_eprob = EnsembleProblem (base_sprob)
105
+ base_esol = solve (base_eprob, ImplicitEM (); seed, trajectories = 2 , saveat = 1.0 )
106
+
107
+ # Simulates problems for all input types, checking that identical solutions are found.
108
+ for u0 in u0_alts, p in p_alts
109
+ sprob = remake (base_sprob; u0, p)
110
+ @test base_sol == solve (sprob, ImplicitEM (); seed, saveat = 1.0 )
111
+ eprob = remake (base_eprob; u0, p)
112
+ @test base_esol == solve (eprob, ImplicitEM (); seed, trajectories = 2 , saveat = 1.0 )
113
+ end
114
+ end
115
+
116
+ # Perform jump simulations (singular and ensemble).
117
+ let
118
+ # Creates normal and ensemble problems.
119
+ base_dprob = DiscreteProblem (model, u0_alts[1 ], tspan, p_alts[1 ])
120
+ base_jprob = JumpProblem (model, base_dprob, Direct (); rng)
121
+ base_sol = solve (base_jprob, SSAStepper (); seed, saveat = 1.0 )
122
+ base_eprob = EnsembleProblem (base_jprob)
123
+ base_esol = solve (base_eprob, SSAStepper (); seed, trajectories = 2 , saveat = 1.0 )
124
+
125
+ # Simulates problems for all input types, checking that identical solutions are found.
126
+ for u0 in u0_alts, p in p_alts
127
+ jprob = remake (base_jprob; u0, p)
128
+ @test base_sol == solve (base_jprob, SSAStepper (); seed, saveat = 1.0 )
129
+ eprob = remake (base_eprob; u0, p)
130
+ @test base_esol == solve (eprob, SSAStepper (); seed, trajectories = 2 , saveat = 1.0 )
131
+ end
132
+ end
133
+
134
+ # Solves a nonlinear problem (EnsembleProblems are not possible for these).
135
+ let
136
+ base_nlprob = NonlinearProblem (model, u0_alts[1 ], p_alts[1 ])
137
+ base_sol = solve (base_nlprob, NewtonRaphson ())
138
+ for u0 in u0_alts, p in p_alts
139
+ nlprob = remake (base_nlprob; u0, p)
140
+ @test base_sol == solve (nlprob, NewtonRaphson ())
141
+ end
142
+ end
143
+
144
+ # Perform steady state simulations (singular and ensemble).
145
+ let
146
+ # Creates normal and ensemble problems.
147
+ base_ssprob = SteadyStateProblem (model, u0_alts[1 ], p_alts[1 ])
148
+ base_sol = solve (base_ssprob, DynamicSS (Tsit5 ()))
149
+ base_eprob = EnsembleProblem (base_ssprob)
150
+ base_esol = solve (base_eprob, DynamicSS (Tsit5 ()); trajectories = 2 )
151
+
152
+ # Simulates problems for all input types, checking that identical solutions are found.
153
+ for u0 in u0_alts, p in p_alts
154
+ ssprob = remake (base_ssprob; u0, p)
155
+ @test base_sol == solve (ssprob, DynamicSS (Tsit5 ()))
156
+ eprob = remake (base_eprob; u0, p)
157
+ @test base_esol == solve (eprob, DynamicSS (Tsit5 ()); trajectories = 2 )
158
+ end
159
+ end
160
+
161
+
162
+ # ## Checks Errors On Faulty Inputs ###
163
+
164
+ # Checks various erroneous problem inputs, ensuring that these throw errors.
165
+ let
166
+ # Declares the model.
167
+ rn = @reaction_network begin
168
+ (k1,k2), X1 <--> X2
169
+ end
170
+ @unpack k1, k2, X1, X2 = rn
171
+ t = default_t ()
172
+ @species X3 (t)
173
+ @parameters k3
174
+
175
+ # Declares valid initial conditions and parameter values
176
+ u0_valid = [X1 => 1 , X2 => 2 ]
177
+ ps_valid = [k1 => 0.5 , k2 => 0.1 ]
178
+
179
+ # Declares invalid initial conditions and parameters. This includes both cases where values are
180
+ # missing, or additional ones are given. Includes vector/Tuple/Dict forms.
181
+ u0s_invalid = [
182
+ # Missing a value.
183
+ [X1 => 1 ],
184
+ [rn. X1 => 1 ],
185
+ [:X1 => 1 ],
186
+ Dict ([X1 => 1 ]),
187
+ Dict ([rn. X1 => 1 ]),
188
+ Dict ([:X1 => 1 ]),
189
+ (X1 => 1 ),
190
+ (rn. X1 => 1 ),
191
+ (:X1 => 1 ),
192
+ # Contain an additional value.
193
+ [X1 => 1 , X2 => 2 , X3 => 3 ],
194
+ [:X1 => 1 , :X2 => 2 , :X3 => 3 ],
195
+ Dict ([X1 => 1 , X2 => 2 , X3 => 3 ]),
196
+ Dict ([:X1 => 1 , :X2 => 2 , :X3 => 3 ]),
197
+ (X1 => 1 , X2 => 2 , X3 => 3 ),
198
+ (:X1 => 1 , :X2 => 2 , :X3 => 3 )
199
+ ]
200
+ ps_invalid = [
201
+ # Missing a value.
202
+ [k1 => 1.0 ],
203
+ [rn. k1 => 1.0 ],
204
+ [:k1 => 1.0 ],
205
+ Dict ([k1 => 1.0 ]),
206
+ Dict ([rn. k1 => 1.0 ]),
207
+ Dict ([:k1 => 1.0 ]),
208
+ (k1 => 1.0 ),
209
+ (rn. k1 => 1.0 ),
210
+ (:k1 => 1.0 ),
211
+ # Contain an additional value.
212
+ [k1 => 1.0 , k2 => 2.0 , k3 => 3.0 ],
213
+ [:k1 => 1.0 , :k2 => 2.0 , :k3 => 3.0 ],
214
+ Dict ([k1 => 1.0 , k2 => 2.0 , k3 => 3.0 ]),
215
+ Dict ([:k1 => 1.0 , :k2 => 2.0 , :k3 => 3.0 ]),
216
+ (k1 => 1.0 , k2 => 2.0 , k3 => 3.0 ),
217
+ (:k1 => 1.0 , :k2 => 2.0 , :k3 => 3.0 )
218
+ ]
219
+
220
+ # Loops through all potential parameter sets, checking their inputs yield errors.
221
+ for ps in [ps_valid; ps_invalid], u0 in [u0_valid; u0s_invalid]
222
+ # Handles problems with/without tspan separately. Special check ensuring that valid inputs passes.
223
+ for XProblem in [ODEProblem, SDEProblem, DiscreteProblem]
224
+ if (ps == ps_valid) && (u0 == u0_valid)
225
+ XProblem (rn, u0, (0.0 , 1.0 ), ps); @test true ;
226
+ else
227
+ # Several of these cases do not throw errors (https://github.com/SciML/ModelingToolkit.jl/issues/2624).
228
+ @test_broken false
229
+ continue
230
+ @test_throws Exception XProblem (rn, u0, (0.0 , 1.0 ), ps)
231
+ end
232
+ end
233
+ for XProblem in [NonlinearProblem, SteadyStateProblem]
234
+ if (ps == ps_valid) && (u0 == u0_valid)
235
+ XProblem (rn, u0, ps); @test true ;
236
+ else
237
+ @test_broken false
238
+ continue
239
+ @test_throws Exception XProblem (rn, u0, ps)
240
+ end
241
+ end
242
+ end
243
+ end
0 commit comments