Merge pull request #1681 from SciML/s/calling-convention-update

Use new dependent array variable convention

Author: shashi

Project.toml

@@ -73,7 +73,7 @@ Setfield = "0.7, 0.8, 1"
 SpecialFunctions = "0.7, 0.8, 0.9, 0.10, 1.0, 2"
 StaticArrays = "0.10, 0.11, 0.12, 1.0"
 SymbolicUtils = "0.19"
-Symbolics = "4.5 - 4.8"
+Symbolics = "4.9"
 UnPack = "0.1, 1.0"
 Unitful = "1.1"
 julia = "1.6"

docs/src/tutorials/spring_mass.md

@@ -13,14 +13,14 @@ D = Differential(t)
 
 function Mass(; name, m = 1.0, xy = [0., 0.], u = [0., 0.])
     ps = @parameters m=m
-    sts = @variables pos[1:2](t)=xy v[1:2](t)=u
+    sts = @variables pos(t)[1:2]=xy v(t)[1:2]=u
     eqs = scalarize(D.(pos) .~ v)
     ODESystem(eqs, t, [pos..., v...], ps; name)
 end
 
 function Spring(; name, k = 1e4, l = 1.)
     ps = @parameters k=k l=l
-    @variables x(t), dir[1:2](t)
+    @variables x(t), dir(t)[1:2]
     ODESystem(Equation[], t, [x, dir...], ps; name)
 end
 
@@ -63,7 +63,7 @@ For each component we use a Julia function that returns an `ODESystem`. At the t
 ```@example component
 function Mass(; name, m = 1.0, xy = [0., 0.], u = [0., 0.])
     ps = @parameters m=m
-    sts = @variables pos[1:2](t)=xy v[1:2](t)=u
+    sts = @variables pos(t)[1:2]=xy v(t)[1:2]=u
     eqs = scalarize(D.(pos) .~ v)
     ODESystem(eqs, t, [pos..., v...], ps; name)
 end
@@ -86,7 +86,7 @@ Next we build the spring component. It is characterised by the spring constant `
 ```@example component
 function Spring(; name, k = 1e4, l = 1.)
     ps = @parameters k=k l=l
-    @variables x(t), dir[1:2](t)
+    @variables x(t), dir(t)[1:2]
     ODESystem(Equation[], t, [x, dir...], ps; name)
 end
 ```

src/parameters.jl

@@ -1,9 +1,19 @@
-import SymbolicUtils: symtype, term, hasmetadata
+import SymbolicUtils: symtype, term, hasmetadata, issym
 struct MTKParameterCtx end
 
-isparameter(x::Num) = isparameter(value(x))
-isparameter(x::Symbolic) = getmetadata(x, MTKParameterCtx, false)
-isparameter(x) = false
+function isparameter(x)
+    x = unwrap(x)
+    if istree(x) && operation(x) isa Symbolic
+        getmetadata(x, MTKParameterCtx, false) ||
+        isparameter(operation(x))
+    elseif istree(x) && operation(x) == (getindex)
+        isparameter(arguments(x)[1])
+    elseif x isa Symbolic
+        getmetadata(x, MTKParameterCtx, false)
+    else
+        false
+    end
+end
 
 """
     toparam(s::Sym)
@@ -15,13 +25,11 @@ function toparam(s)
         Symbolics.wrap(toparam(Symbolics.unwrap(s)))
     elseif s isa AbstractArray
         map(toparam, s)
-    elseif symtype(s) <: AbstractArray
-        Symbolics.recurse_and_apply(toparam, s)
     else
         setmetadata(s, MTKParameterCtx, true)
     end
 end
-toparam(s::Num) = Num(toparam(value(s)))
+toparam(s::Num) = wrap(toparam(value(s)))
 
 """
     tovar(s::Sym)

src/systems/connectors.jl

@@ -282,9 +282,7 @@ function generate_connection_equations_and_stream_connections(csets::AbstractVec
 
     for cset in csets
         v = cset.set[1].v
-        if hasmetadata(v, Symbolics.GetindexParent)
-            v = getparent(v)
-        end
+        v = getparent(v, v)
         vtype = get_connection_type(v)
         if vtype === Stream
             push!(stream_connections, cset)

src/utils.jl

@@ -148,7 +148,7 @@ function check_variables(dvs, iv)
     for dv in dvs
         isequal(iv, dv) &&
             throw(ArgumentError("Independent variable $iv not allowed in dependent variables."))
-        (is_delay_var(iv, dv) || occursin(iv, iv_from_nested_derivative(dv))) ||
+        (is_delay_var(iv, dv) || occursin(iv, dv)) ||
             throw(ArgumentError("Variable $dv is not a function of independent variable $iv."))
     end
 end
@@ -201,19 +201,35 @@ function check_equations(eqs, iv)
     end
 end
 "Get all the independent variables with respect to which differentials/differences are taken."
-function collect_ivs_from_nested_operator!(ivs, x::Term, target_op)
-    op = operation(x)
+function collect_ivs_from_nested_operator!(ivs, x, target_op)
+    if !istree(x)
+        return
+    end
+    op = operation(unwrap(x))
     if op isa target_op
         push!(ivs, get_iv(op))
-        collect_ivs_from_nested_operator!(ivs, arguments(x)[1], target_op)
+        x = if target_op <: Differential
+            op.x
+        elseif target_op <: Difference
+            op.t
+        else
+            error("Unknown target op type in collect_ivs $target_op. Pass Difference or Differential")
+        end
+        collect_ivs_from_nested_operator!(ivs, x, target_op)
     end
 end
 
-function iv_from_nested_derivative(x::Term, op = Differential)
-    operation(x) isa op ? iv_from_nested_derivative(arguments(x)[1], op) : arguments(x)[1]
+function iv_from_nested_derivative(x, op = Differential)
+    if istree(x) && operation(x) == getindex
+        iv_from_nested_derivative(arguments(x)[1], op)
+    elseif istree(x)
+        operation(x) isa op ? iv_from_nested_derivative(arguments(x)[1], op) : arguments(x)[1]
+    elseif issym(x)
+        x
+    else
+        nothing
+    end
 end
-iv_from_nested_derivative(x::Sym, op = Differential) = x
-iv_from_nested_derivative(x, op = Differential) = nothing
 
 hasdefault(v) = hasmetadata(v, Symbolics.VariableDefaultValue)
 getdefault(v) = value(getmetadata(v, Symbolics.VariableDefaultValue))

test/latexify.jl

@@ -30,7 +30,7 @@ eqs = [D(x) ~ σ * (y - x) * D(x - y) / D(z),
 # Latexify.@generate_test latexify(eqs)
 @test_reference "latexify/10.tex" latexify(eqs)
 
-@variables u[1:3](t)
+@variables u(t)[1:3]
 @parameters p[1:3]
 eqs = [D(u[1]) ~ p[3] * (u[2] - u[1]),
     0 ~ p[2] * p[3] * u[1] * (p[1] - u[1]) / 10 - u[2],

test/latexify/20.tex

@@ -1,5 +1,5 @@
 \begin{align}
-\frac{du_1(t)}{dt} =& \left(  - \mathrm{u_1}\left( t \right) + \mathrm{u_2}\left( t \right) \right) p_3 \\
-0 =&  - \mathrm{u_2}\left( t \right) + \frac{1}{10} \left(  - \mathrm{u_1}\left( t \right) + p_1 \right) \mathrm{u_1}\left( t \right) p_2 p_3 \\
-\frac{du_3(t)}{dt} =& \left( \mathrm{u_2}\left( t \right) \right)^{\frac{2}{3}} \mathrm{u_1}\left( t \right) - \mathrm{u_3}\left( t \right) p_3
+\mathrm{\frac{d}{d t}}\left( u(t)_1 \right) =& \left(  - u(t)_1 + u(t)_2 \right) p_3 \\
+0 =&  - u(t)_2 + \frac{1}{10} \left(  - u(t)_1 + p_1 \right) p_2 p_3 u(t)_1 \\
+\mathrm{\frac{d}{d t}}\left( u(t)_3 \right) =& u(t)_2^{\frac{2}{3}} u(t)_1 - p_3 u(t)_3
 \end{align}

test/latexify/30.tex

@@ -1,5 +1,5 @@
 \begin{align}
-\frac{du_1(t)}{dt} =& \left(  - \mathrm{u_1}\left( t \right) + \mathrm{u_2}\left( t \right) \right) p_3 \\
-\frac{du_2(t)}{dt} =&  - \mathrm{u_2}\left( t \right) + \frac{1}{10} \left(  - \mathrm{u_1}\left( t \right) + p_1 \right) \mathrm{u_1}\left( t \right) p_2 p_3 \\
-\frac{du_3(t)}{dt} =& \left( \mathrm{u_2}\left( t \right) \right)^{\frac{2}{3}} \mathrm{u_1}\left( t \right) - \mathrm{u_3}\left( t \right) p_3
+\mathrm{\frac{d}{d t}}\left( u(t)_1 \right) =& \left(  - u(t)_1 + u(t)_2 \right) p_3 \\
+\mathrm{\frac{d}{d t}}\left( u(t)_2 \right) =&  - u(t)_2 + \frac{1}{10} \left(  - u(t)_1 + p_1 \right) p_2 p_3 u(t)_1 \\
+\mathrm{\frac{d}{d t}}\left( u(t)_3 \right) =& u(t)_2^{\frac{2}{3}} u(t)_1 - p_3 u(t)_3
 \end{align}

test/mass_matrix.jl

@@ -1,6 +1,6 @@
 using OrdinaryDiffEq, ModelingToolkit, Test, LinearAlgebra
 @parameters t
-@variables y[1:3](t)
+@variables y(t)[1:3]
 @parameters k[1:3]
 D = Differential(t)
 

test/odesystem.jl

@@ -376,7 +376,7 @@ end
 
 # issue 1109
 let
-    @variables t x[1:3, 1:3](t)
+    @variables t x(t)[1:3, 1:3]
     D = Differential(t)
     @named sys = ODESystem(D.(x) .~ x)
     @test_nowarn structural_simplify(sys)
@@ -386,7 +386,7 @@ end
 using Symbolics: unwrap, wrap
 using LinearAlgebra
 @variables t
-sts = @variables x[1:3](t)=[1, 2, 3.0] y(t)=1.0
+sts = @variables x(t)[1:3]=[1, 2, 3.0] y(t)=1.0
 ps = @parameters p[1:3] = [1, 2, 3]
 D = Differential(t)
 eqs = [collect(D.(x) .~ x)
@@ -487,7 +487,7 @@ function foo(a::Num, ms::AbstractVector)
     wrap(term(foo, a, term(SVector, ms...)))
 end
 foo(a, ms::AbstractVector) = a + sum(ms)
-@variables t x(t) ms[1:3](t)
+@variables t x(t) ms(t)[1:3]
 D = Differential(t)
 ms = collect(ms)
 eqs = [D(x) ~ foo(x, ms); D.(ms) .~ 1]
@@ -589,7 +589,7 @@ end
 let
     @parameters t
     D = Differential(t)
-    @variables x[1:2](t) = zeros(2)
+    @variables x(t)[1:2] = zeros(2)
     @variables y(t) = 0
     @parameters k = 1
     eqs = [D(x[1]) ~ x[2]
@@ -680,7 +680,7 @@ let
     eqs_to_lhs(eqs) = eq_to_lhs.(eqs)
 
     @parameters σ=10 ρ=28 β=8 / 3 sigma rho beta
-    @variables t t2 x(t)=1 y(t)=0 z(t)=0 x2(t2)=1 y2(t2)=0 z2(t2)=0 u[1:3](t2)
+    @variables t t2 x(t)=1 y(t)=0 z(t)=0 x2(t2)=1 y2(t2)=0 z2(t2)=0 u(t2)[1:3]
 
     D = Differential(t)
     D2 = Differential(t2)
@@ -749,7 +749,7 @@ end
 let
     @parameters t
 
-    u = collect(first(@variables u[1:4](t)))
+    u = collect(first(@variables u(t)[1:4]))
     Dt = Differential(t)
 
     eqs = [Differential(t)(u[2]) - 1.1u[1] ~ 0