Merge branch 'master' of github.com:SciML/ModelingToolkit.jl into domain-sets

Author: valentinsulzer

Project.toml

@@ -1,7 +1,7 @@
 name = "ModelingToolkit"
 uuid = "961ee093-0014-501f-94e3-6117800e7a78"
 authors = ["Chris Rackauckas <[email protected]>"]
-version = "5.17.1"
+version = "5.17.4"
 
 [deps]
 AbstractTrees = "1520ce14-60c1-5f80-bbc7-55ef81b5835c"

src/systems/abstractsystem.jl

@@ -665,15 +665,17 @@ AbstractTrees.children(sys::ModelingToolkit.AbstractSystem) = ModelingToolkit.ge
 AbstractTrees.printnode(io::IO, sys::ModelingToolkit.AbstractSystem) = print(io, nameof(sys))
 AbstractTrees.nodetype(::ModelingToolkit.AbstractSystem) = ModelingToolkit.AbstractSystem
 
-function check_eqs_u0(eqs, dvs, u0)
+function check_eqs_u0(eqs, dvs, u0; check_length=true, kwargs...)
     if u0 !== nothing
-        if !(length(eqs) == length(dvs) == length(u0))
-            throw(ArgumentError("Equations ($(length(eqs))), states ($(length(dvs))), and initial conditions ($(length(u0))) are of different lengths."))
-        end
-    else
-        if !(length(eqs) == length(dvs))
-            throw(ArgumentError("Equations ($(length(eqs))), states ($(length(dvs))) are of different lengths."))
+        if check_length
+            if !(length(eqs) == length(dvs) == length(u0))
+                throw(ArgumentError("Equations ($(length(eqs))), states ($(length(dvs))), and initial conditions ($(length(u0))) are of different lengths. To allow a different number of equations than states use kwarg check_length=false."))
+            end
+        elseif length(dvs) != length(u0)
+            throw(ArgumentError("States ($(length(dvs))) and initial conditions ($(length(u0))) are of different lengths."))
         end
+    elseif check_length && (length(eqs) != length(dvs))
+        throw(ArgumentError("Equations ($(length(eqs))) and states ($(length(dvs))) are of different lengths. To allow these to differ use kwarg check_length=false."))
     end
     return nothing
 end
@@ -713,7 +715,7 @@ end
 
 promote_connect_rule(::Type{T}, ::Type{S}) where {T, S} = Union{}
 promote_connect_rule(::Type{T}, ::Type{T}) where {T} = T
-promote_connect_type(t1::Type, t2::Type, ts::Type...) = promote_connect_rule(promote_connect_rule(t1, t2), ts...)
+promote_connect_type(t1::Type, t2::Type, ts::Type...) = promote_connect_type(promote_connect_rule(t1, t2), ts...)
 @inline function promote_connect_type(::Type{T}, ::Type{S}) where {T,S}
     promote_connect_result(
         T,

src/systems/diffeqs/abstractodesystem.jl

@@ -426,7 +426,7 @@ function process_DEProblem(constructor, sys::AbstractODESystem,u0map,parammap;
     end
     p = varmap_to_vars(parammap,ps; defaults=pdefs)
 
-    check_eqs_u0(eqs, dvs, u0)
+    check_eqs_u0(eqs, dvs, u0; kwargs...)
 
     f = constructor(sys,dvs,ps,u0;ddvs=ddvs,tgrad=tgrad,jac=jac,checkbounds=checkbounds,
                     linenumbers=linenumbers,parallel=parallel,simplify=simplify,

src/systems/nonlinear/nonlinearsystem.jl

@@ -226,7 +226,7 @@ function process_NonlinearProblem(constructor, sys::NonlinearSystem,u0map,paramm
     u0 = varmap_to_vars(u0map,dvs; defaults=defs)
     p = varmap_to_vars(parammap,ps; defaults=defs)
 
-    check_eqs_u0(eqs, dvs, u0)
+    check_eqs_u0(eqs, dvs, u0; kwargs...)
 
     f = constructor(sys,dvs,ps,u0;jac=jac,checkbounds=checkbounds,
                     linenumbers=linenumbers,parallel=parallel,simplify=simplify,

src/systems/reaction/reactionsystem.jl

@@ -185,7 +185,7 @@ end
 """
     oderatelaw(rx; combinatoric_ratelaw=true)
 
-Given a [`Reaction`](@ref), return the reaction rate law [`Operation`](@ref) used in
+Given a [`Reaction`](@ref), return the symbolic reaction rate law used in
 generated ODEs for the reaction. Note, for a reaction defined by
 
 `k*X*Y, X+Z --> 2X + Y`
@@ -195,7 +195,7 @@ of the form
 
 `k, 2X+3Y --> Z`
 
-the `Operation` that is returned will be `k * (X(t)^2/2) * (Y(t)^3/6)`.
+the expression that is returned will be `k * (X(t)^2/2) * (Y(t)^3/6)`.
 
 Notes:
 - Allocates
@@ -240,13 +240,13 @@ function assemble_oderhs(rs; combinatoric_ratelaws=true)
     rhsvec
 end
 
-function assemble_drift(rs; combinatoric_ratelaws=true, as_odes=true)
+function assemble_drift(rs; combinatoric_ratelaws=true, as_odes=true, include_zero_odes=true)
     rhsvec = assemble_oderhs(rs; combinatoric_ratelaws=combinatoric_ratelaws)
     if as_odes
         D   = Differential(get_iv(rs))
-        eqs = [Equation(D(x),rhs) for (x,rhs) in zip(get_states(rs),rhsvec) if (!_iszero(rhs))]
+        eqs = [Equation(D(x),rhs) for (x,rhs) in zip(get_states(rs),rhsvec) if (include_zero_odes || (!_iszero(rhs)))]
     else
-        eqs = [Equation(0,rhs) for rhs in rhsvec if (!_iszero(rhs))]
+        eqs = [Equation(0,rhs) for rhs in rhsvec if (include_zero_odes || (!_iszero(rhs)))]
     end
     eqs
 end
@@ -272,7 +272,7 @@ end
 """
     jumpratelaw(rx; rxvars=get_variables(rx.rate), combinatoric_ratelaw=true)
 
-Given a [`Reaction`](@ref), return the reaction rate law [`Operation`](@ref) used in
+Given a [`Reaction`](@ref), return the symbolic reaction rate law used in
 generated stochastic chemical kinetics model SSAs for the reaction. Note,
 for a reaction defined by
 
@@ -283,7 +283,7 @@ the form
 
 `k, 2X+3Y --> Z`
 
-the `Operation` that is returned will be `k * binomial(X,2) *
+the expression that is returned will be `k * binomial(X,2) *
 binomial(Y,3)`.
 
 Notes:
@@ -411,8 +411,8 @@ law, i.e. for `2S -> 0` at rate `k` the ratelaw would be `k*S^2/2!`. If
 ignored.
 """
 function Base.convert(::Type{<:ODESystem}, rs::ReactionSystem; 
-                      name=nameof(rs), combinatoric_ratelaws=true, kwargs...)
-    eqs     = assemble_drift(rs; combinatoric_ratelaws=combinatoric_ratelaws)
+                      name=nameof(rs), combinatoric_ratelaws=true, include_zero_odes=true, kwargs...)
+    eqs     = assemble_drift(rs; combinatoric_ratelaws=combinatoric_ratelaws, include_zero_odes=include_zero_odes)
     systems = map(sys -> (sys isa ODESystem) ? sys : convert(ODESystem, sys), get_systems(rs))
     ODESystem(eqs, get_iv(rs), get_states(rs), get_ps(rs); name=name, systems=systems, kwargs...)
 end
@@ -431,8 +431,8 @@ law, i.e. for `2S -> 0` at rate `k` the ratelaw would be `k*S^2/2!`. If
 ignored.
 """
 function Base.convert(::Type{<:NonlinearSystem},rs::ReactionSystem;
-                      name=nameof(rs), combinatoric_ratelaws=true, kwargs...)
-    eqs     = assemble_drift(rs; combinatoric_ratelaws=combinatoric_ratelaws, as_odes=false)
+                      name=nameof(rs), combinatoric_ratelaws=true, include_zero_odes=true, kwargs...)
+    eqs     = assemble_drift(rs; combinatoric_ratelaws=combinatoric_ratelaws, as_odes=false, include_zero_odes=include_zero_odes)
     systems = convert.(NonlinearSystem, get_systems(rs))
     NonlinearSystem(eqs, get_states(rs), get_ps(rs); name=name, systems=systems, kwargs...)
 end
@@ -449,17 +449,18 @@ Notes:
 law, i.e. for `2S -> 0` at rate `k` the ratelaw would be `k*S^2/2!`. If
 `combinatoric_ratelaws=false` then the ratelaw is `k*S^2`, i.e. the scaling factor is
 ignored.
-- `noise_scaling=nothing::Union{Vector{Operation},Operation,Nothing}` allows for linear
+- `noise_scaling=nothing::Union{Vector{Num},Num,Nothing}` allows for linear
 scaling of the noise in the chemical Langevin equations. If `nothing` is given, the default
-value as in Gillespie 2000 is used. Alternatively, an `Operation` can be given, this is
+value as in Gillespie 2000 is used. Alternatively, a `Num` can be given, this is
 added as a parameter to the system (at the end of the parameter array). All noise terms
 are linearly scaled with this value. The parameter may be one already declared in the `ReactionSystem`.
-Finally, a `Vector{Operation}` can be provided (the length must be equal to the number of reactions).
+Finally, a `Vector{Num}` can be provided (the length must be equal to the number of reactions).
 Here the noise for each reaction is scaled by the corresponding parameter in the input vector.
 This input may contain repeat parameters.
 """
 function Base.convert(::Type{<:SDESystem}, rs::ReactionSystem; 
-                      noise_scaling=nothing, name=nameof(rs), combinatoric_ratelaws=true, kwargs...)
+                      noise_scaling=nothing, name=nameof(rs), combinatoric_ratelaws=true, 
+                      include_zero_odes=true, kwargs...)
 
     if noise_scaling isa Vector
         (length(noise_scaling)!=length(equations(rs))) &&
@@ -470,7 +471,8 @@ function Base.convert(::Type{<:SDESystem}, rs::ReactionSystem;
         noise_scaling = fill(value(noise_scaling),length(equations(rs)))
     end
 
-    eqs      = assemble_drift(rs; combinatoric_ratelaws=combinatoric_ratelaws)
+    eqs      = assemble_drift(rs; combinatoric_ratelaws=combinatoric_ratelaws, 
+                                  include_zero_odes=include_zero_odes)
     noiseeqs = assemble_diffusion(rs,noise_scaling;
                                   combinatoric_ratelaws=combinatoric_ratelaws)
     systems  = convert.(SDESystem, get_systems(rs))
@@ -506,13 +508,13 @@ end
 
 
 # ODEProblem from AbstractReactionNetwork
-function DiffEqBase.ODEProblem(rs::ReactionSystem, u0, tspan, p=DiffEqBase.NullParameters(), args...; kwargs...)
-    return ODEProblem(convert(ODESystem,rs; kwargs...),u0,tspan,p, args...; kwargs...)
+function DiffEqBase.ODEProblem(rs::ReactionSystem, u0, tspan, p=DiffEqBase.NullParameters(), args...; check_length=false, kwargs...)
+    return ODEProblem(convert(ODESystem,rs; kwargs...),u0,tspan,p, args...; check_length, kwargs...)
 end
 
 # NonlinearProblem from AbstractReactionNetwork
-function DiffEqBase.NonlinearProblem(rs::ReactionSystem, u0, p=DiffEqBase.NullParameters(), args...; kwargs...)
-    return NonlinearProblem(convert(NonlinearSystem,rs; kwargs...), u0, p, args...; kwargs...)
+function DiffEqBase.NonlinearProblem(rs::ReactionSystem, u0, p=DiffEqBase.NullParameters(), args...; check_length=false, kwargs...)
+    return NonlinearProblem(convert(NonlinearSystem,rs; kwargs...), u0, p, args...; check_length, kwargs...)
 end
 
 

test/connectors.jl

@@ -13,9 +13,19 @@ end
     ODESystem(Equation[], t, [x], [p], defaults=Dict(x=>1.0, p=>1.0))
 end
 
-ModelingToolkit.connect(::Type{<:Foo}, sys1, sys2) = [sys1.x ~ sys2.x]
+function ModelingToolkit.connect(::Type{<:Foo}, ss...)
+    n = length(ss)-1
+    eqs = Vector{Equation}(undef, n)
+    for i in 1:n
+        eqs[i] = ss[i].x ~ ss[i+1].x
+    end
+    eqs
+end
+
 @named f1 = Foo()
 @named f2 = Foo()
+@named f3 = Foo()
+@named f4 = Foo()
 @named g = Goo()
 
 @test isequal(connect(f1, f2), [f1.x ~ f2.x])
@@ -24,6 +34,9 @@ ModelingToolkit.connect(::Type{<:Foo}, sys1, sys2) = [sys1.x ~ sys2.x]
 # Note that since there're overloadings, these tests are not re-runable.
 ModelingToolkit.promote_connect_rule(::Type{<:Foo}, ::Type{<:Goo}) = Foo
 @test isequal(connect(f1, g), [f1.x ~ g.x])
+@test isequal(connect(f1, f2, g), [f1.x ~ f2.x; f2.x ~ g.x])
+@test isequal(connect(f1, f2, g, f3), [f1.x ~ f2.x; f2.x ~ g.x; g.x ~ f3.x])
+@test isequal(connect(f1, f2, g, f3, f4), [f1.x ~ f2.x; f2.x ~ g.x; g.x ~ f3.x; f3.x ~ f4.x])
 ModelingToolkit.promote_connect_rule(::Type{<:Goo}, ::Type{<:Foo}) = Foo
 @test isequal(connect(f1, g), [f1.x ~ g.x])
 # test conflict

test/odesystem.jl

@@ -330,3 +330,13 @@ sys = ode_order_lowering(sys)
 prob = ODEProblem(sys, [], tspan)
 sol = solve(prob, Tsit5())
 @test sum(abs, sol[end]) < 1
+
+
+# check_eqs_u0 kwarg test
+@parameters t
+@variables x1(t) x2(t)
+D =Differential(t)
+eqs = [D(x1) ~ -x1]
+sys = ODESystem(eqs,t,[x1,x2],[])
+@test_throws ArgumentError ODEProblem(sys, [1.0,1.0], (0.0,1.0))
+prob = ODEProblem(sys, [1.0,1.0], (0.0,1.0), check_length=false)

test/reactionsystem_components.jl

@@ -17,9 +17,9 @@ pars  = [α₀,α,K,n,δ,β,μ]
 @named rs = ReactionSystem(rxs, t, specs, pars)
 
 # using ODESystem components
-@named os₁ = convert(ODESystem, rs)
-@named os₂ = convert(ODESystem, rs)
-@named os₃ = convert(ODESystem, rs)
+@named os₁ = convert(ODESystem, rs; include_zero_odes=false)
+@named os₂ = convert(ODESystem, rs; include_zero_odes=false)
+@named os₃ = convert(ODESystem, rs; include_zero_odes=false)
 connections = [os₁.R ~ os₃.P,
                os₂.R ~ os₁.P,
                os₃.R ~ os₂.P]