Skip to content

Commit 498b4e9

Browse files
TorkelETorkelE
committed
add problem input test
1 parent 9bc10f7 commit 498b4e9

File tree

2 files changed

+257
-1
lines changed

2 files changed

+257
-1
lines changed

‎test/runtests.jl

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -43,7 +43,8 @@ end
4343
@safetestset "PDE Construction Test" include("pde.jl")
4444
@safetestset "JumpSystem Test" include("jumpsystem.jl")
4545
@safetestset "Constraints Test" include("constraints.jl")
46-
@safetestset "Structure Interfacing Test" include("sciml_struct_interfacing.jl")
46+
@safetestset "SciML Problem Input Test" include("sciml_problem_inputs.jl")
47+
@safetestset "Structure Interfacing Test" include("sciml_struct_interfacing.jl")
4748
@safetestset "Reduction Test" include("reduction.jl")
4849
@safetestset "Split Parameters Test" include("split_parameters.jl")
4950
@safetestset "StaticArrays Test" include("static_arrays.jl")

‎test/sciml_problem_inputs.jl

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

0 commit comments

Comments
 (0)