Skip to content

Commit ef878c3

Browse files
AayushSabharwalAayushSabharwal
committed
test: test implementation of SII/SciMLBase discrete saving interface
1 parent 46b493a commit ef878c3

File tree

3 files changed

+213
-89
lines changed

3 files changed

+213
-89
lines changed

‎test/mtkparameters.jl

Lines changed: 62 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -224,3 +224,65 @@ function loss(x)
224224
end
225225

226226
@test_nowarn ForwardDiff.gradient(loss, collect(tunables))
227+
228+
# Parameter timeseries
229+
dt = 0.1
230+
dt2 = 0.2
231+
@variables x(t)=0 y(t)=0 u(t)=0 yd1(t)=0 ud1(t)=0 yd2(t)=0 ud2(t)=0
232+
@parameters kp=1 r=1
233+
234+
eqs = [
235+
# controller (time discrete part `dt=0.1`)
236+
yd1 ~ Sample(t, dt)(y)
237+
ud1 ~ kp * (r - yd1)
238+
# controller (time discrete part `dt=0.2`)
239+
yd2 ~ Sample(t, dt2)(y)
240+
ud2 ~ kp * (r - yd2)
241+
242+
# plant (time continuous part)
243+
u ~ Hold(ud1) + Hold(ud2)
244+
D(x) ~ -x + u
245+
y ~ x]
246+
247+
@mtkbuild cl = ODESystem(eqs, t)
248+
ps = MTKParameters(cl, [kp => 1.0])
249+
with_updated_parameter_timeseries_values(
250+
ps, 1 => ModelingToolkit.NestedGetIndex(([5.0, 10.0],)))
251+
@test ps.discrete[1][1] == [5.0, 10.0]
252+
with_updated_parameter_timeseries_values(
253+
ps, 1 => ModelingToolkit.NestedGetIndex(([3.0, 30.0],)),
254+
2 => ModelingToolkit.NestedGetIndex(([4.0, 40.0],)))
255+
@test ps.discrete[1][1] == [3.0, 30.0]
256+
@test ps.discrete[2][1] == [4.0, 40.0]
257+
@test SciMLBase.get_saveable_values(ps, 1) == ModelingToolkit.NestedGetIndex(ps.discrete[1])
258+
259+
# With multiple types and clocks
260+
@variables x(t) xd1(t) xd2(t) flag(t)::Bool yd1(t) yd2(t) yc1(t) yc2(t)
261+
dt = 0.1
262+
k1 = ShiftIndex(t, dt)
263+
ssc = ModelingToolkit.SolverStepClock(t)
264+
k2 = ShiftIndex(ssc)
265+
266+
eqs = [
267+
flag ~ ~flag(k1 - 1),
268+
xd1 ~ Sample(t, dt)(x),
269+
yd1 ~ ifelse(flag, xd1, yd1(k1 - 1)), xd2 ~ Sample(ssc)(x),
270+
yd2 ~ yd2(k2 - 1) + xd2, yc1 ~ Hold(yd1),
271+
yc2 ~ Hold(yd2),
272+
D(x) ~ yc1 + yc2
273+
]
274+
@mtkbuild sys = ODESystem(eqs, t)
275+
ps = MTKParameters(sys,
276+
[flag => true, yd1 => ifelse(flag, Sample(t, dt)(x), 1.0),
277+
yd2 => 2.0 + Sample(ssc)(x), Sample(t, dt)(x) => x,
278+
Sample(ssc)(x) => x, Hold(yd1) => yd1, Hold(yd2) => yd2],
279+
[x => 3.0])
280+
@test SciMLBase.get_saveable_values(ps, 1).x isa Tuple{Vector{Float64}, Vector{Bool}}
281+
@test length(ps.discrete[1][1]) == 3
282+
@test length(ps.discrete[1][2]) == 1
283+
@test length(ps.discrete[2][1]) == 3
284+
@test length(ps.discrete[2][2]) == 0
285+
with_updated_parameter_timeseries_values(
286+
ps, 1 => ModelingToolkit.NestedGetIndex(([10.0, 11.0, 12.0], [false])))
287+
@test ps.discrete[1][1] == [10.0, 11.0, 12.0]
288+
@test ps.discrete[1][2][] == false

‎test/split_parameters.jl

Lines changed: 6 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -4,6 +4,7 @@ using OrdinaryDiffEq
44
using ModelingToolkit: t_nounits as t, D_nounits as D
55
using ModelingToolkit: MTKParameters, ParameterIndex, DEPENDENT_PORTION, NONNUMERIC_PORTION
66
using SciMLStructures: Tunable, Discrete, Constants
7+
using StaticArrays: SizedVector
78

89
x = [1, 2.0, false, [1, 2, 3], Parameter(1.0)]
910

@@ -194,22 +195,22 @@ S = get_sensitivity(closed_loop, :u)
194195

195196
@testset "Indexing MTKParameters with ParameterIndex" begin
196197
ps = MTKParameters(([1.0, 2.0], [3, 4]),
197-
([true, false], [[1 2; 3 4]]),
198+
SizedVector{2}([([true, false], [[1 2; 3 4]]), ([false, true], [[2 4; 6 8]])]),
198199
([5, 6],),
199200
([7.0, 8.0],),
200201
(["hi", "bye"], [:lie, :die]),
201202
nothing,
202203
nothing)
203204
@test ps[ParameterIndex(Tunable(), (1, 2))] === 2.0
204205
@test ps[ParameterIndex(Tunable(), (2, 2))] === 4
205-
@test ps[ParameterIndex(Discrete(), (2, 1, 2, 2))] === 4
206-
@test ps[ParameterIndex(Discrete(), (2, 1))] == [1 2; 3 4]
206+
@test ps[ParameterIndex(Discrete(), (1, 2, 1, 2, 2))] === 4
207+
@test ps[ParameterIndex(Discrete(), (2, 2, 1))] == [2 4; 6 8]
207208
@test ps[ParameterIndex(Constants(), (1, 1))] === 5
208209
@test ps[ParameterIndex(DEPENDENT_PORTION, (1, 1))] === 7.0
209210
@test ps[ParameterIndex(NONNUMERIC_PORTION, (2, 2))] === :die
210211

211212
ps[ParameterIndex(Tunable(), (1, 2))] = 3.0
212-
ps[ParameterIndex(Discrete(), (2, 1, 2, 2))] = 5
213+
ps[ParameterIndex(Discrete(), (1, 2, 1, 2, 2))] = 5
213214
@test ps[ParameterIndex(Tunable(), (1, 2))] === 3.0
214-
@test ps[ParameterIndex(Discrete(), (2, 1, 2, 2))] === 5
215+
@test ps[ParameterIndex(Discrete(), (1, 2, 1, 2, 2))] === 5
215216
end

‎test/symbolic_indexing_interface.jl

Lines changed: 145 additions & 84 deletions
Original file line numberDiff line numberDiff line change
@@ -1,87 +1,148 @@
11
using ModelingToolkit, SymbolicIndexingInterface, SciMLBase
22
using ModelingToolkit: t_nounits as t, D_nounits as D
33

4-
@parameters a b
5-
@variables x(t)=1.0 y(t)=2.0 xy(t)
6-
eqs = [D(x) ~ a * y + t, D(y) ~ b * t]
7-
@named odesys = ODESystem(eqs, t, [x, y], [a, b]; observed = [xy ~ x + y])
8-
9-
@test all(is_variable.((odesys,), [x, y, 1, 2, :x, :y]))
10-
@test all(.!is_variable.((odesys,), [a, b, t, 3, 0, :a, :b]))
11-
@test variable_index.((odesys,), [x, y, a, b, t, 1, 2, :x, :y, :a, :b]) ==
12-
[1, 2, nothing, nothing, nothing, 1, 2, 1, 2, nothing, nothing]
13-
@test isequal(variable_symbols(odesys), [x, y])
14-
@test all(is_parameter.((odesys,), [a, b, 1, 2, :a, :b]))
15-
@test all(.!is_parameter.((odesys,), [x, y, t, 3, 0, :x, :y]))
16-
@test parameter_index.((odesys,), [x, y, a, b, t, 1, 2, :x, :y, :a, :b]) ==
17-
[nothing, nothing, 1, 2, nothing, 1, 2, nothing, nothing, 1, 2]
18-
@test isequal(parameter_symbols(odesys), [a, b])
19-
@test all(is_independent_variable.((odesys,), [t, :t]))
20-
@test all(.!is_independent_variable.((odesys,), [x, y, a, :x, :y, :a]))
21-
@test isequal(independent_variable_symbols(odesys), [t])
22-
@test is_time_dependent(odesys)
23-
@test constant_structure(odesys)
24-
@test !isempty(default_values(odesys))
25-
@test default_values(odesys)[x] == 1.0
26-
@test default_values(odesys)[y] == 2.0
27-
@test isequal(default_values(odesys)[xy], x + y)
28-
29-
@named odesys = ODESystem(
30-
eqs, t, [x, y], [a, b]; defaults = [xy => 3.0], observed = [xy ~ x + y])
31-
@test default_values(odesys)[xy] == 3.0
32-
33-
@variables x y z
34-
@parameters σ ρ β
35-
36-
eqs = [0 ~ σ * (y - x),
37-
0 ~ x * (ρ - z) - y,
38-
0 ~ x * y - β * z]
39-
@named ns = NonlinearSystem(eqs, [x, y, z], [σ, ρ, β])
40-
41-
@test !is_time_dependent(ns)
42-
43-
@parameters x
44-
@variables u(..)
45-
Dxx = Differential(x)^2
46-
Dtt = Differential(t)^2
47-
Dt = D
48-
49-
#2D PDE
50-
C = 1
51-
eq = Dtt(u(t, x)) ~ C^2 * Dxx(u(t, x))
52-
53-
# Initial and boundary conditions
54-
bcs = [u(t, 0) ~ 0.0,# for all t > 0
55-
u(t, 1) ~ 0.0,# for all t > 0
56-
u(0, x) ~ x * (1.0 - x), #for all 0 < x < 1
57-
Dt(u(0, x)) ~ 0.0] #for all 0 < x < 1]
58-
59-
# Space and time domains
60-
domains = [t ∈ (0.0, 1.0),
61-
x ∈ (0.0, 1.0)]
62-
63-
@named pde_system = PDESystem(eq, bcs, domains, [t, x], [u])
64-
65-
@test pde_system.ps == SciMLBase.NullParameters()
66-
@test parameter_symbols(pde_system) == []
67-
68-
@parameters x
69-
@constants h = 1
70-
@variables u(..)
71-
Dt = D
72-
Dxx = Differential(x)^2
73-
eq = Dt(u(t, x)) ~ h * Dxx(u(t, x))
74-
bcs = [u(0, x) ~ -h * x * (x - 1) * sin(x),
75-
u(t, 0) ~ 0, u(t, 1) ~ 0]
76-
77-
domains = [t ∈ (0.0, 1.0),
78-
x ∈ (0.0, 1.0)]
79-
80-
analytic = [u(t, x) ~ -h * x * (x - 1) * sin(x) * exp(-2 * h * t)]
81-
analytic_function = (ps, t, x) -> -ps[1] * x * (x - 1) * sin(x) * exp(-2 * ps[1] * t)
82-
83-
@named pdesys = PDESystem(eq, bcs, domains, [t, x], [u], [h], analytic = analytic)
84-
85-
@test isequal(pdesys.ps, [h])
86-
@test isequal(parameter_symbols(pdesys), [h])
87-
@test isequal(parameters(pdesys), [h])
4+
@testset "ODESystem" begin
5+
@parameters a b
6+
@variables x(t)=1.0 y(t)=2.0 xy(t)
7+
eqs = [D(x) ~ a * y + t, D(y) ~ b * t]
8+
@named odesys = ODESystem(eqs, t, [x, y], [a, b]; observed = [xy ~ x + y])
9+
10+
@test all(is_variable.((odesys,), [x, y, 1, 2, :x, :y]))
11+
@test all(.!is_variable.((odesys,), [a, b, t, 3, 0, :a, :b]))
12+
@test variable_index.((odesys,), [x, y, a, b, t, 1, 2, :x, :y, :a, :b]) ==
13+
[1, 2, nothing, nothing, nothing, 1, 2, 1, 2, nothing, nothing]
14+
@test isequal(variable_symbols(odesys), [x, y])
15+
@test all(is_parameter.((odesys,), [a, b, 1, 2, :a, :b]))
16+
@test all(.!is_parameter.((odesys,), [x, y, t, 3, 0, :x, :y]))
17+
@test parameter_index.((odesys,), [x, y, a, b, t, 1, 2, :x, :y, :a, :b]) ==
18+
[nothing, nothing, 1, 2, nothing, 1, 2, nothing, nothing, 1, 2]
19+
@test isequal(parameter_symbols(odesys), [a, b])
20+
@test all(is_independent_variable.((odesys,), [t, :t]))
21+
@test all(.!is_independent_variable.((odesys,), [x, y, a, :x, :y, :a]))
22+
@test isequal(independent_variable_symbols(odesys), [t])
23+
@test is_time_dependent(odesys)
24+
@test constant_structure(odesys)
25+
@test !isempty(default_values(odesys))
26+
@test default_values(odesys)[x] == 1.0
27+
@test default_values(odesys)[y] == 2.0
28+
@test isequal(default_values(odesys)[xy], x + y)
29+
30+
@named odesys = ODESystem(
31+
eqs, t, [x, y], [a, b]; defaults = [xy => 3.0], observed = [xy ~ x + y])
32+
@test default_values(odesys)[xy] == 3.0
33+
pobs = parameter_observed(odesys, a + b)
34+
@test pobs.timeseries_idx === nothing
35+
@test pobs.observed_fn(
36+
ModelingToolkit.MTKParameters(odesys, [a => 1.0, b => 2.0]), 0.0) ≈ 3.0
37+
pobs = parameter_observed(odesys, [a + b, a - b])
38+
@test pobs.timeseries_idx === nothing
39+
@test pobs.observed_fn(
40+
ModelingToolkit.MTKParameters(odesys, [a => 1.0, b => 2.0]), 0.0) ≈ [3.0, -1.0]
41+
end
42+
43+
@testset "Clock system" begin
44+
dt = 0.1
45+
dt2 = 0.2
46+
@variables x(t)=0 y(t)=0 u(t)=0 yd1(t)=0 ud1(t)=0 yd2(t)=0 ud2(t)=0
47+
@parameters kp=1 r=1
48+
49+
eqs = [
50+
# controller (time discrete part `dt=0.1`)
51+
yd1 ~ Sample(t, dt)(y)
52+
ud1 ~ kp * (r - yd1)
53+
# controller (time discrete part `dt=0.2`)
54+
yd2 ~ Sample(t, dt2)(y)
55+
ud2 ~ kp * (r - yd2)
56+
57+
# plant (time continuous part)
58+
u ~ Hold(ud1) + Hold(ud2)
59+
D(x) ~ -x + u
60+
y ~ x]
61+
62+
@mtkbuild cl = ODESystem(eqs, t)
63+
ts_params = [
64+
Hold(ud1), Hold(ud2), Sample(t, dt)(y), Sample(t, dt2)(y), ud1, ud2, yd1, yd2]
65+
ts_idxs = [1, 2, 1, 2, 1, 2, 1, 2]
66+
@test all(Base.Fix1(is_timeseries_parameter, cl),
67+
[Hold(ud1), Hold(ud2), Sample(t, dt)(y), Sample(t, dt2)(y), ud1, ud2, yd1, yd2])
68+
for (ts_param, ts_idx) in zip(ts_params[1:4], ts_idxs[1:4])
69+
@test timeseries_parameter_index(cl, ts_param).timeseries_idx == ts_idx
70+
end
71+
ps = ModelingToolkit.MTKParameters(cl, [kp => 1.0, Sample(t, dt)(y) => 1.0])
72+
pobs = parameter_observed(cl, Shift(t, 1)(yd1))
73+
@test pobs.timeseries_idx == 1
74+
@test pobs.observed_fn(ps, 0.0) == 1.0
75+
pobs = parameter_observed(cl, [Shift(t, 1)(yd1), Shift(t, 1)(ud1)])
76+
@test pobs.timeseries_idx == 1
77+
@test pobs.observed_fn(ps, 0.0) == [1.0, 0.0]
78+
pobs = parameter_observed(cl, [Shift(t, 1)(yd1), Shift(t, 1)(ud2)])
79+
@test pobs.timeseries_idx === nothing
80+
@test pobs.observed_fn(ps, 0.0) == [1.0, 1.0]
81+
end
82+
83+
@testset "Nonlinear system" begin
84+
@variables x y z
85+
@parameters σ ρ β
86+
87+
eqs = [0 ~ σ * (y - x),
88+
0 ~ x * (ρ - z) - y,
89+
0 ~ x * y - β * z]
90+
@named ns = NonlinearSystem(eqs, [x, y, z], [σ, ρ, β])
91+
ns = complete(ns)
92+
@test !is_time_dependent(ns)
93+
ps = ModelingToolkit.MTKParameters(ns, [σ => 1.0, ρ => 2.0, β => 3.0])
94+
pobs = parameter_observed(ns, σ + ρ)
95+
@test pobs.timeseries_idx === nothing
96+
@test pobs.observed_fn(ps) == 3.0
97+
pobs = parameter_observed(ns, [σ + ρ, ρ + β])
98+
@test pobs.timeseries_idx === nothing
99+
@test pobs.observed_fn(ps) == [3.0, 5.0]
100+
end
101+
102+
@testset "PDESystem" begin
103+
@parameters x
104+
@variables u(..)
105+
Dxx = Differential(x)^2
106+
Dtt = Differential(t)^2
107+
Dt = D
108+
109+
#2D PDE
110+
C = 1
111+
eq = Dtt(u(t, x)) ~ C^2 * Dxx(u(t, x))
112+
113+
# Initial and boundary conditions
114+
bcs = [u(t, 0) ~ 0.0,# for all t > 0
115+
u(t, 1) ~ 0.0,# for all t > 0
116+
u(0, x) ~ x * (1.0 - x), #for all 0 < x < 1
117+
Dt(u(0, x)) ~ 0.0] #for all 0 < x < 1]
118+
119+
# Space and time domains
120+
domains = [t ∈ (0.0, 1.0),
121+
x ∈ (0.0, 1.0)]
122+
123+
@named pde_system = PDESystem(eq, bcs, domains, [t, x], [u])
124+
125+
@test pde_system.ps == SciMLBase.NullParameters()
126+
@test parameter_symbols(pde_system) == []
127+
128+
@parameters x
129+
@constants h = 1
130+
@variables u(..)
131+
Dt = D
132+
Dxx = Differential(x)^2
133+
eq = Dt(u(t, x)) ~ h * Dxx(u(t, x))
134+
bcs = [u(0, x) ~ -h * x * (x - 1) * sin(x),
135+
u(t, 0) ~ 0, u(t, 1) ~ 0]
136+
137+
domains = [t ∈ (0.0, 1.0),
138+
x ∈ (0.0, 1.0)]
139+
140+
analytic = [u(t, x) ~ -h * x * (x - 1) * sin(x) * exp(-2 * h * t)]
141+
analytic_function = (ps, t, x) -> -ps[1] * x * (x - 1) * sin(x) * exp(-2 * ps[1] * t)
142+
143+
@named pdesys = PDESystem(eq, bcs, domains, [t, x], [u], [h], analytic = analytic)
144+
145+
@test isequal(pdesys.ps, [h])
146+
@test isequal(parameter_symbols(pdesys), [h])
147+
@test isequal(parameters(pdesys), [h])
148+
end

0 commit comments

Comments
 (0)