Merge pull request #268 from SciML/auto_detect

Make variables typed and use types for auto-detection

Author: ChrisRackauckas

src/ModelingToolkit.jl

@@ -78,6 +78,7 @@ Get the set of parameters variables for the given system.
 function parameters end
 
 include("variables.jl")
+include("context_dsl.jl")
 include("operations.jl")
 include("differentials.jl")
 include("equations.jl")

src/build_function.jl

@@ -93,7 +93,7 @@ function _build_function(target::JuliaTarget, rhss, vs, ps = (), args = (),
 end
 
 get_varnumber(varop::Operation,vars::Vector{Operation}) =  findfirst(x->isequal(x,varop),vars)
-get_varnumber(varop::Operation,vars::Vector{Variable})  =  findfirst(x->isequal(x,varop.op),vars)
+get_varnumber(varop::Operation,vars::Vector{<:Variable})  =  findfirst(x->isequal(x,varop.op),vars)
 
 function numbered_expr(O::Equation,args...;kwargs...)
   :($(numbered_expr(O.lhs,args...;kwargs...)) = $(numbered_expr(O.rhs,args...;kwargs...)))
@@ -120,7 +120,7 @@ function numbered_expr(O::Operation,vars,parameters;
                         varname=varname,paramname=paramname) for x in O.args]...)
 end
 
-function numbered_expr(de::ModelingToolkit.Equation,vars::Vector{Variable},parameters;
+function numbered_expr(de::ModelingToolkit.Equation,vars::Vector{<:Variable},parameters;
                        derivname=:du,varname=:u,paramname=:p)
     i = findfirst(x->isequal(x.name,var_from_nested_derivative(de.lhs)[1].name),vars)
     :($derivname[$i] = $(numbered_expr(de.rhs,vars,parameters;

src/context_dsl.jl

@@ -0,0 +1,12 @@
+struct Parameter{T} end
+isparameter(::Variable) = false
+isparameter(::Variable{<:Parameter}) = true
+
+"""
+$(SIGNATURES)
+
+Define one or more known variables.
+"""
+macro parameters(xs...)
+    esc(_parse_vars(:parameters, Parameter{Number}, xs))
+end

src/systems/diffeqs/abstractodesystem.jl

@@ -55,7 +55,7 @@ function calculate_factorized_W(sys::AbstractODESystem, simplify=true)
     isempty(sys.Wfact[]) || return (sys.Wfact[],sys.Wfact_t[])
 
     jac = calculate_jacobian(sys)
-    gam = Variable(:gam; known = true)()
+    gam = Variable(:__MTKWgamma)()
 
     W = - LinearAlgebra.I + gam*jac
     Wfact = lu(W, Val(false), check=false).factors
@@ -76,15 +76,15 @@ function calculate_factorized_W(sys::AbstractODESystem, simplify=true)
 end
 
 function generate_factorized_W(sys::AbstractODESystem, vs = states(sys), ps = parameters(sys), simplify=true, expression = Val{true}; kwargs...)
-    (Wfact,Wfact_t) = calculate_factorized_W(sys,simplify)
+    Wfact,Wfact_t = calculate_factorized_W(sys,simplify)
     siz = size(Wfact)
     constructor = :(x -> begin
                         A = SMatrix{$siz...}(x)
                         StaticArrays.LU(LowerTriangular( SMatrix{$siz...}(UnitLowerTriangular(A)) ), UpperTriangular(A), SVector(ntuple(n->n, max($siz...))))
                     end)
 
-    Wfact_func   = build_function(Wfact  , vs, ps, (:gam,:t), ODEToExpr(sys), expression;constructor=constructor,kwargs...)
-    Wfact_t_func = build_function(Wfact_t, vs, ps, (:gam,:t), ODEToExpr(sys), expression;constructor=constructor,kwargs...)
+    Wfact_func   = build_function(Wfact  , vs, ps, (:__MTKWgamma,:t), ODEToExpr(sys), expression;constructor=constructor,kwargs...)
+    Wfact_t_func = build_function(Wfact_t, vs, ps, (:__MTKWgamma,:t), ODEToExpr(sys), expression;constructor=constructor,kwargs...)
 
     return (Wfact_func, Wfact_t_func)
 end

src/systems/diffeqs/first_order_transform.jl

@@ -1,7 +1,7 @@
 function lower_varname(var::Variable, idv, order)
     order == 0 && return var
     name = Symbol(var.name, :_, string(idv.name)^order)
-    return Variable(name; known = var.known)
+    return Variable{vartype(var)}(name)
 end
 
 function flatten_differential(O::Operation)

src/systems/diffeqs/modelingtoolkitize.jl

@@ -9,7 +9,7 @@ function modelingtoolkitize(prob::DiffEqBase.ODEProblem)
     @parameters t
     vars = [Variable(:x, i)(t) for i in eachindex(prob.u0)]
     params = prob.p isa DiffEqBase.NullParameters ? [] :
-             [Variable(:α,i; known = true)() for i in eachindex(prob.p)]
+             [Variable(:α,i)() for i in eachindex(prob.p)]
     @derivatives D'~t
 
     rhs = [D(var) for var in vars]
@@ -22,7 +22,7 @@ function modelingtoolkitize(prob::DiffEqBase.ODEProblem)
     end
 
     eqs = vcat([rhs[i] ~ lhs[i] for i in eachindex(prob.u0)]...)
-    de = ODESystem(eqs)
+    de = ODESystem(eqs,t,vars,params)
 
     de, vars, params
 end

src/systems/diffeqs/odesystem.jl

@@ -85,7 +85,7 @@ function ODESystem(eqs; kwargs...)
 
     dvs = unique(var_from_nested_derivative(eq.lhs)[1] for eq ∈ eqs)
     ps = filter(vars(eq.rhs for eq ∈ eqs)) do x
-        x.known & !isequal(x, iv)
+        isparameter(x) & !isequal(x, iv)
     end |> collect
     ODESystem(eqs, iv, dvs, ps; kwargs...)
 end

src/utils.jl

@@ -35,7 +35,7 @@ end
 
 function detime_dvs(op::Operation)
   if op.op isa Variable
-    Operation(Variable(op.op.name,known=op.op.known),Expression[])
+    Operation(Variable{vartype(op.op)}(op.op.name),Expression[])
   else
     Operation(op.op,detime_dvs.(op.args))
   end
@@ -44,7 +44,7 @@ detime_dvs(op::Constant) = op
 
 function retime_dvs(op::Operation,dvs,iv)
   if op.op isa Variable && op.op ∈ dvs
-    Operation(Variable(op.op.name),Expression[iv])
+    Operation(Variable{vartype(op.op)}(op.op.name),Expression[iv])
   else
     Operation(op.op,retime_dvs.(op.args,(dvs,),iv))
   end

src/variables.jl

@@ -17,34 +17,35 @@ end
 """
 $(TYPEDEF)
 
-A named variable which represents a numerical value. The variable's value may
-be known (parameters, independent variables) or unknown (dependent variables).
+A named variable which represents a numerical value.
 
 # Fields
 $(FIELDS)
 """
-struct Variable <: Function
+struct Variable{T} <: Function
     """The variable's unique name."""
     name::Symbol
-    """
-    Whether the variable's value is known.
-    """
-    known::Bool
-    Variable(name; known = false) = new(name, known)
+    Variable(name) = new{Number}(name)
+    Variable{T}(name) where T = new{T}(name)
+    function Variable{T}(name, indices...) where T
+        var_name = Symbol("$(name)$(join(map_subscripts.(indices), "ˏ"))")
+        Variable{T}(var_name)
+    end
 end
-function Variable(name, indices...; known = false)
+
+function Variable(name, indices...)
     var_name = Symbol("$(name)$(join(map_subscripts.(indices), "ˏ"))")
-    Variable(var_name; known=known)
+    Variable(var_name)
 end
 
+vartype(::Variable{T}) where T = T
 (x::Variable)(args...) = Operation(x, collect(Expression, args))
 
-Base.isequal(x::Variable, y::Variable) = (x.name, x.known) == (y.name, y.known)
+Base.isequal(x::Variable, y::Variable) = x.name == y.name
 Base.print(io::IO, x::Variable) = show(io, x)
 Base.show(io::IO, x::Variable) = print(io, x.name)
 function Base.show(io::IO, ::MIME"text/plain", x::Variable)
-    known = x.known ? "known" : "unknown"
-    print(io, x.name, " (callable ", known, " variable)")
+    print(io, x.name)
 end
 
 
@@ -79,7 +80,7 @@ Base.convert(::Type{Expr}, c::Constant) = c.value
 
 
 # Build variables more easily
-function _parse_vars(macroname, known, x)
+function _parse_vars(macroname, type, x)
     ex = Expr(:block)
     var_names = Symbol[]
     # if parsing things in the form of
@@ -97,9 +98,9 @@ function _parse_vars(macroname, known, x)
         @assert iscall || isarray || issym "@$macroname expects a tuple of expressions or an expression of a tuple (`@$macroname x y z(t) v[1:3] w[1:2,1:4]` or `@$macroname x, y, z(t) v[1:3] w[1:2,1:4]`)"
 
         if iscall
-            var_name, expr = _construct_vars(_var.args[1], known, _var.args[2:end])
+            var_name, expr = _construct_vars(_var.args[1], type, _var.args[2:end])
         else
-            var_name, expr = _construct_vars(_var, known, nothing)
+            var_name, expr = _construct_vars(_var, type, nothing)
         end
         push!(var_names, var_name)
         push!(ex.args, expr)
@@ -108,45 +109,45 @@ function _parse_vars(macroname, known, x)
     return ex
 end
 
-function _construct_vars(_var, known, call_args)
+function _construct_vars(_var, type, call_args)
     issym  = _var isa Symbol
     isarray = isa(_var, Expr) && _var.head == :ref
     if isarray
         var_name = _var.args[1]
         indices = _var.args[2:end]
-        expr = _construct_array_vars(var_name, known, call_args, indices...)
+        expr = _construct_array_vars(var_name, type, call_args, indices...)
     else
         # Implicit 0-args call
         var_name = _var
-        expr = _construct_var(var_name, known, call_args)
+        expr = _construct_var(var_name, type, call_args)
     end
     var_name, :($var_name = $expr)
 end
 
-function _construct_var(var_name, known, call_args)
+function _construct_var(var_name, type, call_args)
     if call_args === nothing
-        :(Variable($(Meta.quot(var_name)); known = $known)())
+        :(Variable{$type}($(Meta.quot(var_name)))())
     elseif !isempty(call_args) && call_args[end] == :..
-        :(Variable($(Meta.quot(var_name)); known = $known))
+        :(Variable{$type}($(Meta.quot(var_name))))
     else
-        :(Variable($(Meta.quot(var_name)); known = $known)($(call_args...)))
+        :(Variable{$type}($(Meta.quot(var_name)))($(call_args...)))
     end
 end
 
-function _construct_var(var_name, known, call_args, ind)
+function _construct_var(var_name, type, call_args, ind)
     if call_args === nothing
-        :(Variable($(Meta.quot(var_name)), $ind...; known = $known)())
+        :(Variable{$type}($(Meta.quot(var_name)), $ind...)())
     elseif !isempty(call_args) && call_args[end] == :..
-        :(Variable($(Meta.quot(var_name)), $ind...; known = $known))
+        :(Variable{$type}($(Meta.quot(var_name)), $ind...))
     else
-        :(Variable($(Meta.quot(var_name)), $ind...; known = $known)($(call_args...)))
+        :(Variable{$type}($(Meta.quot(var_name)), $ind...)($(call_args...)))
     end
 end
 
 
-function _construct_array_vars(var_name, known, call_args, indices...)
+function _construct_array_vars(var_name, type, call_args, indices...)
     :(map(Iterators.product($(indices...))) do ind
-        $(_construct_var(var_name, known, call_args, :ind))
+        $(_construct_var(var_name, type, call_args, :ind))
     end)
 end
 
@@ -157,20 +158,11 @@ $(SIGNATURES)
 Define one or more unknown variables.
 """
 macro variables(xs...)
-    esc(_parse_vars(:variables, false, xs))
-end
-
-"""
-$(SIGNATURES)
-
-Define one or more known variables.
-"""
-macro parameters(xs...)
-    esc(_parse_vars(:parameters, true, xs))
+    esc(_parse_vars(:variables, Number, xs))
 end
 
 function rename(x::Variable,name::Symbol)
-    Variable(name,known=x.known)
+    Variable{vartype(x)}(name)
 end
 
 TreeViews.hastreeview(x::Variable) = true

test/variable_parsing.jl

@@ -1,6 +1,8 @@
 using ModelingToolkit
 using Test
 
+const PType = ModelingToolkit.Parameter{Number}
+
 @parameters t
 @variables x(t) y(t) # test multi-arg
 @variables z(t) # test single-arg
@@ -17,9 +19,9 @@ z1 = Variable(:z)(t)
 end
 @parameters σ(..)
 
-t1 = Variable(:t; known = true)()
-s1 = Variable(:s; known = true)()
-σ1 = Variable(:σ; known = true)
+t1 = Variable{PType}(:t)()
+s1 = Variable{PType}(:s)()
+σ1 = Variable(:σ)
 @test isequal(t1, t)
 @test isequal(s1, s)
 @test isequal(σ1, σ)
@@ -42,12 +44,12 @@ convert(Expression, :($x == 0 ? $y : $x))
 end
 @parameters σ[1:2](..)
 
-t1 = [Variable(:t, 1; known = true)(),
-      Variable(:t, 2; known = true)()]
-s1 = [Variable(:s, 1, 1; known = true)() Variable(:s, 1, 2; known = true)()
-      Variable(:s, 3, 1; known = true)() Variable(:s, 3, 2; known = true)()]
-σ1 = [Variable(:σ, 1; known = true),
-      Variable(:σ, 2; known = true)]
+t1 = [Variable{PType}(:t, 1)(),
+      Variable{PType}(:t, 2)()]
+s1 = [Variable{PType}(:s, 1, 1)() Variable{PType}(:s, 1, 2)()
+      Variable{PType}(:s, 3, 1)() Variable{PType}(:s, 3, 2)()]
+σ1 = [Variable(:σ, 1),
+      Variable(:σ, 2)]
 @test isequal(t1, t)
 @test isequal(s1, s)
 @test isequal(σ1, σ)