Skip to content

Commit f1344b2

Browse files
AayushSabharwalAayushSabharwal
committed
test: test implementation of SII/SciMLBase discrete saving interface
1 parent ef728cc commit f1344b2

File tree

3 files changed

+190
-90
lines changed

3 files changed

+190
-90
lines changed

‎test/mtkparameters.jl

Lines changed: 26 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -307,3 +307,29 @@ end
307307
newoprob = remake(oprob_scal_scal; p = ps_vec)
308308
@test newoprob.ps[k] == [2.0, 3.0, 4.0, 5.0]
309309
end
310+
311+
# Parameter timeseries
312+
ps = MTKParameters(([1.0, 1.0],), SizedArray{2}([([0.0, 0.0],), ([0.0, 0.0],)]), (), (), (), nothing, nothing)
313+
with_updated_parameter_timeseries_values(
314+
ps, 1 => ModelingToolkit.NestedGetIndex(([5.0, 10.0],)))
315+
@test ps.discrete[1][1] == [5.0, 10.0]
316+
with_updated_parameter_timeseries_values(
317+
ps, 1 => ModelingToolkit.NestedGetIndex(([3.0, 30.0],)),
318+
2 => ModelingToolkit.NestedGetIndex(([4.0, 40.0],)))
319+
@test ps.discrete[1][1] == [3.0, 30.0]
320+
@test ps.discrete[2][1] == [4.0, 40.0]
321+
@test SciMLBase.get_saveable_values(ps, 1).x == ps.discrete[1]
322+
323+
# With multiple types and clocks
324+
ps = MTKParameters((), SizedVector{2}([([1.0, 2.0, 3.0], [false]), ([4.0, 5.0, 6.0], Bool[])]), (), (), (), nothing, nothing)
325+
@test SciMLBase.get_saveable_values(ps, 1).x isa Tuple{Vector{Float64}, Vector{Bool}}
326+
tsidx1 = 1
327+
tsidx2 = 2
328+
@test length(ps.discrete[tsidx1][1]) == 3
329+
@test length(ps.discrete[tsidx1][2]) == 1
330+
@test length(ps.discrete[tsidx2][1]) == 3
331+
@test length(ps.discrete[tsidx2][2]) == 0
332+
with_updated_parameter_timeseries_values(
333+
ps, tsidx1 => ModelingToolkit.NestedGetIndex(([10.0, 11.0, 12.0], [false])))
334+
@test ps.discrete[tsidx1][1] == [10.0, 11.0, 12.0]
335+
@test ps.discrete[tsidx1][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: 158 additions & 85 deletions
Original file line numberDiff line numberDiff line change
@@ -1,90 +1,155 @@
11
using ModelingToolkit, SymbolicIndexingInterface, SciMLBase
2-
using ModelingToolkit: t_nounits as t, D_nounits as D
2+
using ModelingToolkit: t_nounits as t, D_nounits as D, ParameterIndex
3+
using SciMLStructures: Tunable
4+
5+
@testset "ODESystem" begin
6+
@parameters a b
7+
@variables x(t)=1.0 y(t)=2.0 xy(t)
8+
eqs = [D(x) ~ a * y + t, D(y) ~ b * t]
9+
@named odesys = ODESystem(eqs, t, [x, y], [a, b]; observed = [xy ~ x + y])
10+
odesys = complete(odesys)
11+
@test all(is_variable.((odesys,), [x, y, 1, 2, :x, :y]))
12+
@test all(.!is_variable.((odesys,), [a, b, t, 3, 0, :a, :b]))
13+
@test variable_index.((odesys,), [x, y, a, b, t, 1, 2, :x, :y, :a, :b]) ==
14+
[1, 2, nothing, nothing, nothing, 1, 2, 1, 2, nothing, nothing]
15+
@test isequal(variable_symbols(odesys), [x, y])
16+
@test all(is_parameter.((odesys,), [a, b, ParameterIndex(Tunable(), (1, 1)), :a, :b]))
17+
@test all(.!is_parameter.((odesys,), [x, y, t, 3, 0, :x, :y]))
18+
@test parameter_index(odesys, a) == parameter_index(odesys, :a)
19+
@test parameter_index(odesys, a) isa ParameterIndex{Tunable, Tuple{Int, Int}}
20+
@test parameter_index(odesys, b) == parameter_index(odesys, :b)
21+
@test parameter_index(odesys, b) isa ParameterIndex{Tunable, Tuple{Int, Int}}
22+
@test parameter_index.((odesys,), [x, y, t, ParameterIndex(Tunable(), (1, 1)), :x, :y,]) ==
23+
[nothing, nothing, nothing, ParameterIndex(Tunable(), (1, 1)), nothing, nothing]
24+
@test isequal(parameter_symbols(odesys), [a, b])
25+
@test all(is_independent_variable.((odesys,), [t, :t]))
26+
@test all(.!is_independent_variable.((odesys,), [x, y, a, :x, :y, :a]))
27+
@test isequal(independent_variable_symbols(odesys), [t])
28+
@test is_time_dependent(odesys)
29+
@test constant_structure(odesys)
30+
@test !isempty(default_values(odesys))
31+
@test default_values(odesys)[x] == 1.0
32+
@test default_values(odesys)[y] == 2.0
33+
@test isequal(default_values(odesys)[xy], x + y)
34+
35+
@named odesys = ODESystem(
36+
eqs, t, [x, y], [a, b]; defaults = [xy => 3.0], observed = [xy ~ x + y])
37+
odesys = complete(odesys)
38+
@test default_values(odesys)[xy] == 3.0
39+
pobs = parameter_observed(odesys, a + b)
40+
@test pobs.timeseries_idx === nothing
41+
@test pobs.observed_fn(
42+
ModelingToolkit.MTKParameters(odesys, [a => 1.0, b => 2.0]), 0.0) ≈ 3.0
43+
pobs = parameter_observed(odesys, [a + b, a - b])
44+
@test pobs.timeseries_idx === nothing
45+
@test pobs.observed_fn(
46+
ModelingToolkit.MTKParameters(odesys, [a => 1.0, b => 2.0]), 0.0) ≈ [3.0, -1.0]
47+
end
348

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])
49+
# @testset "Clock system" begin
50+
# dt = 0.1
51+
# dt2 = 0.2
52+
# @variables x(t)=0 y(t)=0 u(t)=0 yd1(t)=0 ud1(t)=0 yd2(t)=0 ud2(t)=0
53+
# @parameters kp=1 r=1
54+
55+
# eqs = [
56+
# # controller (time discrete part `dt=0.1`)
57+
# yd1 ~ Sample(t, dt)(y)
58+
# ud1 ~ kp * (r - yd1)
59+
# # controller (time discrete part `dt=0.2`)
60+
# yd2 ~ Sample(t, dt2)(y)
61+
# ud2 ~ kp * (r - yd2)
62+
63+
# # plant (time continuous part)
64+
# u ~ Hold(ud1) + Hold(ud2)
65+
# D(x) ~ -x + u
66+
# y ~ x]
67+
68+
# @mtkbuild cl = ODESystem(eqs, t)
69+
# partition1_params = [Hold(ud1), Sample(t, dt)(y), ud1, yd1]
70+
# partition2_params = [Hold(ud2), Sample(t, dt2)(y), ud2, yd2]
71+
# @test all(
72+
# Base.Fix1(is_timeseries_parameter, cl), vcat(partition1_params, partition2_params))
73+
# @test allequal(timeseries_parameter_index(cl, p).timeseries_idx
74+
# for p in partition1_params)
75+
# @test allequal(timeseries_parameter_index(cl, p).timeseries_idx
76+
# for p in partition2_params)
77+
# tsidx1 = timeseries_parameter_index(cl, partition1_params[1]).timeseries_idx
78+
# tsidx2 = timeseries_parameter_index(cl, partition2_params[1]).timeseries_idx
79+
# @test tsidx1 != tsidx2
80+
# ps = ModelingToolkit.MTKParameters(cl, [kp => 1.0, Sample(t, dt)(y) => 1.0])
81+
# pobs = parameter_observed(cl, Shift(t, 1)(yd1))
82+
# @test pobs.timeseries_idx == tsidx1
83+
# @test pobs.observed_fn(ps, 0.0) == 1.0
84+
# pobs = parameter_observed(cl, [Shift(t, 1)(yd1), Shift(t, 1)(ud1)])
85+
# @test pobs.timeseries_idx == tsidx1
86+
# @test pobs.observed_fn(ps, 0.0) == [1.0, 0.0]
87+
# pobs = parameter_observed(cl, [Shift(t, 1)(yd1), Shift(t, 1)(ud2)])
88+
# @test pobs.timeseries_idx === nothing
89+
# @test pobs.observed_fn(ps, 0.0) == [1.0, 1.0]
90+
# end
91+
92+
@testset "Nonlinear system" begin
93+
@variables x y z
94+
@parameters σ ρ β
95+
96+
eqs = [0 ~ σ * (y - x),
97+
0 ~ x * (ρ - z) - y,
98+
0 ~ x * y - β * z]
99+
@named ns = NonlinearSystem(eqs, [x, y, z], [σ, ρ, β])
100+
ns = complete(ns)
101+
@test !is_time_dependent(ns)
102+
ps = ModelingToolkit.MTKParameters(ns, [σ => 1.0, ρ => 2.0, β => 3.0])
103+
pobs = parameter_observed(ns, σ + ρ)
104+
@test pobs.timeseries_idx === nothing
105+
@test pobs.observed_fn(ps) == 3.0
106+
pobs = parameter_observed(ns, [σ + ρ, ρ + β])
107+
@test pobs.timeseries_idx === nothing
108+
@test pobs.observed_fn(ps) == [3.0, 5.0]
109+
end
110+
111+
@testset "PDESystem" begin
112+
@parameters x
113+
@variables u(..)
114+
Dxx = Differential(x)^2
115+
Dtt = Differential(t)^2
116+
Dt = D
117+
118+
#2D PDE
119+
C = 1
120+
eq = Dtt(u(t, x)) ~ C^2 * Dxx(u(t, x))
121+
122+
# Initial and boundary conditions
123+
bcs = [u(t, 0) ~ 0.0,# for all t > 0
124+
u(t, 1) ~ 0.0,# for all t > 0
125+
u(0, x) ~ x * (1.0 - x), #for all 0 < x < 1
126+
Dt(u(0, x)) ~ 0.0] #for all 0 < x < 1]
127+
128+
# Space and time domains
129+
domains = [t ∈ (0.0, 1.0),
130+
x ∈ (0.0, 1.0)]
131+
132+
@named pde_system = PDESystem(eq, bcs, domains, [t, x], [u])
133+
134+
@test pde_system.ps == SciMLBase.NullParameters()
135+
@test parameter_symbols(pde_system) == []
136+
137+
@parameters x
138+
@constants h = 1
139+
@variables u(..)
140+
Dt = D
141+
Dxx = Differential(x)^2
142+
eq = Dt(u(t, x)) ~ h * Dxx(u(t, x))
143+
bcs = [u(0, x) ~ -h * x * (x - 1) * sin(x),
144+
u(t, 0) ~ 0, u(t, 1) ~ 0]
145+
146+
domains = [t ∈ (0.0, 1.0),
147+
x ∈ (0.0, 1.0)]
148+
149+
@test isequal(pdesys.ps, [h])
150+
@test isequal(parameter_symbols(pdesys), [h])
151+
@test isequal(parameters(pdesys), [h])
152+
end
88153

89154
# Issue#2767
90155
using ModelingToolkit
@@ -113,4 +178,12 @@ get_dep = @test_nowarn getu(prob, 2p1)
113178
@test getu(prob, z)(prob) == getu(prob, :z)(prob)
114179
@test getu(prob, p1)(prob) == getu(prob, :p1)(prob)
115180
@test getu(prob, p2)(prob) == getu(prob, :p2)(prob)
181+
analytic = [u(t, x) ~ -h * x * (x - 1) * sin(x) * exp(-2 * h * t)]
182+
analytic_function = (ps, t, x) -> -ps[1] * x * (x - 1) * sin(x) * exp(-2 * ps[1] * t)
183+
184+
@named pdesys = PDESystem(eq, bcs, domains, [t, x], [u], [h], analytic = analytic)
185+
186+
@test isequal(pdesys.ps, [h])
187+
@test isequal(parameter_symbols(pdesys), [h])
188+
@test isequal(parameters(pdesys), [h])
116189
end

0 commit comments

Comments
 (0)