Merge pull request #740 from SciML/myb/dep

Deprecate `@derivatives` and better `derivatives` function

Author: ChrisRackauckas

Project.toml

@@ -1,7 +1,7 @@
 name = "ModelingToolkit"
 uuid = "961ee093-0014-501f-94e3-6117800e7a78"
 authors = ["Chris Rackauckas <[email protected]>"]
-version = "4.5.0"
+version = "5.0.0"
 
 [deps]
 ArrayInterface = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9"

README.md

@@ -30,7 +30,7 @@ using ModelingToolkit, OrdinaryDiffEq
 
 @parameters t σ ρ β
 @variables x(t) y(t) z(t)
-@derivatives D'~t
+D = Differential(t)
 
 eqs = [D(D(x)) ~ σ*(y-x),
        D(y) ~ x*(ρ-z)-y,
@@ -67,7 +67,7 @@ using ModelingToolkit, OrdinaryDiffEq
 
 @parameters t σ ρ β
 @variables x(t) y(t) z(t)
-@derivatives D'~t
+D = Differential(t)
 
 eqs = [D(x) ~ σ*(y-x),
        D(y) ~ x*(ρ-z)-y,

docs/src/IR.md

@@ -3,17 +3,19 @@
 ModelingToolkit IR mirrors the Julia AST but allows for easy mathematical
 manipulation by itself following mathematical semantics. The base of the IR is
 the `Sym` type, which defines a symbolic variable. Registered (mathematical)
-functions on `Sym`s (or `Term`s) return `Term`s.  For example, `op1 = x+y` is
-one `Term` and `op2 = 2z` is another, and so `op1*op2` is another `Term`. Then,
-at the top, an `Equation`, normally written as `op1 ~ op2`, defines the
-symbolic equality between two operations.
+functions on `Sym`s (or `istree` objects) return an expression that `istree`.
+For example, `op1 = x+y` is one symbolic object and `op2 = 2z` is another, and
+so `op1*op2` is another tree object. Then, at the top, an `Equation`, normally
+written as `op1 ~ op2`, defines the symbolic equality between two operations.
 
 ### Types
 `Sym`, `Term`, and `FnType` are from [SymbolicUtils.jl](https://juliasymbolics.github.io/SymbolicUtils.jl/api/). Note that in
 ModelingToolkit, we always use `Sym{Real}`, `Term{Real}`, and
-`FnType{Tuple{Any}, Real}`. To get the arguments of a `Term` use
-`arguments(t::Term)`, and to get the operation of a `Term` use
-`operation(t::Term)`.
+`FnType{Tuple{Any}, Real}`. To get the arguments of a `istree` object use
+`arguments(t::Term)`, and to get the operation, use `operation(t::Term)`.
+However, note that one should never dispatch on `Term` or test `isa Term`.
+Instead, one needs to use `SymbolicUtils.istree` to check if `arguments` and
+`operation` is defined.
 
 ```@docs
 Equation

docs/src/highlevel.md

@@ -10,14 +10,15 @@ methods.
 ```@docs
 @parameters
 @variables
-@derivatives
+Differential
 Base.:~(::Num, ::Num)
 modelingtoolkitize
 ```
 
 ## Differentiation Functions
 
 ```@docs
+ModelingToolkit.derivative
 ModelingToolkit.gradient
 ModelingToolkit.jacobian
 ModelingToolkit.sparsejacobian

docs/src/tutorials/higher_order.md

@@ -15,7 +15,7 @@ using ModelingToolkit, OrdinaryDiffEq
 
 @parameters t σ ρ β
 @variables x(t) y(t) z(t)
-@derivatives D'~t
+D = Differential(t)
 
 eqs = [D(D(x)) ~ σ*(y-x),
        D(y) ~ x*(ρ-z)-y,
@@ -25,10 +25,11 @@ sys = ODESystem(eqs)
 ```
 
 Note that we could've used an alternative syntax for 2nd order, i.e.
-`@derivatives E''~t` and then `E(x)` would be the second derivative,
-and this syntax extends to Nth order.
+`D = Differential(t)^2` and then `E(x)` would be the second derivative,
+and this syntax extends to `N`-th order. Also, we can use `*` or `∘` to compose
+`Differential`s, like `Differential(t) * Differential(x)`.
 
-Now let's transform this into the ODESystem of first order components.
+Now let's transform this into the `ODESystem` of first order components.
 We do this by simply calling `ode_order_lowering`:
 
 ```julia

docs/src/tutorials/ode_modeling.md

@@ -9,7 +9,7 @@ using ModelingToolkit, OrdinaryDiffEq
 
 @parameters t σ ρ β
 @variables x(t) y(t) z(t)
-@derivatives D'~t
+D = Differential(t)
 
 eqs = [D(x) ~ σ*(y-x),
        D(y) ~ x*(ρ-z)-y,
@@ -58,7 +58,7 @@ using ModelingToolkit
 
 @parameters t σ ρ β
 @variables x(t) y(t) z(t)
-@derivatives D'~t
+D = Differential(t)
 ```
 
 Then we build the system:

docs/src/tutorials/symbolic_functions.md

@@ -7,17 +7,15 @@ The way to define symbolic variables is via the `@variables` macro:
 @variables x y
 ```
 
-After defining variables as symbolic, symbolic expressions, which we
-call a `Term`, can be generated by utilizing Julia expressions.
-For example:
+After defining variables as symbolic, symbolic expressions, which we call a
+`istree` object, can be generated by utilizing Julia expressions. For example:
 
 ```julia
 z = x^2 + y
 ```
 
-Here, `z` is the `Term` for "square `x` and add `y`". To
-make an array of symbolic expressions, simply make an array of
-symbolic expressions:
+Here, `z` is an expression tree for "square `x` and add `y`". To make an array
+of symbolic expressions, simply make an array of symbolic expressions:
 
 ```julia
 A = [x^2+y 0 2x
@@ -30,7 +28,10 @@ A = [x^2+y 0 2x
   y ^ 2 + x  0  0
 ```
 
-Note that by default, `@variables` returns `Sym` or `Term` objects wrapped in `Num` in order to make them behave like subtypes of `Real`. Any operation on these `Num` objects will return a new `Num` object, wrapping the result of computing symbolically on the underlying values.
+Note that by default, `@variables` returns `Sym` or `istree` objects wrapped in
+`Num` in order to make them behave like subtypes of `Real`. Any operation on
+these `Num` objects will return a new `Num` object, wrapping the result of
+computing symbolically on the underlying values.
 
 To better view the results, we can use [Latexify.jl](https://github.com/korsbo/Latexify.jl).
 ModelingToolkit.jl comes with Latexify recipes so it works automatically:
@@ -282,11 +283,11 @@ sufficiently large!)
 One common thing to compute in a symbolic system is derivatives. In
 ModelingToolkit.jl, derivatives are represented lazily via operations,
 just like any other function. To build a differential operator, use
-`@derivatives` like:
+`Differential` like:
 
 ```julia
 @variables t
-@derivatives D'~t
+Differential(t)
 ```
 
 This is the differential operator ``D = \frac{\partial}{\partial t}``: the number of
@@ -450,7 +451,7 @@ z = g(x) + g(y)
 
 One of the benefits of a one-language Julia symbolic stack is that the
 primitives are all written in Julia, and therefore it's trivially
-extendable from Julia itself. By default, new functions are traced
+extendible from Julia itself. By default, new functions are traced
 to the primitives and the symbolic expressions are written on the
 primitives. However, we can expand the allowed primitives by registering
 new functions. For example, let's register a new function `h`:
@@ -479,7 +480,8 @@ ModelingToolkit.derivative(::typeof(h), args::NTuple{2,Any}, ::Val{2}) = 1
 and now it works with the rest of the system:
 
 ```julia
-@derivatives Dx'~x Dy'~y
+Dx = Differential(x)
+Dy = Differential(y)
 expand_derivatives(Dx(h(x,y) + y^2)) # 2x
 expand_derivatives(Dy(h(x,y) + y^2)) # 1 + 2y
 ```

src/ModelingToolkit.jl

@@ -23,7 +23,9 @@ RuntimeGeneratedFunctions.init(@__MODULE__)
 using RecursiveArrayTools
 
 import SymbolicUtils
-import SymbolicUtils: Term, Add, Mul, Pow, Sym, to_symbolic, FnType, @rule, Rewriters, substitute, similarterm
+import SymbolicUtils: Term, Add, Mul, Pow, Sym, to_symbolic, FnType,
+                      @rule, Rewriters, substitute, similarterm,
+                      promote_symtype
 
 import SymbolicUtils.Rewriters: Chain, Postwalk, Prewalk, Fixpoint
 

src/differentials.jl

@@ -16,10 +16,16 @@ julia> using ModelingToolkit
 julia> @variables x y;
 
 julia> D = Differential(x)
-(D'~x())
+(D'~x)
 
-julia> D(y)  # Differentiate y wrt. x
-(D'~x())(y())
+julia> D(y) # Differentiate y wrt. x
+(D'~x)(y)
+
+julia> Dx = Differential(x) * Differential(y) # d^2/dxy operator
+(D'~x(t)) ∘ (D'~y(t))
+
+julia> D3 = Differential(x)^3 # 3rd order differential operator
+(D'~x(t)) ∘ (D'~x(t)) ∘ (D'~x(t))
 ```
 """
 struct Differential <: Function
@@ -31,6 +37,11 @@ end
 (D::Differential)(x::Num) = Num(D(value(x)))
 SymbolicUtils.promote_symtype(::Differential, x) = x
 
+Base.:*(D1, D2::Differential) = D1 ∘ D2
+Base.:*(D1::Differential, D2) = D1 ∘ D2
+Base.:*(D1::Differential, D2::Differential) = D1 ∘ D2
+Base.:^(D::Differential, n::Integer) = _repeat_apply(D, n)
+
 Base.show(io::IO, D::Differential) = print(io, "(D'~", D.x, ")")
 
 Base.:(==)(D1::Differential, D2::Differential) = isequal(D1.x, D2.x)
@@ -229,6 +240,7 @@ end
 _repeat_apply(f, n) = n == 1 ? f : f ∘ _repeat_apply(f, n-1)
 function _differential_macro(x)
     ex = Expr(:block)
+    push!(ex.args,  :(Base.depwarn("`@derivatives D'''~x` is deprecated. Use `Differential(x)^3` instead.", Symbol("@derivatives"), force=true)))
     lhss = Symbol[]
     x = x isa Tuple && first(x).head == :tuple ? first(x).args : x # tuple handling
     x = flatten_expr!(x)

src/direct.jl

@@ -1,3 +1,19 @@
+"""
+```julia
+derivative(O, v; simplify = true)
+```
+
+A helper function for computing the derivative of an expression with respect to
+`var`.
+"""
+function derivative(O, v; simplify = true)
+    if O isa AbstractArray
+        Num[Num(expand_derivatives(Differential(v)(value(o)), simplify)) for o in O]
+    else
+        Num(expand_derivatives(Differential(v)(value(O)), simplify))
+    end
+end
+
 """
 ```julia
 gradient(O, vars::AbstractVector; simplify = true)
@@ -7,7 +23,7 @@ A helper function for computing the gradient of an expression with respect to
 an array of variable expressions.
 """
 function gradient(O, vars::AbstractVector; simplify = true)
-    Num[expand_derivatives(Differential(v)(value(O)),simplify) for v in vars]
+    Num[Num(expand_derivatives(Differential(v)(value(O)),simplify)) for v in vars]
 end
 
 """
@@ -19,7 +35,7 @@ A helper function for computing the Jacobian of an array of expressions with res
 an array of variable expressions.
 """
 function jacobian(ops::AbstractVector, vars::AbstractVector; simplify = true)
-    Num[expand_derivatives(Differential(value(v))(value(O)),simplify) for O in ops, v in vars]
+    Num[Num(expand_derivatives(Differential(value(v))(value(O)),simplify)) for O in ops, v in vars]
 end
 
 """
@@ -41,7 +57,7 @@ function sparsejacobian(ops::AbstractVector, vars::AbstractVector; simplify = tr
     exprs = Num[]
 
     for (i,j) in zip(I, J)
-        push!(exprs, expand_derivatives(Differential(vars[j])(ops[i]), simplify))
+        push!(exprs, Num(expand_derivatives(Differential(vars[j])(ops[i]), simplify)))
     end
     sparse(I,J, exprs, length(ops), length(vars))
 end
@@ -232,7 +248,9 @@ function toexpr(O; canonicalize=true)
     op = operation(O)
     args = arguments(O)
     if op isa Differential
-        return :(derivative($(toexpr(args[1]; canonicalize=canonicalize)),$(toexpr(op.x; canonicalize=canonicalize))))
+        ex = toexpr(args[1]; canonicalize=canonicalize)
+        wrt = toexpr(op.x; canonicalize=canonicalize)
+        return :(_derivative($ex, $wrt))
     elseif op isa Sym
         isempty(args) && return nameof(op)
         return Expr(:call, toexpr(op; canonicalize=canonicalize), toexpr(args; canonicalize=canonicalize)...)