New abstract time dependent/independent types & independent_variables method (#1091)

Author: lamorton

docs/src/basics/AbstractSystem.md

@@ -7,6 +7,13 @@ It establishes a common set of functionality that is used between systems
 from ODEs and chemical reactions, allowing users to have a common framework for
 model manipulation and compilation.
 
+### Subtypes
+
+There are three immediate subtypes of `AbstractSystem`, classified by how many independent variables each type has:
+* `AbstractTimeIndependentSystem`: has no independent variable (eg: `NonlinearSystem`)
+* `AbstractTimeDependentSystem`: has a single independent variable (eg: `ODESystem`)
+* `AbstractMultivariateSystem`: may have multiple independent variables (eg: `PDESystem`)
+
 ## Constructors and Naming
 
 The `AbstractSystem` interface has a consistent method for constructing systems.
@@ -50,12 +57,11 @@ Optionally, a system could have:
 - `observed(sys)`: All observed equations of the system and its subsystems.
 - `get_observed(sys)`: Observed equations of the current-level system.
 - `get_defaults(sys)`: A `Dict` that maps variables into their default values.
-- `independent_variable(sys)`: The independent variable of a system.
+- `independent_variables(sys)`: The independent variables of a system.
 - `get_noiseeqs(sys)`: Noise equations of the current-level system.
 
-Note that there's `get_iv(sys)`, but it is not advised to use, since it errors
-when the system has no field `iv`. `independent_variable(sys)` returns `nothing`
-for `NonlinearSystem`s.
+Note that if you know a system is an `AbstractTimeDependentSystem` you could use `get_iv` to get the 
+unique independent variable directly, rather than using `independenent_variables(sys)[1]`, which is clunky and may cause problems if `sys` is an `AbstractMultivariateSystem` because there may be more than one independent variable. `AbstractTimeIndependentSystem`s do not have a method `get_iv`, and `independent_variables(sys)` will return a size-zero result for such. For an `AbstractMultivariateSystem`, `get_ivs` is equivalent.
 
 A system could also have caches:
 

docs/src/systems/JumpSystem.md

@@ -11,7 +11,7 @@ JumpSystem
 - `get_eqs(sys)` or `equations(sys)`: The equations that define the jump system.
 - `get_states(sys)` or `states(sys)`: The set of states in the jump system.
 - `get_ps(sys)` or `parameters(sys)`: The parameters of the jump system.
-- `independent_variable(sys)`: The independent variable of the jump system.
+- `get_iv(sys)`: The independent variable of the jump system.
 
 ## Transformations
 

docs/src/systems/ODESystem.md

@@ -11,7 +11,7 @@ ODESystem
 - `get_eqs(sys)` or `equations(sys)`: The equations that define the ODE.
 - `get_states(sys)` or `states(sys)`: The set of states in the ODE.
 - `get_ps(sys)` or `parameters(sys)`: The parameters of the ODE.
-- `independent_variable(sys)`: The independent variable of the ODE.
+- `get_iv(sys)`: The independent variable of the ODE.
 
 ## Transformations
 

docs/src/systems/ReactionSystem.md

@@ -49,7 +49,7 @@ ReactionSystem
 - `get_eqs(sys)` or `equations(sys)`: The reactions that define the system.
 - `get_states(sys)` or `states(sys)`: The set of chemical species in the system.
 - `get_ps(sys)` or `parameters(sys)`: The parameters of the system.
-- `independent_variable(sys)`: The independent variable of the
+- `get_iv(sys)`: The independent variable of the
   reaction system, usually time.
 
 ## Query Functions

docs/src/systems/SDESystem.md

@@ -17,7 +17,7 @@ sde = SDESystem(ode, noiseeqs)
 - `get_eqs(sys)` or `equations(sys)`: The equations that define the SDE.
 - `get_states(sys)` or `states(sys)`: The set of states in the SDE.
 - `get_ps(sys)` or `parameters(sys)`: The parameters of the SDE.
-- `independent_variable(sys)`: The independent variable of the SDE.
+- `get_iv(sys)`: The independent variable of the SDE.
 
 ## Transformations
 

src/ModelingToolkit.jl

@@ -77,7 +77,10 @@ $(TYPEDEF)
 TODO
 """
 abstract type AbstractSystem end
-abstract type AbstractODESystem <: AbstractSystem end
+abstract type AbstractTimeDependentSystem <: AbstractSystem end
+abstract type AbstractTimeIndependentSystem <: AbstractSystem end
+abstract type AbstractODESystem <: AbstractTimeDependentSystem end
+abstract type AbstractMultivariateSystem <: AbstractSystem end
 
 """
 $(TYPEDSIGNATURES)
@@ -86,6 +89,8 @@ Get the set of independent variables for the given system.
 """
 function independent_variables end
 
+function independent_variable end
+
 """
 $(TYPEDSIGNATURES)
 
@@ -148,6 +153,7 @@ end
 
 struct Flow end
 
+export AbstractTimeDependentSystem, AbstractTimeIndependentSystem, AbstractMultivariateSystem
 export ODESystem, ODEFunction, ODEFunctionExpr, ODEProblemExpr, convert_system
 export DAEFunctionExpr, DAEProblemExpr
 export SDESystem, SDEFunction, SDEFunctionExpr, SDESystemExpr
@@ -170,7 +176,7 @@ export Differential, expand_derivatives, @derivatives
 export Equation, ConstrainedEquation
 export Term, Sym
 export SymScope, LocalScope, ParentScope, GlobalScope
-export independent_variable, states, parameters, equations, controls, observed, structure
+export independent_variables, independent_variable, states, parameters, equations, controls, observed, structure
 export structural_simplify
 export DiscreteSystem, DiscreteProblem
 

src/structural_transformation/StructuralTransformations.jl

@@ -13,10 +13,10 @@ using SymbolicUtils.Rewriters
 using SymbolicUtils: similarterm, istree
 
 using ModelingToolkit
-using ModelingToolkit: ODESystem, AbstractSystem, var_from_nested_derivative,
-                       Differential, states, equations, vars, Symbolic,
-                       diff2term, value, operation, arguments, Sym, Term,
-                       simplify, solve_for, isdiffeq, isdifferential,
+using ModelingToolkit: ODESystem, AbstractSystem,var_from_nested_derivative, Differential,
+                       states, equations, vars, Symbolic, diff2term, value,
+                       operation, arguments, Sym, Term, simplify, solve_for,
+                       isdiffeq, isdifferential, get_structure, get_iv, independent_variables,
                        get_structure, defaults, InvalidSystemException,
                        ExtraEquationsSystemException,
                        ExtraVariablesSystemException

src/structural_transformation/codegen.jl

@@ -208,7 +208,7 @@ function build_torn_function(
               out
               DestructuredArgs(states, inbounds=!checkbounds)
               DestructuredArgs(parameters(sys), inbounds=!checkbounds)
-              independent_variable(sys)
+              independent_variables(sys)
              ],
              [],
              Let(
@@ -312,7 +312,7 @@ function build_observed_function(
         [
          DestructuredArgs(diffvars, inbounds=!checkbounds)
          DestructuredArgs(parameters(sys), inbounds=!checkbounds)
-         independent_variable(sys)
+         independent_variables(sys)
         ],
         [],
         Let(

src/structural_transformation/pantelides.jl

@@ -2,7 +2,7 @@
 ### Reassemble: structural information -> system
 ###
 
-function pantelides_reassemble(sys, eqassoc, assign)
+function pantelides_reassemble(sys::ODESystem, eqassoc, assign)
     s = structure(sys)
     @unpack fullvars, varassoc = s
     # Step 1: write derivative equations
@@ -15,7 +15,7 @@ function pantelides_reassemble(sys, eqassoc, assign)
     fill!(out_vars, nothing)
     out_vars[1:length(fullvars)] .= fullvars
 
-    D = Differential(independent_variable(sys))
+    D = Differential(get_iv(sys))
 
     for (i, v) in enumerate(varassoc)
         # fullvars[v] = D(fullvars[i])
@@ -75,11 +75,11 @@ end
 
 Perform Pantelides algorithm.
 """
-function pantelides!(sys; maxiters = 8000)
+function pantelides!(sys::ODESystem; maxiters = 8000)
     s = structure(sys)
     # D(j) = assoc[j]
     @unpack graph, fullvars, varassoc = s
-    iv = independent_variable(sys)
+    iv = get_iv(sys)
     neqs = nsrcs(graph)
     nvars = length(varassoc)
     vcolor = falses(nvars)

src/systems/abstractsystem.jl

@@ -1,6 +1,6 @@
 """
 ```julia
-calculate_tgrad(sys::AbstractSystem)
+calculate_tgrad(sys::AbstractTimeDependentSystem)
 ```
 
 Calculate the time gradient of a system.
@@ -72,7 +72,7 @@ function calculate_hessian end
 
 """
 ```julia
-generate_tgrad(sys::AbstractSystem, dvs = states(sys), ps = parameters(sys), expression = Val{true}; kwargs...)
+generate_tgrad(sys::AbstractTimeDependentSystem, dvs = states(sys), ps = parameters(sys), expression = Val{true}; kwargs...)
 ```
 
 Generates a function for the time gradient of a system. Extra arguments control
@@ -138,8 +138,28 @@ function getname(t)
         nameof(t)
     end
 end
+#Deprecated
+function independent_variable(sys::AbstractSystem)
+    Base.depwarn("`independent_variable` is deprecated. Use `get_iv` or `independent_variables` instead.",:independent_variable)
+    isdefined(sys, :iv) ? getfield(sys, :iv) : nothing
+end
 
-independent_variable(sys::AbstractSystem) = isdefined(sys, :iv) ? getfield(sys, :iv) : nothing
+#Treat the result as a vector of symbols always
+function independent_variables(sys::AbstractSystem)
+    systype = typeof(sys)
+    @warn "Please declare ($systype) as a subtype of `AbstractTimeDependentSystem`, `AbstractTimeIndependentSystem` or `AbstractMultivariateSystem`."
+    if isdefined(sys, :iv)
+        return [getfield(sys, :iv)]
+    elseif isdefined(sys, :ivs)
+        return getfield(sys,:ivs)
+    else
+        return []
+    end
+end
+
+independent_variables(sys::AbstractTimeDependentSystem) = [getfield(sys, :iv)]
+independent_variables(sys::AbstractTimeIndependentSystem) = []
+independent_variables(sys::AbstractMultivariateSystem) = getfield(sys, :ivs)
 
 function structure(sys::AbstractSystem)
     s = get_structure(sys)
@@ -333,35 +353,35 @@ namespace_controls(sys::AbstractSystem) = controls(sys, controls(sys))
 
 function namespace_defaults(sys)
     defs = defaults(sys)
-    Dict((isparameter(k) ? parameters(sys, k) : states(sys, k)) => namespace_expr(defs[k], nameof(sys), independent_variable(sys)) for k in keys(defs))
+    Dict((isparameter(k) ? parameters(sys, k) : states(sys, k)) => namespace_expr(defs[k], nameof(sys), independent_variables(sys)) for k in keys(defs))
 end
 
 function namespace_equations(sys::AbstractSystem)
     eqs = equations(sys)
     isempty(eqs) && return Equation[]
-    iv = independent_variable(sys)
-    map(eq->namespace_equation(eq,nameof(sys),iv), eqs)
+    ivs = independent_variables(sys)
+    map(eq -> namespace_equation(eq, nameof(sys), ivs), eqs)
 end
 
-function namespace_equation(eq::Equation,name,iv)
-    _lhs = namespace_expr(eq.lhs,name,iv)
-    _rhs = namespace_expr(eq.rhs,name,iv)
+function namespace_equation(eq::Equation, name, ivs)
+    _lhs = namespace_expr(eq.lhs, name, ivs)
+    _rhs = namespace_expr(eq.rhs, name, ivs)
     _lhs ~ _rhs
 end
 
-function namespace_expr(O::Sym,name,iv)
-    isequal(O, iv) ? O : renamespace(name,O)
+function namespace_expr(O::Sym, name, ivs)
+    any(isequal(O), ivs) ? O : renamespace(name, O)
 end
 
 _symparam(s::Symbolic{T}) where {T} = T
-function namespace_expr(O,name,iv) where {T}
+function namespace_expr(O, name, ivs) where {T}
     O = value(O)
     if istree(O)
-        renamed = map(a->namespace_expr(a,name,iv), arguments(O))
+        renamed = map(a -> namespace_expr(a, name, ivs), arguments(O))
         if operation(O) isa Sym
             renamespace(name, O)
         else
-            similarterm(O,operation(O),renamed)
+            similarterm(O, operation(O), renamed)
         end
     else
         O
@@ -389,13 +409,13 @@ function controls(sys::AbstractSystem)
 end
 
 function observed(sys::AbstractSystem)
-    iv = independent_variable(sys)
+    ivs = independent_variables(sys)
     obs = get_observed(sys)
     systems = get_systems(sys)
     [obs;
      reduce(vcat,
-            (map(o->namespace_equation(o, nameof(s), iv), observed(s)) for s in systems),
-            init=Equation[])]
+            (map(o -> namespace_equation.(o, nameof(s), ivs), observed(s)) for s in systems),
+            init = Equation[])]
 end
 
 Base.@deprecate default_u0(x) defaults(x) false
@@ -514,9 +534,10 @@ function toexpr(sys::AbstractSystem)
     stmt = expr.args
 
     name = Meta.quot(nameof(sys))
-    iv = independent_variable(sys)
+    ivs = independent_variables(sys)
     ivname = gensym(:iv)
-    if iv !== nothing
+    for iv in ivs
+        ivname = gensym(:iv)
         push!(stmt, :($ivname = (@variables $(getname(iv)))[1]))
     end
 
@@ -536,9 +557,12 @@ function toexpr(sys::AbstractSystem)
     defs_name = push_defaults!(stmt, defaults(sys), var2name)
 
     if sys isa ODESystem
-        push!(stmt, :($ODESystem($eqs_name, $ivname, $stsname, $psname; defaults=$defs_name, name=$name)))
+        iv = get_iv(sys)
+        ivname = gensym(:iv)
+        push!(stmt, :($ivname = (@variables $(getname(iv)))[1]))
+        push!(stmt, :($ODESystem($eqs_name, $ivname, $stsname, $psname; defaults = $defs_name, name=$name)))
     elseif sys isa NonlinearSystem
-        push!(stmt, :($NonlinearSystem($eqs_name, $stsname, $psname; defaults=$defs_name, name=$name)))
+        push!(stmt, :($NonlinearSystem($eqs_name, $stsname, $psname; defaults = $defs_name, name=$name)))
     end
 
     striplines(expr) # keeping the line numbers is never helpful
@@ -891,7 +915,7 @@ function Base.hash(sys::AbstractSystem, s::UInt)
     s = foldr(hash, get_ps(sys), init=s)
     s = foldr(hash, get_eqs(sys), init=s)
     s = foldr(hash, get_observed(sys), init=s)
-    s = hash(independent_variable(sys), s)
+    s = hash(independent_variables(sys), s)
     return s
 end
 
@@ -903,11 +927,13 @@ by default.
 """
 function extend(sys::AbstractSystem, basesys::AbstractSystem; name::Symbol=nameof(sys))
     T = SciMLBase.parameterless_type(basesys)
-    iv = independent_variable(basesys)
-    if iv === nothing
+    ivs = independent_variables(basesys)
+    if length(ivs) == 0
         sys = convert_system(T, sys)
+    elseif length(ivs) == 1
+        sys = convert_system(T, sys, ivs[1])
     else
-        sys = convert_system(T, sys, iv)
+        throw("Extending multivariate systems is not supported")
     end
 
     eqs = union(equations(basesys), equations(sys))
@@ -917,10 +943,10 @@ function extend(sys::AbstractSystem, basesys::AbstractSystem; name::Symbol=nameo
     defs = merge(defaults(basesys), defaults(sys)) # prefer `sys`
     syss = union(get_systems(basesys), get_systems(sys))
 
-    if iv === nothing
-        T(eqs, sts, ps, observed=obs, defaults=defs, name=name, systems=syss)
-    else
-        T(eqs, iv, sts, ps, observed=obs, defaults=defs, name=name, systems=syss)
+    if length(ivs) == 0
+        T(eqs, sts, ps, observed = obs, defaults = defs, name=name, systems = syss)
+    elseif length(ivs) == 1
+        T(eqs, ivs[1], sts, ps, observed = obs, defaults = defs, name = name, systems = syss)
     end
 end
 

src/systems/control/controlsystem.jl

@@ -1,4 +1,4 @@
-abstract type AbstractControlSystem <: AbstractSystem end
+abstract type AbstractControlSystem <: AbstractTimeDependentSystem end
 
 function namespace_controls(sys::AbstractControlSystem)
     [rename(x,renamespace(nameof(sys),nameof(x))) for x in controls(sys)]