Merge branch 'master' into hg/feature/varfuns

Author: HarrisonGrodin

REQUIRE

@@ -1,4 +1,4 @@
-julia 1.0
+julia 1.1
 MacroTools
 DiffEqBase
 DiffRules

src/differentials.jl

@@ -61,9 +61,10 @@ _repeat_apply(f, n) = n == 1 ? f : f ∘ _repeat_apply(f, n-1)
 function _differential_macro(x)
     ex = Expr(:block)
     lhss = Symbol[]
+    x = x isa Tuple && first(x).head == :tuple ? first(x).args : x # tuple handling
     x = flatten_expr!(x)
     for di in x
-        @assert di isa Expr && di.args[1] == :~ "@derivatives expects a form that looks like `@derivatives D''~t E'~t`"
+        @assert di isa Expr && di.args[1] == :~ "@derivatives expects a form that looks like `@derivatives D''~t E'~t` or `@derivatives (D''~t), (E'~t)`"
         lhs = di.args[2]
         rhs = di.args[3]
         order, lhs = count_order(lhs)

src/systems/diffeqs/diffeqsystem.jl

@@ -99,10 +99,10 @@ function generate_jacobian(sys::ODESystem; version::FunctionVersion = ArrayFunct
     return build_function(jac, sys.dvs, sys.ps, (sys.iv.name,); version = version)
 end
 
-struct DiffEqToExpr
+struct ODEToExpr
     sys::ODESystem
 end
-function (f::DiffEqToExpr)(O::Operation)
+function (f::ODEToExpr)(O::Operation)
     if isa(O.op, Variable)
         isequal(O.op, f.sys.iv) && return O.op.name  # independent variable
         O.op ∈ f.sys.dvs        && return O.op.name  # dependent variables
@@ -111,41 +111,39 @@ function (f::DiffEqToExpr)(O::Operation)
     end
     return build_expr(:call, Any[O.op; f.(O.args)])
 end
-(f::DiffEqToExpr)(x) = convert(Expr, x)
+(f::ODEToExpr)(x) = convert(Expr, x)
 
 function generate_function(sys::ODESystem, dvs, ps; version::FunctionVersion = ArrayFunction)
     rhss = [deq.rhs for deq ∈ sys.eqs]
     dvs′ = [clean(dv) for dv ∈ dvs]
     ps′ = [clean(p) for p ∈ ps]
-    return build_function(rhss, dvs′, ps′, (sys.iv.name,), DiffEqToExpr(sys); version = version)
+    return build_function(rhss, dvs′, ps′, (sys.iv.name,), ODEToExpr(sys); version = version)
 end
 
 
-function generate_ode_iW(sys::ODESystem, simplify=true; version::FunctionVersion = ArrayFunction)
+function generate_factorized_W(sys::ODESystem, simplify=true; version::FunctionVersion = ArrayFunction)
     jac = calculate_jacobian(sys)
 
     gam = Variable(:gam; known = true)()
 
     W = LinearAlgebra.I - gam*jac
-    W = SMatrix{size(W,1),size(W,2)}(W)
-    iW = inv(W)
+    Wfact = lu(W, Val(false), check=false).factors
 
     if simplify
-        iW = simplify_constants.(iW)
+        Wfact = simplify_constants.(Wfact)
     end
 
-    W = inv(LinearAlgebra.I/gam - jac)
-    W = SMatrix{size(W,1),size(W,2)}(W)
-    iW_t = inv(W)
+    W_t = LinearAlgebra.I/gam - jac
+    Wfact_t = lu(W_t, Val(false), check=false).factors
     if simplify
-        iW_t = simplify_constants.(iW_t)
+        Wfact_t = simplify_constants.(Wfact_t)
     end
 
     vs, ps = sys.dvs, sys.ps
-    iW_func   = build_function(iW  , vs, ps, (:gam,:t); version = version)
-    iW_t_func = build_function(iW_t, vs, ps, (:gam,:t); version = version)
+    Wfact_func   = build_function(Wfact  , vs, ps, (:gam,:t), ODEToExpr(sys); version = version)
+    Wfact_t_func = build_function(Wfact_t, vs, ps, (:gam,:t), ODEToExpr(sys); version = version)
 
-    return (iW_func, iW_t_func)
+    return (Wfact_func, Wfact_t_func)
 end
 
 function DiffEqBase.ODEFunction(sys::ODESystem, dvs, ps; version::FunctionVersion = ArrayFunction)

src/variables.jl

@@ -38,11 +38,12 @@ function _parse_vars(macroname, known, x)
     #     y
     #     z
     # end
+    x = x isa Tuple && first(x) isa Expr && first(x).head == :tuple ? first(x).args : x # tuple handling
     x = flatten_expr!(x)
     for _var in x
         iscall = isa(_var, Expr) && _var.head == :call
         issym    = _var isa Symbol
-        @assert iscall || issym "@$macroname expects a tuple of expressions (`@$macroname x y z(t)`)"
+        @assert iscall || issym "@$macroname expects a tuple of expressions or an expression of a tuple (`@$macroname x y z(t)` or `@$macroname x, y, z(t)`)"
 
         if iscall
             var_name = _var.args[1]
@@ -59,8 +60,8 @@ function _parse_vars(macroname, known, x)
     return ex
 end
 macro variables(xs...)
-    esc(_parse_vars(:Variable, false, xs))
+    esc(_parse_vars(:variables, false, xs))
 end
 macro parameters(xs...)
-    esc(_parse_vars(:Param, true, xs))
+    esc(_parse_vars(:parameters, true, xs))
 end

test/derivatives.jl

@@ -6,6 +6,8 @@ using Test
 @variables x(t) y(t) z(t)
 @derivatives D'~t D2''~t Dx'~x
 
+@test @macroexpand(@derivatives D'~t D2''~t) == @macroexpand(@derivatives (D'~t), (D2''~t))
+
 @test isequal(expand_derivatives(D(t)), 1)
 @test isequal(expand_derivatives(D(D(t))), 0)
 

test/system_construction.jl

@@ -36,7 +36,15 @@ generate_function(de, [x,y,z], [σ,ρ,β]; version=ModelingToolkit.SArrayFunctio
 jac_expr = generate_jacobian(de)
 jac = calculate_jacobian(de)
 f = ODEFunction(de, [x,y,z], [σ,ρ,β])
-ModelingToolkit.generate_ode_iW(de)
+fw, fwt = map(eval, ModelingToolkit.generate_factorized_W(de))
+du = zeros(3)
+u  = collect(1:3)
+p  = collect(4:6)
+f(du, u, p, 0.1)
+@test du == [4, 0, -16]
+FW = zeros(3, 3)
+fw(FW, u, p, 0.2, 0.1)
+fwt(FW, u, p, 0.2, 0.1)
 
 @testset "time-varying parameters" begin
     @parameters σ′(t-1)
@@ -90,6 +98,10 @@ lowered_eqs = [D(u_tt) ~ 2u_tt + u_t + x_t + 1
                D(u)    ~ u_t
                D(x)    ~ x_t]
 @test de1 == ODESystem(lowered_eqs)
+test_diffeq_inference("first-order transform", de1, t, [u_tt, x_t, u_t, u, x], [])
+du = zeros(5)
+ODEFunction(de1, [u_tt, x_t, u_t, u, x], [])(du, ones(5), nothing, 0.1)
+@test du == [5.0, 3.0, 1.0, 1.0, 1.0]
 
 # Internal calculations
 a = y - x

test/variable_parsing.jl

@@ -2,7 +2,8 @@ using ModelingToolkit
 using Test
 
 @parameters t()
-@variables x(t) y(t) z(t)
+@variables x(t) y(t) # test multi-arg
+@variables z(t) # test single-arg
 x1 = Variable(:x)(t)
 y1 = Variable(:y)(t)
 z1 = Variable(:z)(t)
@@ -27,3 +28,6 @@ D1 = Differential(t)
 @test D1 == D
 
 @test isequal(x ≤ y + 1, (x < y + 1) | (x == y + 1))
+
+@test @macroexpand(@parameters x, y, z(t)) == @macroexpand(@parameters x y z(t))
+@test @macroexpand(@variables x, y, z(t)) == @macroexpand(@variables x y z(t))