Merge pull request #119 from JuliaDiffEq/hg/feature/varfuns

Refactor variables as functions

Author: ChrisRackauckas

README.md

@@ -43,11 +43,10 @@ eqs = [D(x) ~ σ*(y-x),
 
 Each operation builds an `Operation` type, and thus `eqs` is an array of
 `Operation` and `Variable`s. This holds a tree of the full system that can be
-analyzed by other programs. We can turn this into a `DiffEqSystem` via:
+analyzed by other programs. We can turn this into a `ODESystem` via:
 
 ```julia
-de = DiffEqSystem(eqs,t,[x,y,z],[σ,ρ,β])
-de = DiffEqSystem(eqs)
+de = ODESystem(eqs)
 ```
 
 where we tell it the variable types and ordering in the first version, or let it
@@ -56,7 +55,7 @@ This can then generate the function. For example, we can see the
 generated code via:
 
 ```julia
-generate_function(de)
+generate_function(de, [x,y,z], [σ,ρ,β])
 
 ## Which returns:
 :((##363, u, p, t)->begin
@@ -71,7 +70,7 @@ generate_function(de)
 and get the generated function via:
 
 ```julia
-f = ODEFunction(de)
+f = ODEFunction(de, [x,y,z], [σ,ρ,β])
 ```
 
 ### Example: Nonlinear System
@@ -88,8 +87,8 @@ derivatives are zero. We use (unknown) variables for our nonlinear system.
 eqs = [0 ~ σ*(y-x),
        0 ~ x*(ρ-z)-y,
        0 ~ x*y - β*z]
-ns = NonlinearSystem(eqs)
-nlsys_func = generate_function(ns)
+ns = NonlinearSystem(eqs, [x,y,z])
+nlsys_func = generate_function(ns, [x,y,z], [ρ,σ,β])
 ```
 
 which generates:
@@ -130,17 +129,49 @@ In this section we define the core pieces of the IR and what they mean.
 
 ### Variables
 
-The most fundamental part of the IR is the `Variable`. The `Variable` is the
+The most fundamental part of the IR is the `Variable`. In order to mirror the
+intention of solving for variables and representing function-like parameters,
+we treat each instance of `Variable` as a function which is called on its
+arguments using the natural syntax. Rather than having additional mechanisms
+for handling constant variables and parameters, we simply represent them as
+constant functions.
+
+The `Variable` is the
 context-aware single variable of the IR. Its fields are described as follows:
 
 - `name`: the name of the `Variable`. Note that this is not necessarily
   the same as the name of the Julia variable. But this symbol itself is considered
   the core identifier of the `Variable` in the sense of equality.
 - `known`: the main denotation of context, storing whether or not the value of
   the variable is known.
-- `dependents`: the vector of variables on which the current variable
-  is dependent. For example, `u(t,x)` has dependents `[t,x]`. Derivatives thus
-  require this information in order to simplify down.
+
+For example, the following code defines an independent variable `t`, a parameter
+`α`, a function parameter `σ`, a variable `x` which depends on `t`, a variable
+`y` with no dependents, and a variable `z` which depends on `t`, `α`, and `x(t)`.
+
+```julia
+t = Variable(:t; known = true)()  # independent variables are treated as known
+α = Variable(:α; known = true)()  # parameters are known
+σ = Variable(:σ; known = true)    # left uncalled, since it is used as a function
+w = Variable(:w; known = false)   # unknown, left uncalled
+x = Variable(:x; known = false)(t)  # unknown, depends on `t`
+y = Variable(:y; known = false)()   # unknown, no dependents
+z = Variable(:z; known = false)(t, α, x)  # unknown, multiple arguments
+
+expr = x + y^α + σ(3) * (z - t) - w(t - 1)
+```
+
+We can rewrite this more concisely using macros. Note the difference between
+including and excluding empty parentheses. When in call format, variables are
+aliased to the given call, allowing implicit use of dependents for convenience.
+
+```julia
+@parameters t() α() σ
+@variables w x(t) y() z(t, α, x)
+
+expr = x + y^α + σ(3) * (z - t) - w(t - 1)
+```
+
 
 ### Constants
 
@@ -243,37 +274,37 @@ is accessible via a function-based interface. This means that all macros are
 syntactic sugar in some form. For example, the variable construction:
 
 ```julia
-@parameters t σ ρ β
+@parameters t() σ ρ() β()
 @variables x(t) y(t) z(t)
 @derivatives D'~t
 ```
 
 is syntactic sugar for:
 
 ```julia
-t = Variable(:t; known = true)
-x = Variable(:x, [t])
-y = Variable(:y, [t])
-z = Variable(:z, [t])
-D = Differential(t)
+t = Variable(:t; known = true)()
 σ = Variable(:σ; known = true)
-ρ = Variable(:ρ; known = true)
-β = Variable(:β; known = true)
+ρ = Variable(:ρ; known = true)()
+β = Variable(:β; known = true)()
+x = Variable(:x)(t)
+y = Variable(:y)(t)
+z = Variable(:z)(t)
+D = Differential(t)
 ```
 
 ### Intermediate Calculations
 
 The system building functions can handle intermediate calculations. For example,
 
 ```julia
-@variables x y z
-@parameters σ ρ β
+@variables x() y() z()
+@parameters σ() ρ() β()
 a = y - x
 eqs = [0 ~ σ*a,
        0 ~ x*(ρ-z)-y,
        0 ~ x*y - β*z]
-ns = NonlinearSystem(eqs,[x,y,z],[σ,ρ,β])
-nlsys_func = generate_function(ns)
+ns = NonlinearSystem(eqs, [x,y,z])
+nlsys_func = generate_function(ns, [x,y,z], [σ,ρ,β])
 ```
 
 expands to:

src/ModelingToolkit.jl

@@ -2,6 +2,7 @@ module ModelingToolkit
 
 export Operation, Expression
 export calculate_jacobian, generate_jacobian, generate_function
+export independent_variables, dependent_variables, parameters
 export @register
 
 
@@ -22,6 +23,10 @@ function calculate_jacobian end
 function generate_jacobian end
 function generate_function end
 
+function independent_variables end
+function dependent_variables   end
+function parameters            end
+
 @enum FunctionVersion ArrayFunction=1 SArrayFunction=2
 
 include("variables.jl")

src/differentials.jl

@@ -4,27 +4,27 @@ export Differential, expand_derivatives, @derivatives
 struct Differential <: Function
     x::Expression
 end
+(D::Differential)(x) = Operation(D, Expression[x])
 
 Base.show(io::IO, D::Differential) = print(io, "(D'~", D.x, ")")
 Base.convert(::Type{Expr}, D::Differential) = D
 
-(D::Differential)(x::Operation) = Operation(D, Expression[x])
-function (D::Differential)(x::Variable)
-    D.x === x             && return Constant(1)
-    has_dependent(x, D.x) || return Constant(0)
-    return Operation(D, Expression[x])
-end
-(::Differential)(::Any) = Constant(0)
 Base.:(==)(D1::Differential, D2::Differential) = isequal(D1.x, D2.x)
 
 function expand_derivatives(O::Operation)
     @. O.args = expand_derivatives(O.args)
 
-    if O.op isa Differential
-        D = O.op
-        o = O.args[1]
-        isa(o, Operation) || return O
-        return simplify_constants(sum(i->derivative(o,i)*expand_derivatives(D(o.args[i])),1:length(o.args)))
+    if isa(O.op, Differential)
+        (D, o) = (O.op, O.args[1])
+
+        isequal(o, D.x)     && return Constant(1)
+        occursin(D.x, o)    || return Constant(0)
+        isa(o, Operation)   || return O
+        isa(o.op, Variable) && return O
+
+        return sum(1:length(o.args)) do i
+            derivative(o, i) * expand_derivatives(D(o.args[i]))
+        end |> simplify_constants
     end
 
     return O
@@ -80,6 +80,6 @@ macro derivatives(x...)
     esc(_differential_macro(x))
 end
 
-function calculate_jacobian(eqs,vars)
-    Expression[Differential(vars[j])(eqs[i]) for i in 1:length(eqs), j in 1:length(vars)]
+function calculate_jacobian(eqs, dvs)
+    Expression[Differential(dv)(eq) for eq ∈ eqs, dv ∈ dvs]
 end

src/equations.jl

@@ -10,39 +10,3 @@ Base.:(==)(a::Equation, b::Equation) = isequal((a.lhs, a.rhs), (b.lhs, b.rhs))
 Base.:~(lhs::Expression, rhs::Expression) = Equation(lhs, rhs)
 Base.:~(lhs::Expression, rhs::Number    ) = Equation(lhs, rhs)
 Base.:~(lhs::Number    , rhs::Expression) = Equation(lhs, rhs)
-
-
-_is_dependent(x::Variable) = !x.known && !isempty(x.dependents)
-_is_parameter(iv) = x -> x.known && !isequal(x, iv)
-_is_known(x::Variable) = x.known
-_is_unknown(x::Variable) = !x.known
-
-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
-
-    return result
-end
-
-get_args(O::Operation) = O.args
-get_args(eq::Equation) = Expression[eq.lhs, eq.rhs]
-function vars!(vars, op)
-    for arg ∈ get_args(op)
-        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/operations.jl

@@ -3,10 +3,8 @@ struct Operation <: Expression
     args::Vector{Expression}
 end
 
-# Recursive ==
-function Base.isequal(x::Operation,y::Operation)
+Base.isequal(x::Operation,y::Operation) =
     x.op == y.op && length(x.args) == length(y.args) && all(isequal.(x.args,y.args))
-end
 Base.isequal(::Operation, ::Number   ) = false
 Base.isequal(::Number   , ::Operation) = false
 Base.isequal(::Operation, ::Variable ) = false