Merge branch 'master' into sm/refactor

Author: YingboMa

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

@@ -78,7 +78,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)
@@ -87,6 +90,8 @@ Get the set of independent variables for the given system.
 """
 function independent_variables end
 
+function independent_variable end
+
 """
 $(TYPEDSIGNATURES)
 
@@ -149,6 +154,7 @@ end
 
 struct Flow end
 
+export AbstractTimeDependentSystem, AbstractTimeIndependentSystem, AbstractMultivariateSystem
 export ODESystem, ODEFunction, ODEFunctionExpr, ODEProblemExpr, convert_system
 export DAEFunctionExpr, DAEProblemExpr
 export SDESystem, SDEFunction, SDEFunctionExpr, SDESystemExpr
@@ -171,7 +177,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
@@ -131,7 +131,28 @@ function generate_function end
 
 Base.nameof(sys::AbstractSystem) = getfield(sys, :name)
 
-independent_variable(sys::AbstractSystem) = isdefined(sys, :iv) ? getfield(sys, :iv) : nothing
+#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
+
+#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)
@@ -365,9 +386,9 @@ function namespace_assignment(eq::Assignment, sys)
 end
 
 function namespace_expr(O, sys) where {T}
-    iv = independent_variable(sys)
+    ivs = independent_variables(sys)
     O = unwrap(O)
-    if isequal(O, iv)
+    if any(isequal(O), ivs)
         return O
     elseif isvariable(O)
         renamespace(sys, O)
@@ -547,9 +568,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
 
@@ -569,9 +591,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
@@ -928,7 +953,7 @@ function Base.hash(sys::AbstractSystem, s::UInt)
         s = foldr(hash, get_eqs(sys), init=s)
     end
     s = foldr(hash, get_observed(sys), init=s)
-    s = hash(independent_variable(sys), s)
+    s = hash(independent_variables(sys), s)
     return s
 end
 
@@ -940,11 +965,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))
@@ -954,10 +981,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)]

src/systems/diffeqs/abstractodesystem.jl

@@ -83,7 +83,7 @@ check_derivative_variables(eq) = check_operator_variables(eq, Differential)
 function generate_function(
         sys::AbstractODESystem, dvs = states(sys), ps = parameters(sys);
         implicit_dae=false,
-        ddvs=implicit_dae ? map(Differential(independent_variable(sys)), dvs) : nothing,
+        ddvs=implicit_dae ? map(Differential(get_iv(sys)), dvs) : nothing,
         kwargs...
     )
     # optimization
@@ -155,15 +155,15 @@ function generate_difference_cb(sys::ODESystem, dvs = states(sys), ps = paramete
     PeriodicCallback(cb_affect!, first(dts))
 end
 
-function time_varying_as_func(x, sys)
+function time_varying_as_func(x, sys::AbstractTimeDependentSystem)
     # if something is not x(t) (the current state)
     # but is `x(t-1)` or something like that, pass in `x` as a callable function rather
     # than pass in a value in place of x(t).
     #
     # This is done by just making `x` the argument of the function.
     if istree(x) &&
         operation(x) isa Sym &&
-        !(length(arguments(x)) == 1 && isequal(arguments(x)[1], independent_variable(sys)))
+        !(length(arguments(x)) == 1 && isequal(arguments(x)[1], get_iv(sys)))
         return operation(x)
     end
     return x
@@ -299,7 +299,7 @@ function DiffEqBase.ODEFunction{iip}(sys::AbstractODESystem, dvs = states(sys),
                      mass_matrix = _M,
                      jac_prototype = jac_prototype,
                      syms = Symbol.(states(sys)),
-                     indepsym = Symbol(independent_variable(sys)),
+                     indepsym = Symbol(get_iv(sys)),
                      observed = observedfun,
                     )
 end
@@ -320,7 +320,7 @@ respectively.
 """
 function DiffEqBase.DAEFunction{iip}(sys::AbstractODESystem, dvs = states(sys),
                                      ps = parameters(sys), u0 = nothing;
-                                     ddvs=map(diff2term ∘ Differential(independent_variable(sys)), dvs),
+                                     ddvs=map(diff2term ∘ Differential(get_iv(sys)), dvs),
                                      version = nothing,
                                      #=
                                      tgrad=false,
@@ -356,7 +356,7 @@ function DiffEqBase.DAEFunction{iip}(sys::AbstractODESystem, dvs = states(sys),
                      f,
                      syms = Symbol.(dvs),
                      # missing fields in `DAEFunction`
-                     #indepsym = Symbol(independent_variable(sys)),
+                     #indepsym = Symbol(get_iv(sys)),
                      #observed = observedfun,
                     )
 end
@@ -453,7 +453,7 @@ function ODEFunctionExpr{iip}(sys::AbstractODESystem, dvs = states(sys),
                           mass_matrix = M,
                           jac_prototype = $jp_expr,
                           syms = $(Symbol.(states(sys))),
-                          indepsym = $(QuoteNode(Symbol(independent_variable(sys)))),
+                          indepsym = $(QuoteNode(Symbol(get_iv(sys)))),
                          )
     end
     !linenumbers ? striplines(ex) : ex
@@ -472,7 +472,7 @@ function process_DEProblem(constructor, sys::AbstractODESystem,u0map,parammap;
     dvs = states(sys)
     ps = parameters(sys)
     defs = defaults(sys)
-    iv = independent_variable(sys)
+    iv = get_iv(sys)
     if parammap isa Dict
         u0defs = merge(parammap, defs)
     elseif eltype(parammap) <: Pair

src/systems/diffeqs/basic_transformations.jl

@@ -45,8 +45,8 @@ Optimal Transport Approach
 Abhishek Halder, Kooktae Lee, and Raktim Bhattacharya
 https://abhishekhalder.bitbucket.io/F16ACC2013Final.pdf
 """
-function liouville_transform(sys)
-      t = independent_variable(sys)
+function liouville_transform(sys::AbstractODESystem)
+      t = get_iv(sys)
       @variables trJ
       D = ModelingToolkit.Differential(t)
       neweq = D(trJ) ~ trJ*-tr(calculate_jacobian(sys))

src/systems/diffeqs/first_order_transform.jl

@@ -5,7 +5,7 @@ Takes a Nth order ODESystem and returns a new ODESystem written in first order
 form by defining new variables which represent the N-1 derivatives.
 """
 function ode_order_lowering(sys::ODESystem)
-    iv = independent_variable(sys)
+    iv = get_iv(sys)
     eqs_lowered, new_vars = ode_order_lowering(equations(sys), iv, states(sys))
     @set! sys.eqs = eqs_lowered
     @set! sys.states = new_vars