Merge 2547a4e into 79b8772

README.md

@@ -136,8 +136,8 @@ context-aware single variable of the IR. Its fields are described as follows:
 - `name`: the name of the `Variable`. Note that this is not necessarily
   the same as the name of the Julia variable. But this symbol itself is considered
   the core identifier of the `Variable` in the sense of equality.
-- `subtype`: the main denotation of context. Variables within systems
-  are grouped according to their `subtype`.
+- `known`: the main denotation of context, storing whether or not the value of
+  the variable is known.
 - `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.
@@ -252,9 +252,9 @@ is syntactic sugar for:
 
 ```julia
 t = Parameter(:t)
-x = Unknown(:x, dependents = [t])
-y = Unknown(:y, dependents = [t])
-z = Unknown(:z, dependents = [t])
+x = Unknown(:x, [t])
+y = Unknown(:y, [t])
+z = Unknown(:z, [t])
 D = Differential(t)
 σ = Parameter(:σ)
 ρ = Parameter(:ρ)

src/ModelingToolkit.jl

@@ -7,31 +7,32 @@ using MacroTools
 import MacroTools: splitdef, combinedef
 
 abstract type Expression <: Number end
-abstract type AbstractOperation <: Expression end
-abstract type AbstractComponent <: Expression end
+abstract type AbstractSystem end
 
-include("variables.jl")
-
-Base.promote_rule(::Type{T},::Type{T2}) where {T<:Number,T2<:Expression} = Expression
+Base.promote_rule(::Type{<:Number},::Type{<:Expression}) = Expression
 Base.zero(::Type{<:Expression}) = Constant(0)
 Base.one(::Type{<:Expression}) = Constant(1)
-Base.convert(::Type{Variable},x::Int64) = Constant(x)
 
-function caclulate_jacobian end
+function calculate_jacobian end
+function generate_jacobian end
+function generate_function end
 
 @enum FunctionVersion ArrayFunction=1 SArrayFunction=2
 
+include("variables.jl")
 include("operations.jl")
 include("differentials.jl")
 include("equations.jl")
-include("systems/systems.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, AbstractComponent, AbstractDomain
+export Operation, Expression, AbstractComponent
+export calculate_jacobian, generate_jacobian, generate_function
+export ArrayFunction, SArrayFunction
 export @register
+
 end # module

src/equations.jl

@@ -5,17 +5,17 @@ struct Equation
     lhs::Expression
     rhs::Expression
 end
-Base.broadcastable(eq::Equation) = Ref(eq)
 Base.:(==)(a::Equation, b::Equation) = (a.lhs, a.rhs) == (b.lhs, b.rhs)
 
 Base.:~(lhs::Expression, rhs::Expression) = Equation(lhs, rhs)
 Base.:~(lhs::Expression, rhs::Number    ) = Equation(lhs, rhs)
 Base.:~(lhs::Number    , rhs::Expression) = Equation(lhs, rhs)
 
 
-_is_dependent(x::Variable) = x.subtype === :Unknown && !isempty(x.dependents)
-_is_parameter(iv) = x -> x.subtype === :Parameter && x ≠ iv
-_subtype(subtype::Symbol) = x -> x.subtype === subtype
+_is_dependent(x::Variable) = !x.known && !isempty(x.dependents)
+_is_parameter(iv) = x -> x.known && x ≠ iv
+_is_known(x::Variable) = x.known
+_is_unknown(x::Variable) = !x.known
 
 function extract_elements(eqs, predicates)
     result = [Variable[] for p ∈ predicates]

src/operations.jl

@@ -1,5 +1,4 @@
-# Parameterize by T so that way it can be Vector{Expression} which is defined after
-struct Operation <: AbstractOperation
+struct Operation <: Expression
     op::Function
     args::Vector{Expression}
 end
@@ -19,22 +18,6 @@ Base.convert(::Type{Expr}, O::Operation) =
     build_expr(:call, Any[Symbol(O.op); convert.(Expr, O.args)])
 Base.show(io::IO, O::Operation) = print(io, convert(Expr, O))
 
-
-"""
-find_replace(O::Operation, x::Expression, y::Expression)
-
-Finds the expression `x` in Operation `O` and replaces it with the Expression `y`
-"""
-function find_replace!(O::Operation, x::Expression, y::Expression)
-    for i in eachindex(O.args)
-        if isequal(O.args[i], x)
-            O.args[i] = y
-        elseif typeof(O.args[i]) <: Operation
-            find_replace!(O.args[i],x,y)
-        end
-    end
-end
-
 # For inv
 Base.convert(::Type{Operation}, x::Number) = Operation(identity, Expression[Constant(x)])
 Base.convert(::Type{Operation}, x::Operation) = x

src/systems/diffeqs/first_order_transform.jl

@@ -1,13 +1,11 @@
-extract_idv(eq::DiffEq) = eq.D.x
-
 function lower_varname(D::Differential, x; lower=false)
     order = lower ? D.order-1 : D.order
     return lower_varname(x, D.x, order)
 end
 function lower_varname(var::Variable, idv, order::Int)
     sym = var.name
     name = order == 0 ? sym : Symbol(sym, :_, string(idv.name)^order)
-    return Variable(name, var.subtype, var.dependents)
+    return Variable(name, var.known, var.dependents)
 end
 
 function ode_order_lowering(sys::DiffEqSystem)
@@ -21,7 +19,7 @@ function ode_order_lowering(eqs, iv)
     new_eqs = similar(eqs, DiffEq)
 
     for (i, eq) ∈ enumerate(eqs)
-        var, maxorder = extract_var_order(eq)
+        var, maxorder = eq.var, eq.D.order
         maxorder == 1 && continue # fast pass
         if maxorder > get(var_order, var, 0)
             var_order[var] = maxorder
@@ -51,6 +49,4 @@ function rename(O::Expression)
     return Operation(O.op, rename.(O.args))
 end
 
-extract_var_order(eq::DiffEq) = (eq.var, eq.D.order)
-
 export ode_order_lowering

src/systems/nonlinear/nonlinear_system.jl

@@ -18,7 +18,7 @@ struct NonlinearSystem <: AbstractSystem
 end
 
 function NonlinearSystem(eqs)
-    vs, ps = extract_elements(eqs, [_subtype(:Unknown), _subtype(:Parameter)])
+    vs, ps = extract_elements(eqs, [_is_unknown, _is_known])
     NonlinearSystem(eqs, vs, ps)
 end
 

src/utils.jl

@@ -17,7 +17,6 @@ function build_expr(head::Symbol, args)
     append!(ex.args, args)
     ex
 end
-expr_arr_to_block(exprs) = build_expr(:block, exprs)
 
 # used in parsing
 isblock(x) = length(x) == 1 && x[1] isa Expr && x[1].head == :block
@@ -31,11 +30,30 @@ function flatten_expr!(x)
     x
 end
 
-function partition(f, xs)
-    idxs = map(f, xs)
-    return (xs[idxs], xs[(!).(idxs)])
+function build_function(rhss, vs, ps, args = (); version::FunctionVersion)
+    var_pairs   = [(u.name, :(u[$i])) for (i, u) ∈ enumerate(vs)]
+    param_pairs = [(p.name, :(p[$i])) for (i, p) ∈ enumerate(ps)]
+    (ls, rs) = zip(var_pairs..., param_pairs...)
+
+    var_eqs = Expr(:(=), build_expr(:tuple, ls), build_expr(:tuple, rs))
+
+    if version === ArrayFunction
+        X = gensym()
+        sys_exprs = [:($X[$i] = $(convert(Expr, rhs))) for (i, rhs) ∈ enumerate(rhss)]
+        let_expr = Expr(:let, var_eqs, build_expr(:block, sys_exprs))
+        :(($X,u,p,$(args...)) -> $let_expr)
+    elseif version === SArrayFunction
+        sys_expr = build_expr(:tuple, [convert(Expr, rhs) for rhs ∈ rhss])
+        let_expr = Expr(:let, var_eqs, sys_expr)
+        :((u,p,$(args...)) -> begin
+            X = $let_expr
+            T = StaticArrays.similar_type(typeof(u), eltype(X))
+            T(X)
+        end)
+    end
 end
 
+
 is_constant(::Constant) = true
 is_constant(::Any) = false
 

src/variables.jl

@@ -3,12 +3,12 @@ export Variable, Unknown, Parameter, @Unknown, @Param
 
 struct Variable <: Expression
     name::Symbol
-    subtype::Symbol
+    known::Bool
     dependents::Vector{Variable}
 end
 
-Parameter(name; dependents = Variable[]) = Variable(name, :Parameter, dependents)
-Unknown(name; dependents = Variable[]) = Variable(name, :Unknown, dependents)
+Parameter(name, dependents = Variable[]) = Variable(name, true, dependents)
+Unknown(name, dependents = Variable[]) = Variable(name, false, dependents)
 
 
 struct Constant <: Expression
@@ -22,7 +22,7 @@ Base.isone(ex::Expression)  = isa(ex, Constant) && isone(ex.value)
 
 
 # Variables use isequal for equality since == is an Operation
-Base.:(==)(x::Variable, y::Variable) = (x.name, x.subtype) == (y.name, y.subtype)
+Base.:(==)(x::Variable, y::Variable) = (x.name, x.known) == (y.name, y.known)
 Base.:(==)(::Variable, ::Number) = false
 Base.:(==)(::Number, ::Variable) = false
 Base.:(==)(::Variable, ::Constant) = false
@@ -32,13 +32,18 @@ Base.:(==)(n::Number, c::Constant) = c.value == n
 Base.:(==)(a::Constant, b::Constant) = a.value == b.value
 
 function Base.convert(::Type{Expr}, x::Variable)
-    x.subtype === :Parameter || return x.name
-    isempty(x.dependents)    && return x.name
+    x.known               || return x.name
+    isempty(x.dependents) && return x.name
     return :($(x.name)($(convert.(Expr, x.dependents)...)))
 end
 Base.convert(::Type{Expr}, c::Constant) = c.value
 
-Base.show(io::IO, x::Variable) = print(io, x.subtype, '(', x.name, ')')
+function Base.show(io::IO, x::Variable)
+    subtype = x.known ? :Parameter : :Unknown
+    print(io, subtype, '(', repr(x.name))
+    isempty(x.dependents) || print(io, ", ", x.dependents)
+    print(io, ')')
+end
 
 # Build variables more easily
 function _parse_vars(macroname, fun, x)
@@ -65,7 +70,7 @@ function _parse_vars(macroname, fun, x)
         end
 
         push!(var_names, var_name)
-        expr = :($var_name = $fun($(Meta.quot(var_name)), dependents = $dependents))
+        expr = :($var_name = $fun($(Meta.quot(var_name)), $dependents))
         push!(ex.args, expr)
     end
     push!(ex.args, build_expr(:tuple, var_names))

test/derivatives.jl

@@ -4,11 +4,14 @@ using Test
 # Derivatives
 @Param t σ ρ β
 @Unknown x(t) y(t) z(t)
-@Deriv D'~t
-dsin = D(sin(t))
-expand_derivatives(dsin)
+@Deriv D'~t D2''~t
+
+@test expand_derivatives(D(t)) == 1
+@test expand_derivatives(D(D(t))) == 0
 
+dsin = D(sin(t))
 @test expand_derivatives(dsin) == cos(t)
+
 dcsch = D(csch(t))
 @test expand_derivatives(dcsch) == simplify_constants(coth(t) * csch(t) * -1)
 
@@ -25,7 +28,7 @@ eqs = [0 ~ σ*(y-x),
        0 ~ x*(ρ-z)-y,
        0 ~ x*y - β*z]
 sys = NonlinearSystem(eqs,[x,y,z],[σ,ρ,β])
-jac = ModelingToolkit.calculate_jacobian(sys)
+jac = calculate_jacobian(sys)
 @test jac[1,1] == σ*-1
 @test jac[1,2] == σ
 @test jac[1,3] == 0