Redesign Differential

README.md

@@ -37,9 +37,9 @@ using ModelingToolkit
 Then we build the system:
 
 ```julia
-eqs = [D*x ~ σ*(y-x),
-       D*y ~ x*(ρ-z)-y,
-       D*z ~ x*y - β*z]
+eqs = [D(x) ~ σ*(y-x),
+       D(y) ~ x*(ρ-z)-y,
+       D(z) ~ x*y - β*z]
 ```
 
 Each operation builds an `Operation` type, and thus `eqs` is an array of
@@ -163,7 +163,7 @@ context-aware single variable of the IR. Its fields are described as follows:
   to set units or denote a variable as being of higher precision.
 - `subtype`: the main denotation of context. Variables within systems
   are grouped according to their `subtype`.
-- `diff`: the operator objects attached to the variable
+- `diff`: the `Differential` object representing the quantity the variable is differentiated with respect to, or `nothing`
 - `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.
@@ -182,19 +182,13 @@ context-aware single variable of the IR. Its fields are described as follows:
 ### Operations
 
 Operations are the basic composition of variables and puts together the pieces
-with a function. The operator `~` is a special operator which denotes equality
-between the arguments.
+with a function. The `~` function denotes equality between the arguments.
 
-### Operators
+### Differentials
 
-An operator is an object which modifies variables via `*`. It adds the operator
-to the `diff` field of the variable and changes the interpretation of the variable.
-The current operators are:
-
-- `Differential`: a differential denotes the derivative with respect to a given
-  variable. It can be expanded via `expand_derivatives` which symbolically
-  differentiates expressions recursively and cancels out appropriate constant
-  variables.
+A `Differential` denotes the derivative with respect to a given variable. It can
+be expanded via `expand_derivatives`, which symbolically differentiates
+expressions recursively and cancels out appropriate constant variables.
 
 ### Systems
 
@@ -255,7 +249,7 @@ to better scale to larger systems. You can define derivatives for your own
 function via the dispatch:
 
 ```julia
-ModelingToolkit.Derivative(::typeof(my_function),args,::Type{Val{i}})
+ModelingToolkit.Derivative(::typeof(my_function),args,::Val{i})
 ```
 
 where `i` means that it's the derivative of the `i`th argument. `args` is the
@@ -265,7 +259,7 @@ You should return an `Operation` for the derivative of your function.
 For example, `sin(t)`'s derivative (by `t`) is given by the following:
 
 ```julia
-ModelingToolkit.Derivative(::typeof(sin),args,::Type{Val{1}}) = cos(args[1])
+ModelingToolkit.Derivative(::typeof(sin),args,::Val{1}) = cos(args[1])
 ```
 
 ### Macro-free Usage

src/ModelingToolkit.jl

@@ -8,7 +8,6 @@ import MacroTools: splitdef, combinedef
 
 abstract type Expression <: Number end
 abstract type AbstractOperation <: Expression end
-abstract type AbstractOperator <: Expression end
 abstract type AbstractComponent <: Expression end
 abstract type AbstractSystem end
 abstract type AbstractDomain end
@@ -26,14 +25,14 @@ function caclulate_jacobian end
 @enum FunctionVersions ArrayFunction=1 SArrayFunction=2
 
 include("operations.jl")
-include("operators.jl")
+include("differentials.jl")
 include("systems/diffeqs/diffeqsystem.jl")
 include("systems/diffeqs/first_order_transform.jl")
 include("systems/nonlinear/nonlinear_system.jl")
 include("function_registration.jl")
 include("simplify.jl")
 include("utils.jl")
 
-export Operation, Expression, AbstractOperator, AbstractComponent, AbstractDomain
+export Operation, Expression, AbstractComponent, AbstractDomain
 export @register
 end # module

src/differentials.jl

@@ -0,0 +1,91 @@
+struct Differential <: Function
+    x::Expression
+    order::Int
+end
+Differential(x) = Differential(x,1)
+
+Base.show(io::IO, D::Differential) = print(io,"($(D.x),$(D.order))")
+Base.Expr(D::Differential) = D
+
+function Derivative end
+(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 Variable(x,D)
+end
+Base.:(==)(D1::Differential, D2::Differential) = D1.order == D2.order && D1.x == D2.x
+
+Variable(x::Variable, D::Differential) = Variable(x.name,x.value,x.value_type,
+                x.subtype,D,x.dependents,x.description,x.flow,x.domain,
+                x.size,x.context)
+
+function expand_derivatives(O::Operation)
+    @. O.args = expand_derivatives(O.args)
+
+    if O.op isa Differential
+        D = O.op
+        o = O.args[1]
+        return simplify_constants(sum(i->Derivative(o,i)*expand_derivatives(D(o.args[i])),1:length(o.args)))
+    end
+
+    return O
+end
+expand_derivatives(x::Variable) = x
+
+# Don't specialize on the function here
+function Derivative(O::Operation,idx)
+    # This calls the Derivative dispatch from the user or pre-defined code
+    Derivative(O.op, O.args, Val(idx))
+end
+Derivative(op, args, idx) = Derivative(op, (args...,), idx)
+
+# Pre-defined derivatives
+import DiffRules, SpecialFunctions, NaNMath
+for (modu, fun, arity) ∈ DiffRules.diffrules()
+    for i ∈ 1:arity
+        @eval function Derivative(::typeof($modu.$fun), args::NTuple{$arity,Any}, ::Val{$i})
+            M, f = $(modu, fun)
+            partials = DiffRules.diffrule(M, f, args...)
+            dx = @static $arity == 1 ? partials : partials[$i]
+            parse(Operation,dx)
+        end
+    end
+end
+
+function count_order(x)
+    @assert !(x isa Symbol) "The variable $x must have an order of differentiation that is greater or equal to 1!"
+    n = 1
+    while !(x.args[1] isa Symbol)
+        n = n+1
+        x = x.args[1]
+    end
+    n, x.args[1]
+end
+
+function _differential_macro(x)
+    ex = Expr(:block)
+    lhss = Symbol[]
+    x = flatten_expr!(x)
+    for di in x
+        @assert di isa Expr && di.args[1] == :~ "@Deriv expects a form that looks like `@Deriv D''~t E'~t`"
+        lhs = di.args[2]
+        rhs = di.args[3]
+        order, lhs = count_order(lhs)
+        push!(lhss, lhs)
+        expr = :($lhs = Differential($rhs, $order))
+        push!(ex.args,  expr)
+    end
+    push!(ex.args, Expr(:tuple, lhss...))
+    ex
+end
+
+macro Deriv(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)]
+end
+
+export Differential, Derivative, expand_derivatives, @Deriv, calculate_jacobian

src/systems/diffeqs/first_order_transform.jl

@@ -36,7 +36,7 @@ function ode_order_lowering!(eqs, naming_scheme)
     for sym in keys(sym_order)
         order = sym_order[sym]
         for o in (order-1):-1:1
-            lhs = D*lower_varname(sym, idv, o-1, dv_name, naming_scheme)
+            lhs = D(lower_varname(sym, idv, o-1, dv_name, naming_scheme))
             rhs = lower_varname(sym, idv, o, dv_name, naming_scheme)
             eq = Operation(==, [lhs, rhs])
             push!(eqs, eq)
@@ -46,7 +46,7 @@ function ode_order_lowering!(eqs, naming_scheme)
 end
 
 function lhs_renaming!(eq, D, naming_scheme)
-    eq.args[1] = D*lower_varname(eq.args[1], naming_scheme, lower=true)
+    eq.args[1] = D(lower_varname(eq.args[1], naming_scheme, lower=true))
     return eq
 end
 function rhs_renaming!(eq, naming_scheme)

src/utils.jl

@@ -18,3 +18,7 @@ function flatten_expr!(x)
 end
 
 toexpr(ex) = MacroTools.postwalk(x->x isa Union{Expression,Operation} ? Expr(x) : x, ex)
+
+has_dependent(t::Variable) = Base.Fix2(has_dependent, t)
+has_dependent(x::Variable, t::Variable) =
+    t ∈ x.dependents || any(has_dependent(t), x.dependents)

src/variables.jl

@@ -3,7 +3,7 @@ mutable struct Variable <: Expression
     value
     value_type::DataType
     subtype::Symbol
-    diff::Union{AbstractOperator,Nothing}
+    diff::Union{Function,Nothing}  # FIXME
     dependents::Vector{Variable}
     description::String
     flow::Bool

test/basic_variables_and_operations.jl

@@ -17,8 +17,8 @@ s = JumpVariable(:s,3,dependents=[t])
 n = NoiseVariable(:n,dependents=[t])
 
 σ*(y-x)
-D*x
-D*x ~ -σ*(y-x)
-D*y ~ x*(ρ-z)-sin(y)
+D(x)
+D(x) ~ -σ*(y-x)
+D(y) ~ x*(ρ-z)-sin(y)
 
-@test D*t == Constant(1)
+@test D(t) == Constant(1)

test/components.jl

@@ -12,9 +12,9 @@ struct Lorenz <: AbstractComponent
 end
 function generate_lorenz_eqs(t,x,y,z,σ,ρ,β)
   D = Differential(t)
-  [D*x ~ σ*(y-x)
-   D*y ~ x*(ρ-z)-y
-   D*z ~ x*y - β*z]
+  [D(x) ~ σ*(y-x)
+   D(y) ~ x*(ρ-z)-y
+   D(z) ~ x*y - β*z]
 end
 Lorenz(t) = Lorenz(first(@DVar(x(t))),first(@DVar(y(t))),first(@DVar(z(t))),first(@Param(σ)),first(@Param(ρ)),first(@Param(β)),generate_lorenz_eqs(t,x,y,z,σ,ρ,β))
 

test/derivatives.jl

@@ -6,24 +6,24 @@ using Test
 @Var x(t) y(t) z(t)
 @Param σ ρ β
 @Deriv D'~t
-dsin = D*sin(t)
+dsin = D(sin(t))
 expand_derivatives(dsin)
 
 @test expand_derivatives(dsin) == cos(t)
-dcsch = D*csch(t)
+dcsch = D(csch(t))
 @test expand_derivatives(dcsch) == simplify_constants(Operation(coth(t)*csch(t)*-1))
 
 # Chain rule
-dsinsin = D*sin(sin(t))
+dsinsin = D(sin(sin(t)))
 @test expand_derivatives(dsinsin) == cos(sin(t))*cos(t)
 # Binary
-dpow1 = Derivative(^,[x, y],Val{1})
-dpow2 = Derivative(^,[x, y],Val{2})
+dpow1 = Derivative(^,[x, y],Val(1))
+dpow2 = Derivative(^,[x, y],Val(2))
 @test dpow1 == y*x^(y-1)
 @test dpow2 == x^y*log(x)
 
-d1 = D*(sin(t)*t)
-d2 = D*(sin(t)*cos(t))
+d1 = D(sin(t)*t)
+d2 = D(sin(t)*cos(t))
 @test expand_derivatives(d1) == t*cos(t)+sin(t)
 @test expand_derivatives(d2) == simplify_constants(cos(t)*cos(t)+sin(t)*(-1*sin(t)))
 
@@ -41,3 +41,7 @@ jac = ModelingToolkit.calculate_jacobian(sys)
 @test jac[3,1] == y
 @test jac[3,2] == x
 @test jac[3,3] == -1*β
+
+# Variable dependence checking in differentiation
+@Var a(t) b(a)
+@test D(b) ≠ Constant(0)