vector valued inputs

Author: TorkelE

test/sciml_problem_inputs.jl

@@ -252,4 +252,185 @@ let
             end
         end
     end
+end
+
+### Vector Parameter/Variable Inputs ###
+
+begin 
+    # Declares the model (with vector species/parameters, with/without default values, and observables).
+    @variables X(t)[1:2] Y(t)[1:2] = [10.0, 20.0] XY(t)[1:2]
+    @parameters p[1:2] d[1:2] = [0.2, 0.5]
+    alg_eqs = [
+        0 ~ p[1] - d[1]*X[1],
+        0 ~ p[2] - d[2]*X[2],
+        0 ~ p[1] - d[1]*Y[1],
+        0 ~ p[2] - d[2]*Y[2],
+    ]
+    diff_eqs = [
+        D(X[1]) ~ p[1] - d[1]*X[1],
+        D(X[2]) ~ p[2] - d[2]*X[2],
+        D(Y[1]) ~ p[1] - d[1]*Y[1],
+        D(Y[2]) ~ p[2] - d[2]*Y[2],
+    ]
+    noise_eqs = fill(0.01, 4, 8)
+    jumps = [
+        MassActionJump(p[1], Pair{Symbolics.BasicSymbolic{Real}, Int64}[], [X[1] => 1]),
+        MassActionJump(p[2], Pair{Symbolics.BasicSymbolic{Real}, Int64}[], [X[2] => 1]),
+        MassActionJump(d[1], [X[1] => 1], [X[1] => -1]),
+        MassActionJump(d[2], [X[2] => 1], [X[2] => -1]),
+        MassActionJump(p[1], Pair{Symbolics.BasicSymbolic{Real}, Int64}[], [Y[1] => 1]),
+        MassActionJump(p[2], Pair{Symbolics.BasicSymbolic{Real}, Int64}[], [Y[2] => 1]),
+        MassActionJump(d[1], [Y[1] => 1], [Y[1] => -1]),
+        MassActionJump(d[2], [Y[2] => 1], [Y[2] => -1])
+    ]
+    observed = [XY[1] ~ X[1] + Y[1], XY[2] ~ X[2] + Y[2]]
+
+    # Create systems (without structural_simplify, since that might modify systems to affect intended tests).
+    osys = complete(ODESystem(diff_eqs, t; observed, name = :osys))
+    ssys = complete(SDESystem(diff_eqs, noise_eqs, t, [X[1], X[2], Y[1], Y[2]], [p, d]; observed, name = :ssys))
+    jsys = complete(JumpSystem(jumps, t, [X[1], X[2], Y[1], Y[2]], [p, d]; observed, name = :jsys))
+    nsys = complete(NonlinearSystem(alg_eqs; observed, name = :nsys))
+
+    # Declares various u0 versions (scalarised and vector forms).
+    u0_alts_vec = [
+        # Vectors not providing default values.
+        [X => [1.0, 2.0]],
+        [X[1] => 1.0, X[2] => 2.0],
+        [osys.X => [1.0, 2.0]],
+        [osys.X[1] => 1.0, osys.X[2] => 2.0],
+        # Vectors providing default values.
+        [X => [1.0, 2.0], Y => [10.0, 20.0]],
+        [X[1] => 1.0, X[2] => 2.0, Y[1] => 10.0, Y[2] => 20.0],
+        [osys.X => [1.0, 2.0], osys.Y => [10.0, 20.0]],
+        [osys.X[1] => 1.0, osys.X[2] => 2.0, osys.Y[1] => 10.0, osys.Y[2] => 20.0],
+        # Dicts not providing default values.
+        Dict([X => [1.0, 2.0]]),
+        Dict([X[1] => 1.0, X[2] => 2.0]),
+        Dict([osys.X => [1.0, 2.0]]),
+        Dict([osys.X[1] => 1.0, osys.X[2] => 2.0]),
+        # Dicts providing default values.
+        Dict([X => [1.0, 2.0], Y => [10.0, 20.0]]),
+        Dict([X[1] => 1.0, X[2] => 2.0, Y[1] => 10.0, Y[2] => 20.0]),
+        Dict([osys.X => [1.0, 2.0], osys.Y => [10.0, 20.0]]),
+        Dict([osys.X[1] => 1.0, osys.X[2] => 2.0, osys.Y[1] => 10.0, osys.Y[2] => 20.0]),
+        # Tuples not providing default values.
+        (X => [1.0, 2.0]),
+        (X[1] => 1.0, X[2] => 2.0),
+        (osys.X => [1.0, 2.0]),
+        (osys.X[1] => 1.0, osys.X[2] => 2.0),
+        # Tuples providing default values.
+        (X => [1.0, 2.0], Y => [10.0, 20.0]),
+        (X[1] => 1.0, X[2] => 2.0, Y[1] => 10.0, Y[2] => 20.0),
+        (osys.X => [1.0, 2.0], osys.Y => [10.0, 20.0]),
+        (osys.X[1] => 1.0, osys.X[2] => 2.0, osys.Y[1] => 10.0, osys.Y[2] => 20.0)
+    ]
+
+    # Declares various ps versions (vector forms only).
+    p_alts_vec = [
+        # Vectors not providing default values.
+        [p => [1.0, 2.0]],
+        [osys.p => [1.0, 2.0]],
+        # Vectors providing default values.
+        [p => [4.0, 5.0], d => [0.2, 0.5]],
+        [osys.p => [4.0, 5.0], osys.d => [0.2, 0.5]],
+        # Dicts not providing default values.
+        Dict([p => [1.0, 2.0]]),
+        Dict([osys.p => [1.0, 2.0]]),
+        # Dicts providing default values.
+        Dict([p => [4.0, 5.0], d => [0.2, 0.5]]),
+        Dict([osys.p => [4.0, 5.0], osys.d => [0.2, 0.5]]),
+        # Tuples not providing default values.
+        (p => [1.0, 2.0]),
+        (osys.p => [1.0, 2.0]),
+        # Tuples providing default values.
+        (p => [4.0, 5.0], d => [0.2, 0.5]),
+        (osys.p => [4.0, 5.0], osys.d => [0.2, 0.5])
+    ]
+
+    # Declares a timespan.
+    tspan = (0.0, 10.0)
+end
+
+# Perform ODE simulations (singular and ensemble).
+# Fails. At least partially related to https://github.com/SciML/ModelingToolkit.jl/issues/2804.
+@test_broken let 
+    # Creates normal and ensemble problems.
+    base_oprob = ODEProblem(osys, u0_alts_vec[1], tspan, p_alts_vec[1])
+    base_sol = solve(base_oprob, Tsit5(); saveat = 1.0)
+    base_eprob = EnsembleProblem(base_oprob)
+    base_esol = solve(base_eprob, Tsit5(); trajectories = 2, saveat = 1.0)
+
+    # Simulates problems for all input types, checking that identical solutions are found.
+    for u0 in u0_alts_vec, p in p_alts_vec
+        oprob = remake(base_oprob; u0, p)
+        @test base_sol == solve(oprob, Tsit5(); saveat = 1.0)
+        eprob = remake(base_eprob; u0, p)
+        @test base_esol == solve(eprob, Tsit5(); trajectories = 2, saveat = 1.0)
+    end
+end
+
+# Perform SDE simulations (singular and ensemble).
+# Fails. At least partially related to https://github.com/SciML/ModelingToolkit.jl/issues/2804.
+@test_broken let 
+    # Creates normal and ensemble problems.
+    base_sprob = SDEProblem(ssys, u0_alts_vec[1], tspan, p_alts_vec[1])
+    base_sol = solve(base_sprob, ImplicitEM(); seed, saveat = 1.0)
+    base_eprob = EnsembleProblem(base_sprob)
+    base_esol = solve(base_eprob, ImplicitEM(); seed, trajectories = 2, saveat = 1.0)
+
+    # Simulates problems for all input types, checking that identical solutions are found.
+    for u0 in u0_alts_vec, p in p_alts_vec
+        sprob = remake(base_sprob; u0, p)
+        @test base_sol == solve(sprob, ImplicitEM(); seed, saveat = 1.0)
+        eprob = remake(base_eprob; u0, p)
+        @test base_esol == solve(eprob, ImplicitEM(); seed, trajectories = 2, saveat = 1.0)
+    end
+end
+
+# Perform jump simulations (singular and ensemble).
+# Fails. At least partially related to https://github.com/SciML/ModelingToolkit.jl/issues/2804.
+@test_broken let 
+    # Creates normal and ensemble problems.
+    base_dprob = DiscreteProblem(jsys, u0_alts_vec[1], tspan, p_alts_vec[1])
+    base_jprob = JumpProblem(jsys, base_dprob, Direct(); rng)
+    base_sol = solve(base_jprob, SSAStepper(); seed, saveat = 1.0)
+    base_eprob = EnsembleProblem(base_jprob)
+    base_esol = solve(base_eprob, SSAStepper(); seed, trajectories = 2, saveat = 1.0)
+
+    # Simulates problems for all input types, checking that identical solutions are found.
+    for u0 in u0_alts_vec, p in p_alts_vec
+        jprob = remake(base_jprob; u0, p)
+        @test base_sol == solve(base_jprob, SSAStepper(); seed, saveat = 1.0)
+        eprob = remake(base_eprob; u0, p)
+        @test base_esol == solve(eprob, SSAStepper(); seed, trajectories = 2, saveat = 1.0)
+    end
+end
+
+# Solves a nonlinear problem (EnsembleProblems are not possible for these).
+# Fails. At least partially related to https://github.com/SciML/ModelingToolkit.jl/issues/2804.
+@test_broken let
+    base_nlprob = NonlinearProblem(nsys, u0_alts_vec[1], p_alts_vec[1])
+    base_sol = solve(base_nlprob, NewtonRaphson())
+    for u0 in u0_alts_vec, p in p_alts_vec
+        nlprob = remake(base_nlprob; u0, p)
+        @test base_sol == solve(nlprob, NewtonRaphson())
+    end
+end
+
+# Perform steady state simulations (singular and ensemble).
+# Fails. At least partially related to https://github.com/SciML/ModelingToolkit.jl/issues/2804.
+@test_broken let 
+    # Creates normal and ensemble problems.
+    base_ssprob = SteadyStateProblem(osys, u0_alts_vec[1], p_alts_vec[1])
+    base_sol = solve(base_ssprob, DynamicSS(Tsit5()))
+    base_eprob = EnsembleProblem(base_ssprob)
+    base_esol = solve(base_eprob, DynamicSS(Tsit5()); trajectories = 2)
+
+    # Simulates problems for all input types, checking that identical solutions are found.
+    for u0 in u0_alts_vec, p in p_alts_vec
+        ssprob = remake(base_ssprob; u0, p)
+        @test base_sol == solve(ssprob, DynamicSS(Tsit5()))
+        eprob = remake(base_eprob; u0, p)
+        @test base_esol == solve(eprob, DynamicSS(Tsit5()); trajectories = 2)
+    end
 end