finish optimizationsystem

Author: ChrisRackauckas

src/ModelingToolkit.jl

@@ -99,7 +99,7 @@ include("systems/diffeqs/validation.jl")
 
 include("systems/nonlinear/nonlinearsystem.jl")
 
-#include("systems/optimization/optimizationsystem.jl")
+include("systems/optimization/optimizationsystem.jl")
 
 include("systems/pde/pdesystem.jl")
 
@@ -110,7 +110,7 @@ include("build_function.jl")
 
 export ODESystem, ODEFunction
 export SDESystem, SDEFunction
-export NonlinearSystem
+export NonlinearSystem, OptimizationSystem
 export ode_order_lowering
 export PDESystem
 export Reaction, ReactionSystem
@@ -119,12 +119,14 @@ export IntervalDomain, ProductDomain, ⊗, CircleDomain
 export Equation, ConstrainedEquation
 export simplify_constants
 
-export Operation, Expression
+export Operation, Expression, Variable
 export calculate_jacobian, generate_jacobian, generate_function
+export calculate_tgrad, generate_tgrad
+export calculate_hessian, generate_hessian
 export calculate_massmatrix, generate_diffusion_function
 export independent_variable, states, parameters, equations
-export simplified_expr, eval_function
-export @register, @I
+export simplified_expr
+export @register
 export modelingtoolkitize
-export Variable, @variables, @parameters
+export @variables, @parameters
 end # module

src/build_function.jl

@@ -8,6 +8,42 @@ function build_function(args...;target = JuliaTarget(),kwargs...)
   _build_function(target,args...;kwargs...)
 end
 
+# Scalar output
+function _build_function(target::JuliaTarget, op::Operation, vs, ps = (), args = (),
+                         conv = simplified_expr, expression = Val{true};
+                         checkbounds = false, constructor=nothing,
+                         linenumbers = true)
+    _vs = map(x-> x isa Operation ? x.op : x, vs)
+    _ps = map(x-> x isa Operation ? x.op : x, ps)
+    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))
+
+    fname = gensym(:ModelingToolkitFunction)
+    out_expr = conv(op)
+    let_expr = Expr(:let, var_eqs, out_expr)
+    bounds_block = checkbounds ? let_expr : :(@inbounds begin $let_expr end)
+
+    fargs = ps == () ? :(u,$(args...)) : :(u,p,$(args...))
+
+    oop_ex = :(
+        ($(fargs.args...),) -> begin
+            $bounds_block
+        end
+    )
+
+    if !linenumbers
+        oop_ex = striplines(oop_ex)
+    end
+
+    if expression == Val{true}
+        return oop_ex
+    else
+        return GeneralizedGenerated.mk_function(@__MODULE__,oop_ex)
+    end
+end
+
 function _build_function(target::JuliaTarget, rhss, vs, ps = (), args = (),
                          conv = simplified_expr, expression = Val{true};
                          checkbounds = false, constructor=nothing,

src/systems/abstractsystem.jl

@@ -10,7 +10,7 @@ function Base.getproperty(sys::AbstractSystem, name::Symbol)
     i = findfirst(x->x.name==name,sys.states)
     if i !== nothing
         x = rename(sys.states[i],renamespace(sys.name,name))
-        if iv ∈ fieldnames(typeof(sys))
+        if :iv ∈ fieldnames(typeof(sys))
             return x(getfield(sys,:iv)())
         else
             return x()
@@ -22,7 +22,7 @@ function Base.getproperty(sys::AbstractSystem, name::Symbol)
             return rename(sys.ps[i],renamespace(sys.name,name))()
         end
     end
-    throw(error("Variable name does not exist"))
+    throw(error("Variable $name does not exist"))
 end
 
 renamespace(namespace,name) = Symbol(string(namespace)*"′"*string(name))

src/systems/diffeqs/abstractodesystem.jl

@@ -46,7 +46,7 @@ end
 
 function generate_function(sys::AbstractODESystem, dvs = states(sys), ps = parameters(sys), expression = Val{true}; kwargs...)
     rhss = [deq.rhs for deq ∈ equations(sys)]
-    dvs′ = convert.(Variable,states)
+    dvs′ = convert.(Variable,dvs)
     ps′ = convert.(Variable,ps)
     return build_function(rhss, dvs′, ps′, (sys.iv.name,), ODEToExpr(sys), expression; kwargs...)
 end

src/systems/diffeqs/odesystem.jl

@@ -64,9 +64,9 @@ end
 function ODESystem(deqs::AbstractVector{<:Equation}, iv, dvs, ps;
                    systems = ODESystem[],
                    name=gensym(:ODESystem))
-    iv′ = clean(iv)
-    dvs′ = [clean(dv) for dv ∈ dvs]
-    ps′ = [clean(p) for p ∈ ps]
+    iv′ = convert(Variable,iv)
+    dvs′ = convert.(Variable,dvs)
+    ps′ = convert.(Variable,ps)
     tgrad = RefValue(Vector{Expression}(undef, 0))
     jac = RefValue(Matrix{Expression}(undef, 0, 0))
     Wfact   = RefValue(Matrix{Expression}(undef, 0, 0))
@@ -92,7 +92,7 @@ end
 
 Base.:(==)(sys1::ODESystem, sys2::ODESystem) =
     _eq_unordered(sys1.eqs, sys2.eqs) && isequal(sys1.iv, sys2.iv) &&
-    _eq_unordered(sys1.dvs, sys2.dvs) && _eq_unordered(sys1.ps, sys2.ps)
+    _eq_unordered(sys1.states, sys2.states) && _eq_unordered(sys1.ps, sys2.ps)
 # NOTE: equality does not check cached Jacobian
 
 function rename(sys::ODESystem,name)

src/systems/diffeqs/sdesystem.jl

@@ -42,9 +42,9 @@ end
 function SDESystem(deqs::AbstractVector{<:Equation}, neqs, iv, dvs, ps;
                    systems = SDESystem[],
                    name = gensym(:SDESystem))
-    dvs′ = [clean(dv) for dv ∈ dvs]
-    ps′ = [clean(p) for p ∈ ps]
-    iv′ = clean(iv)
+    iv′ = convert(Variable,iv)
+    dvs′ = convert.(Variable,dvs)
+    ps′ = convert.(Variable,ps)
     tgrad = RefValue(Vector{Expression}(undef, 0))
     jac = RefValue(Matrix{Expression}(undef, 0, 0))
     Wfact   = RefValue(Matrix{Expression}(undef, 0, 0))
@@ -53,8 +53,8 @@ function SDESystem(deqs::AbstractVector{<:Equation}, neqs, iv, dvs, ps;
 end
 
 function generate_diffusion_function(sys::SDESystem, dvs = sys.states, ps = sys.ps, expression = Val{true}; kwargs...)
-    dvs′ = [clean(dv) for dv ∈ dvs]
-    ps′ = [clean(p) for p ∈ ps]
+    dvs′ = convert.(Variable,dvs)
+    ps′ = convert.(Variable,ps)
     return build_function(sys.noiseeqs, dvs′, ps′, (sys.iv.name,), ODEToExpr(sys), expression; kwargs...)
 end
 

src/systems/nonlinear/nonlinearsystem.jl

@@ -22,7 +22,7 @@ struct NonlinearSystem <: AbstractSystem
     """Vector of equations defining the system."""
     eqs::Vector{Equation}
     """Unknown variables."""
-    states::Vector{Expression}
+    states::Vector{Variable}
     """Parameters."""
     ps::Vector{Variable}
     """
@@ -38,16 +38,16 @@ end
 function NonlinearSystem(eqs, states, ps;
                          name = gensym(:NonlinearSystem),
                          systems = NonlinearSystem[])
-    NonlinearSystem(eqs, states, convert.(Variable,ps), name, systems)
+    NonlinearSystem(eqs, convert.(Variable,states), convert.(Variable,ps), name, systems)
 end
 
 function calculate_jacobian(sys::NonlinearSystem)
     rhs = [eq.rhs for eq in sys.eqs]
-    jac = expand_derivatives.(calculate_jacobian(rhs, sys.states))
+    jac = expand_derivatives.(calculate_jacobian(rhs, [dv() for dv in states(sys)]))
     return jac
 end
 
-function generate_jacobian(sys::NonlinearSystem, vs = sys.states, ps = sys.ps, expression = Val{true}; kwargs...)
+function generate_jacobian(sys::NonlinearSystem, vs = states(sys), ps = parameters(sys), expression = Val{true}; kwargs...)
     jac = calculate_jacobian(sys)
     return build_function(jac, convert.(Variable,vs), convert.(Variable,ps), (), NLSysToExpr(sys))
 end
@@ -65,7 +65,7 @@ function (f::NLSysToExpr)(O::Operation)
 end
 (f::NLSysToExpr)(x) = convert(Expr, x)
 
-function generate_function(sys::NonlinearSystem, vs = sys.states, ps = sys.ps, expression = Val{true}; kwargs...)
+function generate_function(sys::NonlinearSystem, vs = states(sys), ps = parameters(sys), expression = Val{true}; kwargs...)
     rhss = [eq.rhs for eq ∈ sys.eqs]
     vs′ = convert.(Variable,vs)
     ps′ = convert.(Variable,ps)

src/systems/optimization/optimizationsystem.jl

@@ -20,7 +20,7 @@ struct OptimizationSystem <: AbstractSystem
     """Vector of equations defining the system."""
     op::Operation
     """Unknown variables."""
-    states::Vector{Expression}
+    states::Vector{Variable}
     """Parameters."""
     ps::Vector{Variable}
     """
@@ -30,42 +30,29 @@ struct OptimizationSystem <: AbstractSystem
     """
     systems: The internal systems
     """
-    systems::Vector{NonlinearSystem}
+    systems::Vector{OptimizationSystem}
 end
 
-function NonlinearSystem(eqs, states, ps;
-                         name = gensym(:NonlinearSystem),
-                         systems = NonlinearSystem[])
-    NonlinearSystem(eqs, states, convert.(Variable,ps), name, systems)
+function OptimizationSystem(op, states, ps;
+                            name = gensym(:OptimizationSystem),
+                            systems = OptimizationSystem[])
+    OptimizationSystem(op, convert.(Variable,states), convert.(Variable,ps), name, systems)
 end
 
-function calculate_jacobian(sys::NonlinearSystem)
-    rhs = [eq.rhs for eq in sys.eqs]
-    jac = expand_derivatives.(calculate_jacobian(rhs, sys.states))
-    return jac
+function calculate_hessian(sys::OptimizationSystem)
+    expand_derivatives.(hessian(equations(sys), [dv() for dv in states(sys)]))
 end
 
-function generate_jacobian(sys::NonlinearSystem, vs = sys.states, ps = sys.ps, expression = Val{true}; kwargs...)
-    jac = calculate_jacobian(sys)
-    return build_function(jac, convert.(Variable,vs), convert.(Variable,ps), (), NLSysToExpr(sys))
+function generate_hessian(sys::OptimizationSystem, vs = states(sys), ps = parameters(sys), expression = Val{true}; kwargs...)
+    hes = calculate_hessian(sys)
+    return build_function(hes, convert.(Variable,vs), convert.(Variable,ps), (), x->convert(Expr, x))
 end
 
-struct NLSysToExpr
-    sys::NonlinearSystem
-end
-function (f::NLSysToExpr)(O::Operation)
-    any(isequal(O), f.sys.states) && return O.op.name  # variables
-    if isa(O.op, Variable)
-        isempty(O.args) && return O.op.name  # 0-ary parameters
-        return build_expr(:call, Any[O.op.name; f.(O.args)])
-    end
-    return build_expr(:call, Any[O.op; f.(O.args)])
-end
-(f::NLSysToExpr)(x) = convert(Expr, x)
-
-function generate_function(sys::NonlinearSystem, vs = sys.states, ps = sys.ps, expression = Val{true}; kwargs...)
-    rhss = [eq.rhs for eq ∈ sys.eqs]
+function generate_function(sys::OptimizationSystem, vs = states(sys), ps = parameters(sys), expression = Val{true}; kwargs...)
     vs′ = convert.(Variable,vs)
     ps′ = convert.(Variable,ps)
-    return build_function(rhss, vs′, ps′, (), NLSysToExpr(sys), expression; kwargs...)
+    return build_function(equations(sys), vs′, ps′, (), x->convert(Expr, x), expression; kwargs...)
 end
+
+equations(sys::OptimizationSystem) = isempty(sys.systems) ? sys.op : sys.op + reduce(+,namespace_operation.(sys.systems))
+namespace_operation(sys::OptimizationSystem) = namespace_operation(sys.op,sys.name,nothing)

test/build_targets.jl

@@ -4,7 +4,7 @@ using ModelingToolkit, Test
 @derivatives D'~t
 eqs = [D(x) ~ a*x - x*y,
        D(y) ~ -3y + x*y]
-@test ModelingToolkit.build_function(eqs,ModelingToolkit.clean.([x,y]),ModelingToolkit.clean.([a]),t,target = ModelingToolkit.StanTarget()) ==
+@test ModelingToolkit.build_function(eqs,convert.(Variable,[x,y]),convert.(Variable,[a]),t,target = ModelingToolkit.StanTarget()) ==
     """
     real[] diffeqf(real t,real[] internal_var___u,real[] internal_var___p,real[] x_r,int[] x_i) {
       real internal_var___du[2];
@@ -38,10 +38,10 @@ sys = ODESystem(eqs,t,[x,y],[a])
       ModelingToolkit.build_function(sys.eqs,[x,y],[a],t,target = ModelingToolkit.MATLABTarget())
 
 @test ModelingToolkit.build_function(eqs,[x,y],[a],t,target = ModelingToolkit.CTarget()) ==
-      ModelingToolkit.build_function(sys.eqs,sys.dvs,sys.ps,sys.iv,target = ModelingToolkit.CTarget())
+      ModelingToolkit.build_function(sys.eqs,sys.states,sys.ps,sys.iv,target = ModelingToolkit.CTarget())
 
 @test ModelingToolkit.build_function(eqs,[x,y],[a],t,target = ModelingToolkit.StanTarget()) ==
-      ModelingToolkit.build_function(sys.eqs,sys.dvs,sys.ps,sys.iv,target = ModelingToolkit.StanTarget())
+      ModelingToolkit.build_function(sys.eqs,sys.states,sys.ps,sys.iv,target = ModelingToolkit.StanTarget())
 
 @test ModelingToolkit.build_function(eqs,[x,y],[a],t,target = ModelingToolkit.MATLABTarget()) ==
-      ModelingToolkit.build_function(sys.eqs,sys.dvs,sys.ps,sys.iv,target = ModelingToolkit.MATLABTarget())
+      ModelingToolkit.build_function(sys.eqs,sys.states,sys.ps,sys.iv,target = ModelingToolkit.MATLABTarget())

test/derivatives.jl

@@ -3,7 +3,7 @@ using Test
 
 # Derivatives
 @parameters t σ ρ β
-@variables x(t) y(t) z(t)
+@variables x y z
 @derivatives D'~t D2''~t Dx'~x
 
 @test @macroexpand(@derivatives D'~t D2''~t) == @macroexpand(@derivatives (D'~t), (D2''~t))
@@ -36,7 +36,7 @@ d2 = D(sin(t)*cos(t))
 eqs = [0 ~ σ*(y-x),
        0 ~ x*(ρ-z)-y,
        0 ~ x*y - β*z]
-sys = NonlinearSystem(eqs, [x,y,z])
+sys = NonlinearSystem(eqs, [x,y,z], [σ,ρ,β])
 jac = calculate_jacobian(sys)
 @test isequal(jac[1,1], σ*-1)
 @test isequal(jac[1,2], σ)
@@ -54,6 +54,8 @@ jac = calculate_jacobian(sys)
 @test isequal(expand_derivatives(D(t)), 1)
 @test isequal(expand_derivatives(Dx(x)), 1)
 
+@variables x(t) y(t) z(t)
+
 @test isequal(expand_derivatives(D(x * y)), simplify_constants(y*D(x) + x*D(y)))
 @test_broken isequal(expand_derivatives(D(x * y)), simplify_constants(D(x)*y + x*D(y)))