refactor tests

Author: vyudu

src/systems/diffeqs/abstractodesystem.jl

@@ -943,9 +943,10 @@ function SciMLBase.BVProblem{iip, specialize}(sys::AbstractODESystem, u0map = []
 
     sts = unknowns(sys)
     ps = parameters(sys)
+    constraintsys = get_constraintsystem(sys)
 
-    if !isnothing(constraints)
-        (length(constraints) + length(u0map) > length(sts)) && 
+    if !isnothing(constraintsys)
+        (length(constraints(constraintsys)) + length(u0map) > length(sts)) && 
         @warn "The BVProblem is overdetermined. The total number of conditions (# constraints + # fixed initial values given by u0map) exceeds the total number of states. The BVP solvers will default to doing a nonlinear least-squares optimization."
     end
 
@@ -976,33 +977,35 @@ function generate_function_bc(sys::ODESystem, u0, u0_idxs, tspan, iip)
     ps = get_ps(sys)
     np = length(ps)
     ns = length(sts)
-    conssys = get_constraintsystem(sys)
-    cons = constraints(conssys)
-
     stidxmap = Dict([v => i for (i, v) in enumerate(sts)])
     pidxmap = Dict([v => i for (i, v) in enumerate(ps)])
 
     @variables sol(..)[1:ns] p[1:np]
-    exprs = Any[]
 
-    for st in get_unknowns(cons)
-        x = operation(st)
-        t = first(arguments(st))
-        idx = stidxmap[x(iv)]
+    conssys = get_constraintsystem(sys)
+    cons = Any[]
+    if !isnothing(conssys)
+        cons = [con.lhs - con.rhs for con in constraints(conssys)]
 
-        cons = Symbolics.substitute(cons, Dict(x(t) => sol(t)[idx]))
-    end
+        for st in get_unknowns(conssys)
+            x = operation(st)
+            t = only(arguments(st))
+            idx = stidxmap[x(iv)]
 
-    for var in get_parameters(cons) 
-        if iscall(var)
-            x = operation(var)
-            t = arguments(var)[1]
-            idx = pidxmap[x]
+            cons = map(c -> Symbolics.substitute(c, Dict(x(t) => sol(t)[idx])), cons)
+        end
 
-            cons = Symbolics.substitute(cons, Dict(x(t) => p[idx]))
-        else
-            idx = pidxmap[var]
-            cons = Symbolics.substitute(cons, Dict(var => p[idx]))
+        for var in parameters(conssys) 
+            if iscall(var)
+                x = operation(var)
+                t = only(arguments(var))
+                idx = pidxmap[x]
+
+                cons = map(c -> Symbolics.substitute(c, Dict(x(t) => p[idx])), cons)
+            else
+                idx = pidxmap[var]
+                cons = map(c -> Symbolics.substitute(c, Dict(var => p[idx])), cons)
+            end
         end
     end
 

src/systems/diffeqs/odesystem.jl

@@ -802,13 +802,11 @@ function validate_constraint_syms(constraintsts, constraintps, sts, ps, iv)
             throw(ArgumentError("Too many arguments for variable $var."))
         elseif length(arguments(var)) == 1
             arg = first(arguments(var))
-            operation(var)(iv) ∈ sts || throw(ArgumentError("Variable $var is not a variable of the ODESystem. Called variables must be variables of the ODESystem."))
+            operation(var)(iv) ∈ sts || 
+                throw(ArgumentError("Variable $var is not a variable of the ODESystem. Called variables must be variables of the ODESystem."))
 
-            isequal(arg, iv) || 
-                isparameter(arg) || 
-                arg isa Integer || 
-                    arg isa AbstractFloat || 
-                        throw(ArgumentError("Invalid argument specified for variable $var. The argument of the variable should be either $iv, a parameter, or a value specifying the time that the constraint holds."))
+            isequal(arg, iv) || isparameter(arg) || arg isa Integer || arg isa AbstractFloat || 
+                throw(ArgumentError("Invalid argument specified for variable $var. The argument of the variable should be either $iv, a parameter, or a value specifying the time that the constraint holds."))
         else
             var ∈ sts && @warn "Variable $var has no argument. It will be interpreted as $var($iv), and the constraint will apply to the entire interval."
         end
@@ -824,7 +822,6 @@ function validate_constraint_syms(constraintsts, constraintps, sts, ps, iv)
             operation(var) ∈ ps || throw(ArgumentError("Parameter $var is not a parameter of the ODESystem. Called parameters must be parameters of the ODESystem."))
 
             isequal(arg, iv) || 
-                isparameter(arg) ||
                 arg isa Integer || 
                     arg isa AbstractFloat || 
                         throw(ArgumentError("Invalid argument specified for callable parameter $var. The argument of the parameter should be either $iv, a parameter, or a value specifying the time that the constraint holds."))

test/bvproblem.jl

@@ -5,7 +5,6 @@ using BenchmarkTools
 using ModelingToolkit
 using SciMLBase
 using ModelingToolkit: t_nounits as t, D_nounits as D
-import ModelingToolkit: process_constraints
 
 ### Test Collocation solvers on simple problems 
 solvers = [MIRK4]
@@ -113,8 +112,8 @@ let
     end
 
     u0 = [1., 2.]; p = [1.5, 1., 1., 3.]
-    genbc_iip = ModelingToolkit.generate_function_bc(lksys, nothing, u0, [1, 2], tspan, true)
-    genbc_oop = ModelingToolkit.generate_function_bc(lksys, nothing, u0, [1, 2], tspan, false)
+    genbc_iip = ModelingToolkit.generate_function_bc(lksys, u0, [1, 2], tspan, true)
+    genbc_oop = ModelingToolkit.generate_function_bc(lksys, u0, [1, 2], tspan, false)
 
     bvpi1 = SciMLBase.BVProblem(lotkavolterra!, bc!, [1.,2.], tspan, p)
     bvpi2 = SciMLBase.BVProblem(lotkavolterra!, genbc_iip, [1.,2.], tspan, p)
@@ -131,7 +130,8 @@ let
     @test sol1 ≈ sol2
 
     # Test with a constraint.
-    constraints = [y(0.5) ~ 2.]
+    constr = [y(0.5) ~ 2.]
+    @mtkbuild lksys = ODESystem(eqs, t; constraints = constr)
 
     function bc!(resid, u, p, t) 
         resid[1] = u(0.0)[1] - 1.
@@ -142,13 +142,13 @@ let
     end
 
     u0 = [1, 1.]
-    genbc_iip = ModelingToolkit.generate_function_bc(lksys, constraints, u0, [1], tspan, true)
-    genbc_oop = ModelingToolkit.generate_function_bc(lksys, constraints, u0, [1], tspan, false)
+    genbc_iip = ModelingToolkit.generate_function_bc(lksys, u0, [1], tspan, true)
+    genbc_oop = ModelingToolkit.generate_function_bc(lksys, u0, [1], tspan, false)
 
     bvpi1 = SciMLBase.BVProblem(lotkavolterra!, bc!, u0, tspan, p)
     bvpi2 = SciMLBase.BVProblem(lotkavolterra!, genbc_iip, u0, tspan, p)
-    bvpi3 = SciMLBase.BVProblem{true, SciMLBase.AutoSpecialize}(lksys, [x(t) => 1.], tspan; guesses = [y(t) => 1.], constraints)
-    bvpi4 = SciMLBase.BVProblem{true, SciMLBase.FullSpecialize}(lksys, [x(t) => 1.], tspan; guesses = [y(t) => 1.], constraints)
+    bvpi3 = SciMLBase.BVProblem{true, SciMLBase.AutoSpecialize}(lksys, [x(t) => 1.], tspan; guesses = [y(t) => 1.])
+    bvpi4 = SciMLBase.BVProblem{true, SciMLBase.FullSpecialize}(lksys, [x(t) => 1.], tspan; guesses = [y(t) => 1.])
 
     sol1 = @btime solve($bvpi1, MIRK4(), dt = 0.01)
     sol2 = @btime solve($bvpi2, MIRK4(), dt = 0.01)
@@ -158,7 +158,7 @@ let
 
     bvpo1 = BVProblem(lotkavolterra, bc, u0, tspan, p)
     bvpo2 = BVProblem(lotkavolterra, genbc_oop, u0, tspan, p)
-    bvpo3 = SciMLBase.BVProblem{false, SciMLBase.FullSpecialize}(lksys, [x(t) => 1.], tspan; guesses = [y(t) => 1.], constraints)
+    bvpo3 = SciMLBase.BVProblem{false, SciMLBase.FullSpecialize}(lksys, [x(t) => 1.], tspan; guesses = [y(t) => 1.])
 
     sol1 = @btime solve($bvpo1, MIRK4(), dt = 0.05)
     sol2 = @btime solve($bvpo2, MIRK4(), dt = 0.05)
@@ -197,12 +197,6 @@ function test_solvers(solvers, prob, u0map, constraints, equations = []; dt = 0.
     end
 end
 
-solvers = [RadauIIa3, RadauIIa5, RadauIIa7,
-           LobattoIIIa2, LobattoIIIa4, LobattoIIIa5,
-           LobattoIIIb2, LobattoIIIb3, LobattoIIIb4, LobattoIIIb5,
-           LobattoIIIc2, LobattoIIIc3, LobattoIIIc4, LobattoIIIc5]
-weird = [MIRK2, MIRK5, RadauIIa2]
-daesolvers = []
 # Simple ODESystem with BVP constraints.
 let
     @parameters α=1.5 β=1.0 γ=3.0 δ=1.0
@@ -222,24 +216,14 @@ let
 
     # Testing that more complicated constr give correct solutions.
     constr = [y(.2) + x(.8) ~ 3., y(.3) ~ 2.]
+    @mtkbuild lksys = ODESystem(eqs, t; constraints = constr)
     bvp = SciMLBase.BVProblem{false, SciMLBase.FullSpecialize}(lksys, u0map, tspan; guesses)
     test_solvers(solvers, bvp, u0map, constr; dt = 0.05)
 
     constr = [α * β - x(.6) ~ 0.0, y(.2) ~ 3.]
+    @mtkbuild lksys = ODESystem(eqs, t; constraints = constr)
     bvp = SciMLBase.BVProblem{true, SciMLBase.AutoSpecialize}(lksys, u0map, tspan; guesses)
     test_solvers(solvers, bvp, u0map, constr)
-
-    # Testing that errors are properly thrown when malformed constr are given.
-    @variables bad(..)
-    constr = [x(1.) + bad(3.) ~ 10]
-    @test_throws ErrorException lksys = ODESystem(eqs, t; constraints = constr)
-
-    constr = [x(t) + y(t) ~ 3]
-    @test_throws ErrorException lksys = ODESystem(eqs, t; constraints = constr)
-
-    @parameters bad2
-    constr = [bad2 + x(0.) ~ 3]
-    @test_throws ErrorException lksys = ODESystem(eqs, t; constraints = constr)
 end
 
 # Cartesian pendulum from the docs.