Merge pull request #400 from SciML/io

[WIP] Causal Input/Output variables, explicit algebraic variables, aliases, and open connections

Author: ChrisRackauckas

src/ModelingToolkit.jl

@@ -132,7 +132,7 @@ export Differential, expand_derivatives, @derivatives
 export IntervalDomain, ProductDomain, ⊗, CircleDomain
 export Equation, ConstrainedEquation
 export Operation, Expression, Variable
-export independent_variable, states, parameters, equations
+export independent_variable, states, parameters, equations, pins, observed
 
 export calculate_jacobian, generate_jacobian, generate_function
 export calculate_tgrad, generate_tgrad

src/symutils.jl

@@ -118,6 +118,7 @@ Generate a function to evaluate the system's equations.
 function generate_function end
 
 function Base.getproperty(sys::AbstractSystem, name::Symbol)
+
     if name ∈ fieldnames(typeof(sys))
         return getfield(sys,name)
     elseif !isempty(sys.systems)
@@ -126,6 +127,7 @@ function Base.getproperty(sys::AbstractSystem, name::Symbol)
             return rename(sys.systems[i],renamespace(sys.name,name))
         end
     end
+
     i = findfirst(x->x.name==name,sys.states)
     if i !== nothing
         x = rename(sys.states[i],renamespace(sys.name,name))
@@ -135,12 +137,21 @@ function Base.getproperty(sys::AbstractSystem, name::Symbol)
             return x()
         end
     end
+
     if :ps ∈ fieldnames(typeof(sys))
         i = findfirst(x->x.name==name,sys.ps)
         if i !== nothing
             return rename(sys.ps[i],renamespace(sys.name,name))()
         end
     end
+
+    if :observed ∈ fieldnames(typeof(sys))
+        i = findfirst(x->convert(Variable,x.lhs).name==name,sys.observed)
+        if i !== nothing
+            return rename(convert(Variable,sys.observed[i].lhs),renamespace(sys.name,name))(getfield(sys,:iv)())
+        end
+    end
+
     throw(error("Variable $name does not exist"))
 end
 
@@ -154,8 +165,13 @@ function namespace_parameters(sys::AbstractSystem)
     [rename(x,renamespace(sys.name,x.name)) for x in parameters(sys)]
 end
 
+function namespace_pins(sys::AbstractSystem)
+    [rename(x,renamespace(sys.name,x.name)) for x in pins(sys)]
+end
+
 namespace_equations(sys::AbstractSystem) = namespace_equation.(equations(sys),sys.name,sys.iv.name)
 
+
 function namespace_equation(eq::Equation,name,ivname)
     _lhs = namespace_operation(eq.lhs,name,ivname)
     _rhs = namespace_operation(eq.rhs,name,ivname)
@@ -172,11 +188,14 @@ end
 namespace_operation(O::Constant,name,ivname) = O
 
 independent_variable(sys::AbstractSystem) = sys.iv
-states(sys::AbstractSystem) = isempty(sys.systems) ? sys.states : [sys.states;reduce(vcat,namespace_variables.(sys.systems))]
+states(sys::AbstractSystem) = isempty(sys.systems) ? setdiff(sys.states, convert.(Variable,sys.pins)) : [sys.states;reduce(vcat,namespace_variables.(sys.systems))]
 parameters(sys::AbstractSystem) = isempty(sys.systems) ? sys.ps : [sys.ps;reduce(vcat,namespace_parameters.(sys.systems))]
-
-function equations(sys::AbstractSystem)
-    isempty(sys.systems) ? sys.eqs : [sys.eqs;reduce(vcat,namespace_equations.(sys.systems))]
+pins(sys::AbstractSystem) = isempty(sys.systems) ? sys.pins : [sys.pins;reduce(vcat,namespace_pins.(sys.systems))]
+function observed(sys::AbstractSystem)
+    [sys.observed;
+     reduce(vcat,
+            (namespace_equation.(s.observed, s.name, s.iv.name) for s in sys.systems),
+            init=Equation[])]
 end
 
 function states(sys::AbstractSystem,name::Symbol)
@@ -189,6 +208,34 @@ function parameters(sys::AbstractSystem,name::Symbol)
     rename(x,renamespace(sys.name,x.name))()
 end
 
+function pins(sys::AbstractSystem,name::Symbol)
+    x = sys.pins[findfirst(x->x.name==name,sys.ps)]
+    rename(x,renamespace(sys.name,x.name))(sys.iv())
+end
+
+lhss(xs) = map(x->x.lhs, xs)
+rhss(xs) = map(x->x.rhs, xs)
+
+function equations(sys::ModelingToolkit.AbstractSystem; remove_aliases = true)
+    if isempty(sys.systems)
+        return sys.eqs
+    else
+        eqs = [sys.eqs;
+               reduce(vcat,
+                      namespace_equations.(sys.systems);
+                      init=Equation[])]
+
+        if !remove_aliases
+            return eqs
+        end
+        aliases = observed(sys)
+        dict = Dict(lhss(aliases) .=> rhss(aliases))
+
+        # Substitute aliases
+        return Equation.(lhss(eqs), Rewriters.Fixpoint(x->substitute(x, dict)).(rhss(eqs)))
+    end
+end
+
 function states(sys::AbstractSystem,args...)
     name = last(args)
     extra_names = reduce(Symbol,[Symbol(:₊,x.name) for x in args[1:end-1]])
@@ -212,6 +259,13 @@ function islinear(sys::AbstractSystem)
     all(islinear(r, dvs) for r in rhs)
 end
 
+function pins(sys::AbstractSystem,args...)
+    name = last(args)
+    extra_names = reduce(Symbol,[Symbol(:₊,x.name) for x in args[1:end-1]])
+    newname = renamespace(extra_names,name)
+    rename(x,renamespace(sys.name,newname))(sys.iv())
+end
+
 struct AbstractSysToExpr
     sys::AbstractSystem
     states::Vector{Variable}

src/systems/abstractsystem.jl

@@ -31,6 +31,8 @@ struct ODESystem <: AbstractODESystem
     states::Vector{Variable}
     """Parameter variables."""
     ps::Vector{Variable}
+    pins::Vector{Variable}
+    observed::Vector{Equation}
     """
     Time-derivative matrix. Note: this field will not be defined until
     [`calculate_tgrad`](@ref) is called on the system.
@@ -62,6 +64,8 @@ struct ODESystem <: AbstractODESystem
 end
 
 function ODESystem(deqs::AbstractVector{<:Equation}, iv, dvs, ps;
+                   pins = Variable[],
+                   observed = Operation[],
                    systems = ODESystem[],
                    name=gensym(:ODESystem))
     iv′ = convert(Variable,iv)
@@ -71,7 +75,7 @@ function ODESystem(deqs::AbstractVector{<:Equation}, iv, dvs, ps;
     jac = RefValue{Any}(Matrix{Expression}(undef, 0, 0))
     Wfact   = RefValue(Matrix{Expression}(undef, 0, 0))
     Wfact_t = RefValue(Matrix{Expression}(undef, 0, 0))
-    ODESystem(deqs, iv′, dvs′, ps′, tgrad, jac, Wfact, Wfact_t, name, systems)
+    ODESystem(deqs, iv′, dvs′, ps′, pins, observed, tgrad, jac, Wfact, Wfact_t, name, systems)
 end
 
 var_from_nested_derivative(x::Constant) = (missing, missing)

src/systems/diffeqs/odesystem.jl

@@ -37,6 +37,8 @@ struct SDESystem <: AbstractODESystem
     states::Vector{Variable}
     """Parameter variables."""
     ps::Vector{Variable}
+    pins::Vector{Variable}
+    observed::Vector{Equation}
     """
     Time-derivative matrix. Note: this field will not be defined until
     [`calculate_tgrad`](@ref) is called on the system.
@@ -68,6 +70,8 @@ struct SDESystem <: AbstractODESystem
 end
 
 function SDESystem(deqs::AbstractVector{<:Equation}, neqs, iv, dvs, ps;
+                   pins = Variable[],
+                   observed = Operation[],
                    systems = SDESystem[],
                    name = gensym(:SDESystem))
     iv′ = convert(Variable,iv)
@@ -77,7 +81,7 @@ function SDESystem(deqs::AbstractVector{<:Equation}, neqs, iv, dvs, ps;
     jac = RefValue(Matrix{Expression}(undef, 0, 0))
     Wfact   = RefValue(Matrix{Expression}(undef, 0, 0))
     Wfact_t = RefValue(Matrix{Expression}(undef, 0, 0))
-    SDESystem(deqs, neqs, iv′, dvs′, ps′, tgrad, jac, Wfact, Wfact_t, name, systems)
+    SDESystem(deqs, neqs, iv′, dvs′, ps′, pins, observed, tgrad, jac, Wfact, Wfact_t, name, systems)
 end
 
 function generate_diffusion_function(sys::SDESystem, dvs = sys.states, ps = sys.ps; kwargs...)

src/systems/diffeqs/sdesystem.jl

@@ -37,14 +37,20 @@ struct JumpSystem{U <: ArrayPartition} <: AbstractSystem
     states::Vector{Variable}
     """The parameters of the system."""
     ps::Vector{Variable}
+    pins::Vector{Variable}
+    observed::Vector{Equation}
     """The name of the system."""
     name::Symbol
     """The internal systems."""
     systems::Vector{JumpSystem}
 end
 
-function JumpSystem(eqs, iv, states, ps; systems = JumpSystem[],
-                                          name = gensym(:JumpSystem))
+function JumpSystem(eqs, iv, states, ps;
+                    pins = Variable[],
+                    observed = Equation[],
+                    systems = JumpSystem[],
+                    name = gensym(:JumpSystem))
+
     ap = ArrayPartition(MassActionJump[], ConstantRateJump[], VariableRateJump[])
     for eq in eqs
         if eq isa MassActionJump
@@ -58,13 +64,9 @@ function JumpSystem(eqs, iv, states, ps; systems = JumpSystem[],
         end
     end
 
-    JumpSystem{typeof(ap)}(ap, convert(Variable,iv), convert.(Variable, states), convert.(Variable, ps), name, systems)
+    JumpSystem{typeof(ap)}(ap, convert(Variable,iv), convert.(Variable, states), convert.(Variable, ps), pins, observed, name, systems)
 end
 
-JumpSystem(eqs::ArrayPartition, iv, states, ps; systems = JumpSystem[], name = gensym(:JumpSystem)) =
-    JumpSystem{typeof(eqs)}(eqs, convert(Variable,iv), convert.(Variable, states), convert.(Variable, ps), name, systems)
-
-
 generate_rate_function(js, rate) = build_function(rate, states(js), parameters(js),
                                         independent_variable(js),
                                         expression=Val{true})

src/systems/jumps/jumpsystem.jl

@@ -25,6 +25,8 @@ struct NonlinearSystem <: AbstractSystem
     states::Vector{Variable}
     """Parameters."""
     ps::Vector{Variable}
+    pins::Vector{Variable}
+    observed::Vector{Equation}
     """
     Name: the name of the system
     """
@@ -36,9 +38,11 @@ struct NonlinearSystem <: AbstractSystem
 end
 
 function NonlinearSystem(eqs, states, ps;
+                         pins = Variable[],
+                         observed = Operation[],
                          name = gensym(:NonlinearSystem),
                          systems = NonlinearSystem[])
-    NonlinearSystem(eqs, convert.(Variable,states), convert.(Variable,ps), name, systems)
+    NonlinearSystem(eqs, convert.(Variable,states), convert.(Variable,ps), pins, observed, name, systems)
 end
 
 function calculate_jacobian(sys::NonlinearSystem;sparse=false,simplify=true)

src/systems/nonlinear/nonlinearsystem.jl

@@ -23,6 +23,8 @@ struct OptimizationSystem <: AbstractSystem
     states::Vector{Variable}
     """Parameters."""
     ps::Vector{Variable}
+    pins::Vector{Variable}
+    observed::Vector{Equation}
     """
     Name: the name of the system
     """
@@ -34,9 +36,11 @@ struct OptimizationSystem <: AbstractSystem
 end
 
 function OptimizationSystem(op, states, ps;
+                            pins = Variable[],
+                            observed = Operation[],
                             name = gensym(:OptimizationSystem),
                             systems = OptimizationSystem[])
-    OptimizationSystem(op, convert.(Variable,states), convert.(Variable,ps), name, systems)
+    OptimizationSystem(op, convert.(Variable,states), convert.(Variable,ps), pins, observed, name, systems)
 end
 
 function calculate_gradient(sys::OptimizationSystem)

src/systems/optimization/optimizationsystem.jl

@@ -129,20 +129,25 @@ struct ReactionSystem <: AbstractSystem
     states::Vector{Variable}
     """Parameter variables."""
     ps::Vector{Variable}
+    pins::Vector{Variable}
+    observed::Vector{Equation}
     """The name of the system"""
     name::Symbol
     """systems: The internal systems"""
     systems::Vector{ReactionSystem}
 end
 
-function ReactionSystem(eqs, iv, species, params; systems = ReactionSystem[],
-                                                  name = gensym(:ReactionSystem))
-
+function ReactionSystem(eqs, iv, species, params;
+                        pins = Variable[],
+                        observed = Operation[],
+                        systems = ReactionSystem[],
+                        name = gensym(:ReactionSystem))
 
     isempty(species) && error("ReactionSystems require at least one species.")
     paramvars = map(v -> convert(Variable,v), params)
     specvars  = map(s -> convert(Variable,s), species)
-    ReactionSystem(eqs, convert(Variable,iv), specvars, paramvars, name, systems)
+    ReactionSystem(eqs, convert(Variable,iv), specvars, paramvars,
+                   pins, observed, name, systems)
 end
 
 """
@@ -294,6 +299,7 @@ explicitly on the independent variable (usually time).
 function ismassaction(rx, rs; rxvars = get_variables(rx.rate),
                               haveivdep = any(var -> isequal(rs.iv,convert(Variable,var)), rxvars),
                               stateset = Set(states(rs)))
+    return !(haveivdep || rx.only_use_rate || any(convert(Variable,rxv) in states(rs) for rxv in rxvars))
     # if no dependencies must be zero order
     (length(rxvars)==0) && return true
     (haveivdep || rx.only_use_rate) && return false

src/systems/reaction/reactionsystem.jl

@@ -0,0 +1,39 @@
+using ModelingToolkit, OrdinaryDiffEq, Test
+
+@parameters t σ ρ β
+@variables x(t) y(t) z(t) F(t) u(t)
+@derivatives D'~t
+
+eqs = [D(x) ~ σ*(y-x) + F,
+       D(y) ~ x*(ρ-z)-y,
+       D(z) ~ x*y - β*z]
+
+aliases = [u ~ x + y - z]
+lorenz1 = ODESystem(eqs,pins=[F],observed=aliases,name=:lorenz1)
+lorenz2 = ODESystem(eqs,pins=[F],observed=aliases,name=:lorenz2)
+
+connections = [lorenz1.F ~ lorenz2.u,
+               lorenz2.F ~ lorenz1.u]
+connected = ODESystem(Equation[],t,[],[],observed=connections,systems=[lorenz1,lorenz2])
+
+sys = connected
+
+@variables lorenz1₊F lorenz2₊F
+@test pins(connected) == Variable[lorenz1₊F, lorenz2₊F]
+@test isequal(observed(connected),
+              [connections...,
+               lorenz1.u ~ lorenz1.x + lorenz1.y - lorenz1.z,
+               lorenz2.u ~ lorenz2.x + lorenz2.y - lorenz2.z])
+
+collapsed_eqs = [D(lorenz1.x) ~ (lorenz1.σ * (lorenz1.y - lorenz1.x) +
+                                 (lorenz2.x + lorenz2.y - lorenz2.z)),
+                 D(lorenz1.y) ~ lorenz1.x * (lorenz1.ρ - lorenz1.z) - lorenz1.y,
+                 D(lorenz1.z) ~ lorenz1.x * lorenz1.y - (lorenz1.β * lorenz1.z),
+                 D(lorenz2.x) ~ (lorenz2.σ * (lorenz2.y - lorenz2.x) +
+                                 (lorenz1.x + lorenz1.y - lorenz1.z)),
+                 D(lorenz2.y) ~ lorenz2.x * (lorenz2.ρ - lorenz2.z) - lorenz2.y,
+                 D(lorenz2.z) ~ lorenz2.x * lorenz2.y - (lorenz2.β * lorenz2.z)]
+
+simplifyeqs(eqs) = Equation.((x->x.lhs).(eqs), simplify.((x->x.rhs).(eqs)))
+
+@test isequal(simplifyeqs(equations(connected)), simplifyeqs(collapsed_eqs))