Merge pull request #85 from JuliaDiffEq/hg/refactor/variable

Refactor :DependentVariable into :Unknown

Author: ChrisRackauckas

README.md

@@ -21,15 +21,15 @@ to manipulate.
 ### Example: ODE
 
 Let's build an ODE. First we define some variables. In a differential equation
-system, we need to differentiate between our dependent variables, independent
+system, we need to differentiate between our unknown (dependent) variables, independent
 variables, and parameters. Therefore we label them as follows:
 
 ```julia
 using ModelingToolkit
 
 # Define some variables
 @IVar t
-@DVar x(t) y(t) z(t)
+@Unknown x(t) y(t) z(t)
 @Deriv D'~t
 @Param σ ρ β
 ```
@@ -47,7 +47,7 @@ Each operation builds an `Operation` type, and thus `eqs` is an array of
 analyzed by other programs. We can turn this into a `DiffEqSystem` via:
 
 ```julia
-de = DiffEqSystem(eqs,[t],[x,y,z],Variable[],[σ,ρ,β])
+de = DiffEqSystem(eqs,[t],[x,y,z],[σ,ρ,β])
 de = DiffEqSystem(eqs)
 ```
 
@@ -87,12 +87,10 @@ f = ODEFunction(de)
 
 We can also build nonlinear systems. Let's say we wanted to solve for the steady
 state of the previous ODE. This is the nonlinear system defined by where the
-derivatives are zero. We could use dependent variables for our nonlinear system
-(for direct compatibility with the above ODE code), or we can use non-tagged
-variables. Here we will show the latter. We write:
+derivatives are zero. We use unknown variables for our nonlinear system.
 
 ```julia
-@Var x y z
+@Unknown x y z
 @Param σ ρ β
 
 # Define a nonlinear system
@@ -191,7 +189,7 @@ structure is as follows:
   the system of equations.
 - Name to subtype mappings: these describe how variable `subtype`s are mapped
   to the contexts of the system. For example, for a differential equation,
-  the dependent variable corresponds to given subtypes and then the `eqs` can
+  the unknown variable corresponds to given subtypes and then the `eqs` can
   be analyzed knowing what the state variables are.
 - Variable names which do not fall into one of the system's core subtypes are
   treated as intermediates which can be used for holding subcalculations and
@@ -262,7 +260,7 @@ syntactic sugar in some form. For example, the variable construction:
 
 ```julia
 @IVar t
-@DVar x(t) y(t) z(t)
+@Unknown x(t) y(t) z(t)
 @Deriv D'~t
 @Param σ ρ β
 ```
@@ -271,9 +269,9 @@ is syntactic sugar for:
 
 ```julia
 t = IndependentVariable(:t)
-x = DependentVariable(:x,dependents = [t])
-y = DependentVariable(:y,dependents = [t])
-z = DependentVariable(:z,dependents = [t])
+x = Unknown(:x, dependents = [t])
+y = Unknown(:y, dependents = [t])
+z = Unknown(:z, dependents = [t])
 D = Differential(t) # Default of first derivative, Derivative(t,1)
 σ = Parameter(:σ)
 ρ = Parameter(:ρ)
@@ -285,10 +283,10 @@ D = Differential(t) # Default of first derivative, Derivative(t,1)
 The system building functions can handle intermediate calculations. For example,
 
 ```julia
-@Var a x y z
+@Unknown x y z
 @Param σ ρ β
-eqs = [a ~ y-x,
-       0 ~ σ*a,
+a = y - x
+eqs = [0 ~ σ*a,
        0 ~ x*(ρ-z)-y,
        0 ~ x*y - β*z]
 ns = NonlinearSystem(eqs,[x,y,z],[σ,ρ,β])

src/equations.jl

@@ -12,40 +12,36 @@ Base.:~(lhs::Expression, rhs::Number    ) = Equation(lhs, rhs)
 Base.:~(lhs::Number    , rhs::Expression) = Equation(lhs, rhs)
 
 
-function extract_elements(eqs, targetmap, default = nothing)
-    elems = Dict{Symbol, Vector{Variable}}()
-    names = Dict{Symbol, Set{Symbol}}()
-    if default == nothing
-        targets =  unique(collect(values(targetmap)))
-    else
-        targets = [unique(collect(values(targetmap))), default]
-    end
-    for target in targets
-        elems[target] = Vector{Variable}()
-        names[target] = Set{Symbol}()
-    end
-    for eq in eqs
-        extract_elements!(eq, elems, names, targetmap, default)
+_is_derivative(x::Variable) = x.diff !== nothing
+_is_dependent(x::Variable) = x.subtype === :Unknown && !isempty(x.dependents)
+_subtype(subtype::Symbol) = x -> x.subtype === subtype
+
+function extract_elements(eqs, predicates)
+    result = [Variable[] for p ∈ predicates]
+    vars = foldl(vars!, eqs; init=Set{Variable}())
+
+    for var ∈ vars
+        for (i, p) ∈ enumerate(predicates)
+            p(var) && (push!(result[i], var); break)
+        end
     end
-    Tuple(elems[target] for target in targets)
+
+    return result
 end
-# Walk the tree recursively and push variables into the right set
-function extract_elements!(op, elems, names, targetmap, default)
+
+function vars!(vars, op)
     args = isa(op, Equation) ? Expression[op.lhs, op.rhs] : op.args
 
-    for arg in args
-        if arg isa Operation
-            extract_elements!(arg, elems, names, targetmap, default)
-        elseif arg isa Variable
-            if default == nothing
-                target = haskey(targetmap, arg.subtype) ? targetmap[arg.subtype] : continue
-            else
-                target = haskey(targetmap, arg.subtype) ? targetmap[arg.subtype] : default
-            end
-            if !in(arg.name, names[target])
-                push!(names[target], arg.name)
-                push!(elems[target], arg)
+    for arg ∈ args
+        if isa(arg, Operation)
+            vars!(vars, arg)
+        elseif isa(arg, Variable)
+            push!(vars, arg)
+            for dep ∈ arg.dependents
+                push!(vars, dep)
             end
         end
     end
+
+    return vars
 end

src/systems/diffeqs/diffeqsystem.jl

@@ -2,51 +2,35 @@ mutable struct DiffEqSystem <: AbstractSystem
     eqs::Vector{Equation}
     ivs::Vector{Variable}
     dvs::Vector{Variable}
-    vs::Vector{Variable}
     ps::Vector{Variable}
-    iv_name::Symbol
-    dv_name::Symbol
-    p_name::Symbol
     jac::Matrix{Expression}
 end
 
-function DiffEqSystem(eqs, ivs, dvs, vs, ps)
-    iv_name = ivs[1].subtype
-    dv_name = dvs[1].subtype
-    p_name = isempty(ps) ? :Parameter : ps[1].subtype
-    DiffEqSystem(eqs, ivs, dvs, vs, ps, iv_name, dv_name, p_name, Matrix{Expression}(undef,0,0))
-end
+DiffEqSystem(eqs, ivs, dvs, ps) = DiffEqSystem(eqs, ivs, dvs, ps, Matrix{Expression}(undef,0,0))
 
-function DiffEqSystem(eqs; iv_name = :IndependentVariable,
-                           dv_name = :DependentVariable,
-                           v_name = :Variable,
-                           p_name = :Parameter)
-    targetmap =  Dict(iv_name => iv_name, dv_name => dv_name, v_name => v_name,
-                       p_name => p_name)
-    ivs, dvs, vs, ps = extract_elements(eqs, targetmap)
-    DiffEqSystem(eqs, ivs, dvs, vs, ps, iv_name, dv_name, p_name, Matrix{Expression}(0,0))
+function DiffEqSystem(eqs)
+    predicates = [_is_derivative, _subtype(:IndependentVariable), _is_dependent, _subtype(:Parameter)]
+    _, ivs, dvs, ps = extract_elements(eqs, predicates)
+    DiffEqSystem(eqs, ivs, dvs, ps, Matrix{Expression}(undef,0,0))
 end
 
-function DiffEqSystem(eqs, ivs;
-                      dv_name = :DependentVariable,
-                      v_name = :Variable,
-                      p_name = :Parameter)
-    targetmap =  Dict(dv_name => dv_name, v_name => v_name, p_name => p_name)
-    dvs, vs, ps = extract_elements(eqs, targetmap)
-    DiffEqSystem(eqs, ivs, dvs, vs, ps, ivs[1].subtype, dv_name, p_name, Matrix{Expression}(undef,0,0))
+function DiffEqSystem(eqs, ivs)
+    predicates = [_is_derivative, _is_dependent, _subtype(:Parameter)]
+    _, dvs, ps = extract_elements(eqs, predicates)
+    DiffEqSystem(eqs, ivs, dvs, ps, Matrix{Expression}(undef,0,0))
 end
 
 function generate_ode_function(sys::DiffEqSystem;version = ArrayFunction)
-    var_exprs = [:($(sys.dvs[i].name) = u[$i]) for i in 1:length(sys.dvs)]
-    param_exprs = [:($(sys.ps[i].name) = p[$i]) for i in 1:length(sys.ps)]
+    var_exprs = [:($(sys.dvs[i].name) = u[$i]) for i in eachindex(sys.dvs)]
+    param_exprs = [:($(sys.ps[i].name) = p[$i]) for i in eachindex(sys.ps)]
     sys_exprs = build_equals_expr.(sys.eqs)
     if version == ArrayFunction
-        dvar_exprs = [:(du[$i] = $(Symbol("$(sys.dvs[i].name)_$(sys.ivs[1].name)"))) for i in 1:length(sys.dvs)]
+        dvar_exprs = [:(du[$i] = $(Symbol("$(sys.dvs[i].name)_$(sys.ivs[1].name)"))) for i in eachindex(sys.dvs)]
         exprs = vcat(var_exprs,param_exprs,sys_exprs,dvar_exprs)
         block = expr_arr_to_block(exprs)
         :((du,u,p,t)->$(toexpr(block)))
     elseif version == SArrayFunction
-        dvar_exprs = [:($(Symbol("$(sys.dvs[i].name)_$(sys.ivs[1].name)"))) for i in 1:length(sys.dvs)]
+        dvar_exprs = [:($(Symbol("$(sys.dvs[i].name)_$(sys.ivs[1].name)"))) for i in eachindex(sys.dvs)]
         svector_expr = quote
             E = eltype(tuple($(dvar_exprs...)))
             T = StaticArrays.similar_type(typeof(u), E)
@@ -84,8 +68,8 @@ function calculate_jacobian(sys::DiffEqSystem, simplify=true)
 end
 
 function generate_ode_jacobian(sys::DiffEqSystem, simplify=true)
-    var_exprs = [:($(sys.dvs[i].name) = u[$i]) for i in 1:length(sys.dvs)]
-    param_exprs = [:($(sys.ps[i].name) = p[$i]) for i in 1:length(sys.ps)]
+    var_exprs = [:($(sys.dvs[i].name) = u[$i]) for i in eachindex(sys.dvs)]
+    param_exprs = [:($(sys.ps[i].name) = p[$i]) for i in eachindex(sys.ps)]
     diff_exprs = filter(!isintermediate, sys.eqs)
     jac = calculate_jacobian(sys, simplify)
     sys.jac = jac
@@ -96,12 +80,12 @@ function generate_ode_jacobian(sys::DiffEqSystem, simplify=true)
 end
 
 function generate_ode_iW(sys::DiffEqSystem, simplify=true)
-    var_exprs = [:($(sys.dvs[i].name) = u[$i]) for i in 1:length(sys.dvs)]
-    param_exprs = [:($(sys.ps[i].name) = p[$i]) for i in 1:length(sys.ps)]
+    var_exprs = [:($(sys.dvs[i].name) = u[$i]) for i in eachindex(sys.dvs)]
+    param_exprs = [:($(sys.ps[i].name) = p[$i]) for i in eachindex(sys.ps)]
     diff_exprs = filter(!isintermediate, sys.eqs)
     jac = sys.jac
 
-    gam = Variable(:gam)
+    gam = Unknown(:gam)
 
     W = LinearAlgebra.I - gam*jac
     W = SMatrix{size(W,1),size(W,2)}(W)

src/systems/nonlinear/nonlinear_system.jl

@@ -2,25 +2,11 @@ struct NonlinearSystem <: AbstractSystem
     eqs::Vector{Equation}
     vs::Vector{Variable}
     ps::Vector{Variable}
-    v_name::Vector{Symbol}
-    p_name::Symbol
 end
 
-function NonlinearSystem(eqs, vs, ps;
-                         v_name = :Variable,
-                         dv_name = :DependentVariable,
-                         p_name = :Parameter)
-    NonlinearSystem(eqs, vs, ps, [v_name,dv_name], p_name)
-end
-
-function NonlinearSystem(eqs;
-                         v_name = :Variable,
-                         dv_name = :DependentVariable,
-                         p_name = :Parameter)
-    # Allow the use of :DependentVariable to make it seamless with DE use
-    targetmap = Dict(v_name => v_name, dv_name => v_name, p_name => p_name)
-    vs, ps = extract_elements(eqs, targetmap)
-    NonlinearSystem(eqs, vs, ps, [v_name,dv_name], p_name)
+function NonlinearSystem(eqs)
+    vs, ps = extract_elements(eqs, [_subtype(:Unknown), _subtype(:Parameter)])
+    NonlinearSystem(eqs, vs, ps)
 end
 
 function generate_nlsys_function(sys::NonlinearSystem)

src/variables.jl

@@ -5,20 +5,15 @@ mutable struct Variable <: Expression
     dependents::Vector{Variable}
 end
 
-Variable(name; subtype::Symbol=:Variable, dependents::Vector{Variable} = Variable[]) =
+Variable(name; subtype::Symbol, dependents::Vector{Variable} = Variable[]) =
     Variable(name, subtype, nothing, dependents)
-Variable(name, args...; kwargs...) = Variable(name, args...; subtype=:Variable, kwargs...)
 
 Parameter(name,args...;kwargs...) = Variable(name,args...;subtype=:Parameter,kwargs...)
 IndependentVariable(name,args...;kwargs...) = Variable(name,args...;subtype=:IndependentVariable,kwargs...)
+Unknown(name,args...;kwargs...) = Variable(name,args...;subtype=:Unknown,kwargs...)
 
-function DependentVariable(name,args...;dependents = [],kwargs...)
-    @assert !isempty(dependents)
-    Variable(name,args...;subtype=:DependentVariable,dependents=dependents,kwargs...)
-end
-
-export Variable,Parameter,Constant,DependentVariable,IndependentVariable,
-       @Var, @Param, @Const, @DVar, @IVar
+export Variable,Parameter,Constant,Unknown,IndependentVariable,
+       @Param, @Const, @Unknown, @IVar
 
 
 Base.copy(x::Variable) = Variable(x.name, x.subtype, x.diff, x.dependents)
@@ -72,7 +67,7 @@ function _parse_vars(macroname, fun, x)
         @assert iscall || issym "@$macroname expects a tuple of expressions!\nE.g. `@$macroname x y z`"
 
         if iscall
-            dependents = :([$(_var.args[2:end]...)])
+            dependents = :(Variable[$(_var.args[2:end]...)])
             lhs = _var.args[1]
         else
             dependents = Variable[]
@@ -88,9 +83,7 @@ function _parse_vars(macroname, fun, x)
     return ex
 end
 
-for funs in ((:DVar, :DependentVariable), (:IVar, :IndependentVariable),
-             (:Var, :Variable),
-             (:Param, :Parameter))
+for funs in [(:Unknown, :Unknown), (:IVar, :IndependentVariable), (:Param, :Parameter)]
     @eval begin
         macro ($(funs[1]))(x...)
             esc(_parse_vars(String($funs[1]), $funs[2], x))

test/ambiguity.jl

@@ -2,7 +2,7 @@ using ModelingToolkit
 using Test
 
 @IVar t
-@DVar x(t) y(t) z(t)
+@Unknown x(t) y(t) z(t)
 
 struct __MyType__ end
 Base.:~(::__MyType__,::Number) = 2

test/basic_variables_and_operations.jl

@@ -2,16 +2,15 @@ using ModelingToolkit
 using Test
 
 @IVar t
-@DVar x(t) y(t) z(t)
-@test_throws AssertionError @DVar w
+@Unknown x(t) y(t) z(t)
 
 @Deriv D'~t # Default of first derivative, Derivative(t,1)
 @Param σ ρ β
 @Const c=0
 
 # Default values
 p = Parameter(:p)
-u = DependentVariable(:u, dependents = [t])
+u = Unknown(:u, dependents = [t])
 
 σ*(y-x)
 D(x)

test/components.jl

@@ -17,7 +17,7 @@ function generate_lorenz_eqs(t,x,y,z,σ,ρ,β)
      D(z) ~ x*y - β*z]
 end
 function Lorenz(t)
-    @DVar x(t) y(t) z(t)
+    @Unknown x(t) y(t) z(t)
     @Param σ ρ β
     Lorenz(x, y, z, σ, ρ, β, generate_lorenz_eqs(t, x, y, z, σ, ρ, β))
 end

test/derivatives.jl

@@ -3,7 +3,7 @@ using Test
 
 # Derivatives
 @IVar t
-@Var x(t) y(t) z(t)
+@Unknown x(t) y(t) z(t)
 @Param σ ρ β
 @Deriv D'~t
 dsin = D(sin(t))
@@ -43,5 +43,5 @@ jac = ModelingToolkit.calculate_jacobian(sys)
 @test jac[3,3] == -1*β
 
 # Variable dependence checking in differentiation
-@Var a(t) b(a)
+@Unknown a(t) b(a)
 @test D(b) ≠ Constant(0)

test/internal.jl

@@ -4,8 +4,7 @@ using Test
 # `Expr`, `Number` -> `Operation`
 @IVar a
 @Param b
-@DVar x(t)
-@Var y
+@Unknown x(t) y()
 @test convert(Expression, 2) == 2
 expr = :(-inv(2sqrt(+($a, $b))))
 op   = Operation(-, [Operation(inv,