test: account for new @constants behavior in tests

Author: AayushSabharwal

test/components.jl

@@ -230,7 +230,7 @@ end
         eqs = [
             v ~ i * R
         ]
-        extend(System(eqs, t, [], []; name = name), oneport)
+        extend(System(eqs, t, [], [R]; name = name), oneport)
     end
     capacitor = Capacitor(; name = :c1, C = 1.0)
     resistor = FixedResistor(; name = :r1)

test/constants.jl

@@ -4,7 +4,8 @@ MT = ModelingToolkit
 UMT = ModelingToolkit.UnitfulUnitCheck
 
 @constants a = 1
-@test_throws MT.ArgumentError @constants b
+@test isconstant(a)
+@test !istunable(a)
 
 @independent_variables t
 @variables x(t) w(t)
@@ -14,9 +15,6 @@ eqs = [D(x) ~ a]
 prob = ODEProblem(complete(sys), [0], [0.0, 1.0], [])
 sol = solve(prob, Tsit5())
 
-newsys = MT.eliminate_constants(sys)
-@test isequal(equations(newsys), [D(x) ~ 1])
-
 # Test structural_simplify substitutions & observed values
 eqs = [D(x) ~ 1,
     w ~ a]
@@ -29,24 +27,23 @@ simp = structural_simplify(sys)
 @constants β=1 [unit = u"m/s"]
 UMT.get_unit(β)
 @test MT.isconstant(β)
+@test !MT.istunable(β)
 @independent_variables t [unit = u"s"]
 @variables x(t) [unit = u"m"]
 D = Differential(t)
 eqs = [D(x) ~ β]
 @named sys = System(eqs, t)
 simp = structural_simplify(sys)
 
-@test isempty(MT.collect_constants(nothing))
-
 @testset "Issue#3044" begin
-    @constants h = 1
+    @constants h
     @parameters τ = 0.5 * h
     @variables x(MT.t_nounits) = h
     eqs = [MT.D_nounits(x) ~ (h - x) / τ]
 
     @mtkbuild fol_model = System(eqs, MT.t_nounits)
 
-    prob = ODEProblem(fol_model, [], (0.0, 10.0))
+    prob = ODEProblem(fol_model, [], (0.0, 10.0), [h => 1])
     @test prob[x] ≈ 1
     @test prob.ps[τ] ≈ 0.5
 end

test/discrete_system.jl

@@ -30,7 +30,7 @@ eqs = [S ~ S(k - 1) - infection * h,
     R ~ R(k - 1) + recovery]
 
 # System
-@named sys = System(eqs, t, [S, I, R], [c, nsteps, δt, β, γ])
+@named sys = System(eqs, t, [S, I, R], [c, nsteps, δt, β, γ, h])
 syss = structural_simplify(sys)
 @test syss == syss
 

test/dq_units.jl

@@ -233,8 +233,7 @@ end
     L(t), [unit = u"m"]
     L_out(t), [unit = u"1"]
 end
-@test to_m in ModelingToolkit.vars(ModelingToolkit.fold_constants(Symbolics.unwrap(L_out *
-                                                                                   -to_m)))
+@test to_m in ModelingToolkit.vars(Symbolics.unwrap(L_out * -to_m))
 
 # test units for registered functions
 let

test/funcaffect.jl

@@ -282,7 +282,7 @@ function bb_affect!(integ, u, p, ctx)
     integ.u[u.v] = -integ.u[u.v]
 end
 
-@named bb_model = System(bb_eqs, t, sts, par,
+@named bb_model = System(bb_eqs, t, sts, [par; zr],
     continuous_events = [
         [y ~ zr] => (bb_affect!, [v], [], [], nothing)
     ])

test/input_output_handling.jl

@@ -447,7 +447,7 @@ end
     @constants c = 2.0
     @variables x(t)
     eqs = [D(x) ~ c * x]
-    @mtkbuild sys = System(eqs, t, [x], [])
+    @mtkbuild sys = System(eqs, t, [x], [c])
 
     f, dvs, ps, io_sys = ModelingToolkit.generate_control_function(sys)
     @test f[1]([0.5], nothing, MTKParameters(io_sys, []), 0.0) ≈ [1.0]

test/jumpsystem.jl

@@ -1,4 +1,5 @@
 using ModelingToolkit, DiffEqBase, JumpProcesses, Test, LinearAlgebra
+using SymbolicIndexingInterface
 using Random, StableRNGs, NonlinearSolve
 using OrdinaryDiffEq
 using ModelingToolkit: t_nounits as t, D_nounits as D
@@ -17,7 +18,7 @@ rate₂ = γ * I + t
 affect₂ = [I ~ Pre(I) - 1, R ~ Pre(R) + 1]
 j₁ = ConstantRateJump(rate₁, affect₁)
 j₂ = VariableRateJump(rate₂, affect₂)
-@named js = JumpSystem([j₁, j₂], t, [S, I, R], [β, γ])
+@named js = JumpSystem([j₁, j₂], t, [S, I, R], [β, γ, h])
 unknowntoid = Dict(MT.value(unknown) => i for (i, unknown) in enumerate(unknowns(js)))
 mtjump1 = MT.assemble_crj(js, j₁, unknowntoid)
 mtjump2 = MT.assemble_vrj(js, j₂, unknowntoid)
@@ -38,7 +39,7 @@ jump2 = VariableRateJump(rate2, affect2!)
 
 # test crjs
 u = [100, 9, 5]
-p = (0.1 / 1000, 0.01)
+p = (0.1 / 1000, 0.01, 1)
 tf = 1.0
 mutable struct TestInt{U, V, T}
     u::U
@@ -62,15 +63,15 @@ jump2.affect!(integrator)
 rate₃ = γ * I * h
 affect₃ = [I ~ Pre(I) * h - 1, R ~ Pre(R) + 1]
 j₃ = ConstantRateJump(rate₃, affect₃)
-@named js2 = JumpSystem([j₁, j₃], t, [S, I, R], [β, γ])
+@named js2 = JumpSystem([j₁, j₃], t, [S, I, R], [β, γ, h])
 js2 = complete(js2)
 u₀ = [999, 1, 0];
-p = (0.1 / 1000, 0.01);
 tspan = (0.0, 250.0);
 u₀map = [S => 999, I => 1, R => 0]
 parammap = [β => 0.1 / 1000, γ => 0.01]
 jprob = JumpProblem(js2, u₀map, tspan, parammap; aggregator = Direct(),
     save_positions = (false, false), rng)
+p = parameter_values(jprob)
 @test jprob.prob isa DiscreteProblem
 Nsims = 30000
 function getmean(jprob, Nsims; use_stepper = true)
@@ -90,7 +91,7 @@ mb = getmean(jprobb, Nsims; use_stepper = false)
 
 @variables S2(t)
 obs = [S2 ~ 2 * S]
-@named js2b = JumpSystem([j₁, j₃], t, [S, I, R], [β, γ], observed = obs)
+@named js2b = JumpSystem([j₁, j₃], t, [S, I, R], [β, γ, h], observed = obs)
 js2b = complete(js2b)
 jprob = JumpProblem(js2b, u₀map, tspan, parammap; aggregator = Direct(),
     save_positions = (false, false), rng)
@@ -110,6 +111,8 @@ jump2 = ConstantRateJump(rate2, affect2!)
 mtjumps = jprob.discrete_jump_aggregation
 @test abs(mtjumps.rates[1](u, p, tf) - jump1.rate(u, p, tf)) < 10 * eps()
 @test abs(mtjumps.rates[2](u, p, tf) - jump2.rate(u, p, tf)) < 10 * eps()
+
+ModelingToolkit.@set! mtintegrator.p = (mtintegrator.p, (1,))
 mtjumps.affects![1](mtintegrator)
 jump1.affect!(integrator)
 @test all(integrator.u .== mtintegrator.u)

test/model_parsing.jl

@@ -511,7 +511,7 @@ using ModelingToolkit: getdefault, scalarize
 
     @test eval(ModelWithComponentArray.structure[:parameters][:r][:unit]) ==
           eval(u"Ω")
-    @test lastindex(parameters(model_with_component_array)) == 3
+    @test lastindex(parameters(model_with_component_array)) == 4
 
     # Test the constant `k`. Manually k's value should be kept in sync here
     # and the ModelParsingPrecompile.

test/nonlinearsystem.jl

@@ -26,13 +26,13 @@ end
 eqs = [0 ~ σ * (y - x) * h,
     0 ~ x * (ρ - z) - y,
     0 ~ x * y - β * z]
-@named ns = System(eqs, [x, y, z], [σ, ρ, β], defaults = Dict(x => 2))
+@named ns = System(eqs, [x, y, z], [σ, ρ, β, h], defaults = Dict(x => 2))
 @test eval(toexpr(ns)) == ns
-test_nlsys_inference("standard", ns, (x, y, z), (σ, ρ, β))
+test_nlsys_inference("standard", ns, (x, y, z), (σ, ρ, β, h))
 @test begin
-    f = generate_rhs(ns, [x, y, z], [σ, ρ, β], expression = Val{false})[2]
+    f = generate_rhs(ns, [x, y, z], [σ, ρ, β, h], expression = Val{false})[2]
     du = [0.0, 0.0, 0.0]
-    f(du, [1, 2, 3], [1, 2, 3])
+    f(du, [1, 2, 3], [1, 2, 3, 1])
     du ≈ [1, -3, -7]
 end
 
@@ -64,9 +64,9 @@ a = y - x
 eqs = [0 ~ σ * a * h,
     0 ~ x * (ρ - z) - y,
     0 ~ x * y - β * z]
-@named ns = System(eqs, [x, y, z], [σ, ρ, β])
+@named ns = System(eqs, [x, y, z], [σ, ρ, β, h])
 ns = complete(ns)
-nlsys_func = generate_rhs(ns, [x, y, z], [σ, ρ, β])
+nlsys_func = generate_rhs(ns, [x, y, z], [σ, ρ, β, h])
 nf = NonlinearFunction(ns)
 jac = calculate_jacobian(ns)
 
@@ -99,7 +99,7 @@ eqs1 = [
     0 ~ x + y - z - u
 ]
 
-lorenz = name -> System(eqs1, [x, y, z, u, F], [σ, ρ, β], name = name)
+lorenz = name -> System(eqs1, [x, y, z, u, F], [σ, ρ, β, h], name = name)
 lorenz1 = lorenz(:lorenz1)
 @test_throws ArgumentError NonlinearProblem(complete(lorenz1), zeros(5), zeros(3))
 lorenz2 = lorenz(:lorenz2)
@@ -132,7 +132,7 @@ sol = solve(prob, FBDF(), reltol = 1e-7, abstol = 1e-7)
 eqs = [0 ~ σ * (y - x),
     0 ~ x * (ρ - z) - y,
     0 ~ x * y - β * z * h]
-@named ns = System(eqs, [x, y, z], [σ, ρ, β])
+@named ns = System(eqs, [x, y, z], [σ, ρ, β, h])
 np = NonlinearProblem(
     complete(ns), [0, 0, 0], [σ => 1, ρ => 2, β => 3], jac = true, sparse = true)
 @test calculate_jacobian(ns, sparse = true) isa SparseMatrixCSC
@@ -214,7 +214,7 @@ testdict = Dict([:test => 1])
     eqs = [0 ~ a * (y - x) * h,
         0 ~ x * (b - z) - y,
         0 ~ x * y - c * z]
-    @named sys = System(eqs, [x, y, z], [a, b, c], defaults = Dict(x => 2.0))
+    @named sys = System(eqs, [x, y, z], [a, b, c, h], defaults = Dict(x => 2.0))
     sys = complete(sys)
     prob = NonlinearProblem(sys, ones(length(unknowns(sys))))
 

test/odesystem.jl

@@ -27,7 +27,7 @@ ModelingToolkit.toexpr.(eqs)[1]
 @named de = System(eqs, t; defaults = Dict(x => 1))
 subed = substitute(de, [σ => k])
 ssort(eqs) = sort(eqs, by = string)
-@test isequal(ssort(parameters(subed)), [k, β, ρ])
+@test isequal(ssort(parameters(subed)), [k, β, κ, ρ])
 @test isequal(equations(subed),
     [D(x) ~ k * (y - x)
      D(y) ~ (ρ - z) * x - y
@@ -47,7 +47,7 @@ function test_diffeq_inference(name, sys, iv, dvs, ps)
     end
 end
 
-test_diffeq_inference("standard", de, t, [x, y, z], [ρ, σ, β])
+test_diffeq_inference("standard", de, t, [x, y, z], [ρ, σ, β, κ])
 jac_expr = generate_jacobian(de)
 jac = calculate_jacobian(de)
 jacfun = eval(jac_expr[2])
@@ -138,11 +138,11 @@ tgrad_iip(du, u, p, t)
 eqs = [D(x) ~ σ(t - 1) * (y - x),
     D(y) ~ x * (ρ - z) - y,
     D(z) ~ x * y - β * z * κ]
-@named de = System(eqs, t, [x, y, z], [σ, ρ, β])
-test_diffeq_inference("single internal iv-varying", de, t, (x, y, z), (σ, ρ, β))
+@named de = System(eqs, t, [x, y, z], [σ, ρ, β, κ])
+test_diffeq_inference("single internal iv-varying", de, t, (x, y, z), (σ, ρ, β, κ))
 f = generate_rhs(de, expression = Val{false}, wrap_gfw = Val{true})
 du = [0.0, 0.0, 0.0]
-f(du, [1.0, 2.0, 3.0], [x -> x + 7, 2, 3], 5.0)
+f(du, [1.0, 2.0, 3.0], [x -> x + 7, 2, 3, 1], 5.0)
 @test du ≈ [11, -3, -7]
 
 eqs = [D(x) ~ x + 10σ(t - 1) + 100σ(t - 2) + 1000σ(t^2)]
@@ -1254,7 +1254,8 @@ end
     prob = ODEProblem(sys, u0, (0.0, 1.0), p)
 
     # evaluate
-    u0_v, p_v, _ = ModelingToolkit.get_u0_p(sys, u0, p)
+    u0_v = prob.u0
+    p_v = prob.p
     @test prob.f(u0_v, p_v, 0.0) == [c_b, c_a]
 end
 

test/parameter_dependencies.jl

@@ -296,9 +296,9 @@ end
     j₁ = ConstantRateJump(rate₁, affect₁)
     j₃ = ConstantRateJump(rate₃, affect₃)
     @named js2 = JumpSystem(
-        [j₃], t, [S, I, R], [γ]; parameter_dependencies = [β => 0.01γ])
-    @test isequal(only(parameters(js2)), γ)
-    @test Set(full_parameters(js2)) == Set([γ, β])
+        [j₃], t, [S, I, R], [γ, h]; parameter_dependencies = [β => 0.01γ])
+    @test issetequal(parameters(js2), [γ, h])
+    @test Set(full_parameters(js2)) == Set([γ, β, h])
     js2 = complete(js2)
     tspan = (0.0, 250.0)
     u₀map = [S => 999, I => 1, R => 0]
@@ -310,7 +310,7 @@ end
     @test_nowarn solve(jprob, SSAStepper())
 
     @named js2 = JumpSystem(
-        [j₁, j₃], t, [S, I, R], [γ]; parameter_dependencies = [β => 0.01γ],
+        [j₁, j₃], t, [S, I, R], [γ, h]; parameter_dependencies = [β => 0.01γ],
         discrete_events = [SymbolicDiscreteCallback(
             [10.0] => [γ ~ 0.02], discrete_parameters = [γ])])
     js2 = complete(js2)

test/structural_transformation/tearing.jl

@@ -18,7 +18,7 @@ eqs = [
     0 ~ u4 - hypot(u2, u3),
     0 ~ u5 - hypot(u4, u1)
 ]
-@named sys = System(eqs, [u1, u2, u3, u4, u5], [])
+@named sys = System(eqs, [u1, u2, u3, u4, u5], [h])
 state = TearingState(sys)
 StructuralTransformations.find_solvables!(state)
 
@@ -149,19 +149,17 @@ eqs = [D(x) ~ z * h
        0 ~ sin(z) + y - p * t]
 @named daesys = System(eqs, t)
 newdaesys = structural_simplify(daesys)
-@test equations(newdaesys) == [D(x) ~ z; 0 ~ y + sin(z) - p * t]
-@test equations(tearing_substitution(newdaesys)) == [D(x) ~ z; 0 ~ x + sin(z) - p * t]
+@test equations(newdaesys) == [D(x) ~ h * z; 0 ~ y + sin(z) - p * t]
+@test equations(tearing_substitution(newdaesys)) == [D(x) ~ h * z; 0 ~ x + sin(z) - p * t]
 @test isequal(unknowns(newdaesys), [x, z])
-@test isequal(unknowns(newdaesys), [x, z])
-@test_deprecated ODAEProblem(newdaesys, [x => 1.0, z => -0.5π], (0, 1.0), [p => 0.2])
 prob = ODEProblem(newdaesys, [x => 1.0, z => -0.5π], (0, 1.0), [p => 0.2])
 du = [0.0, 0.0];
 u = [1.0, -0.5π];
-pr = 0.2;
+pr = prob.p;
 tt = 0.1;
 @test (@ballocated $(prob.f)($du, $u, $pr, $tt)) == 0
 prob.f(du, u, pr, tt)
-@test du≈[u[2], u[1] + sin(u[2]) - pr * tt] atol=1e-5
+@test du≈[u[2], u[1] + sin(u[2]) - prob.ps[p] * tt] atol=1e-5
 
 # test the initial guess is respected
 @named sys = System(eqs, t, defaults = Dict(z => NaN))