Merge pull request #788 from SciML/myb/set

ChrisRackauckas · web-flow · commit 0176fef5ad40 · 2021-02-12T17:30:00.000-05:00
Depwarn fixes
diff --git a/Project.toml b/Project.toml
@@ -1,10 +1,11 @@
 name = "ModelingToolkit"
 uuid = "961ee093-0014-501f-94e3-6117800e7a78"
 authors = ["Chris Rackauckas <accounts@chrisrackauckas.com>"]
-version = "5.6.1"
+version = "5.6.2"
 
 [deps]
 ArrayInterface = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9"
+ConstructionBase = "187b0558-2788-49d3-abe0-74a17ed4e7c9"
 DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
 DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e"
 DiffEqJump = "c894b116-72e5-5b58-be3c-e6d8d4ac2b12"
@@ -37,6 +38,7 @@ Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 
 [compat]
 ArrayInterface = "2.8, 3.0"
+ConstructionBase = "1"
 DataStructures = "0.17, 0.18"
 DiffEqBase = "6.54.0"
 DiffEqJump = "6.7.5"
diff --git a/src/ModelingToolkit.jl b/src/ModelingToolkit.jl
@@ -6,7 +6,7 @@ using StaticArrays, LinearAlgebra, SparseArrays, LabelledArrays
 using Latexify, Unitful, ArrayInterface
 using MacroTools
 using UnPack: @unpack
-using Setfield
+using Setfield, ConstructionBase
 using DiffEqJump
 using DataStructures
 using SpecialFunctions, NaNMath
diff --git a/src/systems/abstractsystem.jl b/src/systems/abstractsystem.jl
@@ -149,6 +149,11 @@ for prop in [
              :Wfact_t
              :systems
              :structure
+             :op
+             :equality_constraints
+             :inequality_constraints
+             :controls
+             :loss
             ]
     fname1 = Symbol(:get_, prop)
     fname2 = Symbol(:has_, prop)
@@ -158,6 +163,25 @@ for prop in [
     end
 end
 
+Setfield.get(obj::AbstractSystem, l::Setfield.PropertyLens{field}) where {field} = getfield(obj, field)
+@generated function ConstructionBase.setproperties(obj::AbstractSystem, patch::NamedTuple)
+    if issubset(fieldnames(patch), fieldnames(obj))
+        args = map(fieldnames(obj)) do fn
+            if fn in fieldnames(patch)
+                :(patch.$fn)
+            else
+                :(getfield(obj, $(Meta.quot(fn))))
+            end
+        end
+        return Expr(:block,
+            Expr(:meta, :inline),
+            Expr(:call,:(constructorof($obj)), args...)
+        )
+    else
+        error("This should never happen. Trying to set $(typeof(obj)) with $patch.")
+    end
+end
+
 function Base.getproperty(sys::AbstractSystem, name::Symbol)
     sysname = nameof(sys)
     systems = get_systems(sys)
diff --git a/src/systems/control/controlsystem.jl b/src/systems/control/controlsystem.jl
@@ -1,22 +1,22 @@
 abstract type AbstractControlSystem <: AbstractSystem end
 
 function namespace_controls(sys::AbstractSystem)
-    [rename(x,renamespace(sys.name,x.name)) for x in controls(sys)]
+    [rename(x,renamespace(nameof(sys),nameof(x))) for x in controls(sys)]
 end
 
 function controls(sys::AbstractControlSystem,args...)
     name = last(args)
-    extra_names = reduce(Symbol,[Symbol(:₊,x.name) for x in args[1:end-1]])
+    extra_names = reduce(Symbol,[Symbol(:₊,nameof(x)) for x in args[1:end-1]])
     newname = renamespace(extra_names,name)
-    rename(x,renamespace(sys.name,newname))(sys.iv())
+    rename(x,renamespace(nameof(sys),newname))(get_iv(sys))
 end
 
 function controls(sys::AbstractControlSystem,name::Symbol)
-    x = sys.controls[findfirst(x->x.name==name,sys.ps)]
-    rename(x,renamespace(sys.name,x.name))()
+    x = get_controls(sys)[findfirst(x->nameof(x)==name,sys.ps)]
+    rename(x,renamespace(nameof(sys),nameof(x)))
 end
 
-controls(sys::AbstractControlSystem) = isempty(sys.systems) ? sys.controls : [sys.controls;reduce(vcat,namespace_controls.(sys.systems))]
+controls(sys::AbstractControlSystem) = isempty(get_systems(sys)) ? get_controls(sys) : [get_controls(sys);reduce(vcat,namespace_controls.(get_systems(sys)))]
 
 """
 $(TYPEDEF)
@@ -107,13 +107,13 @@ function (f::ControlToExpr)(O)
     !istree(O) && return O
     res = if isa(operation(O), Sym)
         # normal variables and control variables
-        (any(isequal(O), f.states) || any(isequal(O), f.controls)) && return tosymbol(O)
+        (any(isequal(O), f.states) || any(isequal(O), controls(f))) && return tosymbol(O)
         build_expr(:call, Any[operation(O).name; f.(arguments(O))])
     else
         build_expr(:call, Any[Symbol(operation(O)); f.(arguments(O))])
     end
 end
-(f::ControlToExpr)(x::Sym) = x.name
+(f::ControlToExpr)(x::Sym) = nameof(x)
 
 function constructRadauIIA5(T::Type = Float64)
   sq6 = sqrt(convert(T, 6))
@@ -146,19 +146,24 @@ function runge_kutta_discretize(sys::ControlSystem,dt,tspan;
     n = length(tspan[1]:dt:tspan[2]) - 1
     m = length(tab.α)
 
-    f = @RuntimeGeneratedFunction(build_function([x.rhs for x in equations(sys)],sys.states,sys.controls,sys.ps,sys.iv,conv = ModelingToolkit.ControlToExpr(sys))[1])
-    L = @RuntimeGeneratedFunction(build_function(sys.loss,sys.states,sys.controls,sys.ps,sys.iv,conv = ModelingToolkit.ControlToExpr(sys)))
+    sts = states(sys)
+    ctr = controls(sys)
+    ps = parameters(sys)
+    lo = get_loss(sys)
+    iv = get_iv(sys)
+    f = @RuntimeGeneratedFunction(build_function([x.rhs for x in equations(sys)],sts,ctr,ps,iv,conv = ModelingToolkit.ControlToExpr(sys))[1])
+    L = @RuntimeGeneratedFunction(build_function(lo,sts,ctr,ps,iv,conv = ModelingToolkit.ControlToExpr(sys)))
 
     var(n, i...) = var(nameof(n), i...)
-    var(n::Symbol, i...) = Sym{FnType{Tuple{symtype(sys.iv)}, Real}}(nameof(Variable(n, i...)))
+    var(n::Symbol, i...) = Sym{FnType{Tuple{symtype(iv)}, Real}}(nameof(Variable(n, i...)))
     # Expand out all of the variables in time and by stages
-    timed_vars = [[var(operation(x),i)(sys.iv) for i in 1:n+1] for x in states(sys)]
-    k_vars = [[var(Symbol(:ᵏ,nameof(operation(x))),i,j)(sys.iv) for i in 1:m, j in 1:n] for x in states(sys)]
+    timed_vars = [[var(operation(x),i)(iv) for i in 1:n+1] for x in states(sys)]
+    k_vars = [[var(Symbol(:ᵏ,nameof(operation(x))),i,j)(iv) for i in 1:m, j in 1:n] for x in states(sys)]
     states_timeseries = [getindex.(timed_vars,j) for j in 1:n+1]
     k_timeseries = [[Num.(getindex.(k_vars,i,j)) for i in 1:m] for j in 1:n]
-    control_timeseries = [[[var(operation(x),i,j)(sys.iv) for x in controls(sys)] for i in 1:m] for j in 1:n]
+    control_timeseries = [[[var(operation(x),i,j)(iv) for x in sts] for i in 1:m] for j in 1:n]
     ps = parameters(sys)
-    iv = sys.iv
+    iv = iv
 
     # Calculate all of the update and stage equations
     mult = [tab.A * k_timeseries[i] for i in 1:n]
diff --git a/src/systems/diffeqs/abstractodesystem.jl b/src/systems/diffeqs/abstractodesystem.jl
@@ -244,8 +244,13 @@ function process_DEProblem(constructor, sys::AbstractODESystem,u0map,parammap;
                            kwargs...)
     dvs = states(sys)
     ps = parameters(sys)
-    u0map′ = lower_mapnames(u0map,get_iv(sys))
-    u0 = varmap_to_vars(u0map′,dvs; defaults=default_u0(sys))
+
+    if u0map !== nothing
+        u0map′ = lower_mapnames(u0map,get_iv(sys))
+        u0 = varmap_to_vars(u0map′,dvs; defaults=default_u0(sys))
+    else
+        u0 = nothing
+    end
 
     if !(parammap isa DiffEqBase.NullParameters)
         parammap′ = lower_mapnames(parammap)
diff --git a/src/systems/diffeqs/first_order_transform.jl b/src/systems/diffeqs/first_order_transform.jl
@@ -5,8 +5,9 @@ Takes a Nth order ODESystem and returns a new ODESystem written in first order
 form by defining new variables which represent the N-1 derivatives.
 """
 function ode_order_lowering(sys::ODESystem)
-    eqs_lowered, new_vars = ode_order_lowering(equations(sys), sys.iv, states(sys))
-    return ODESystem(eqs_lowered, sys.iv, new_vars, sys.ps)
+    iv = independent_variable(sys)
+    eqs_lowered, new_vars = ode_order_lowering(equations(sys), iv, states(sys))
+    return ODESystem(eqs_lowered, iv, new_vars, parameters(sys))
 end
 
 function ode_order_lowering(eqs, iv, states)
diff --git a/src/systems/diffeqs/odesystem.jl b/src/systems/diffeqs/odesystem.jl
@@ -182,9 +182,15 @@ function collect_var!(states, parameters, var, iv)
 end
 
 # NOTE: equality does not check cached Jacobian
-Base.:(==)(sys1::ODESystem, sys2::ODESystem) =
-    _eq_unordered(sys1.eqs, sys2.eqs) && isequal(sys1.iv, sys2.iv) &&
-    _eq_unordered(sys1.states, sys2.states) && _eq_unordered(sys1.ps, sys2.ps)
+function Base.:(==)(sys1::ODESystem, sys2::ODESystem)
+    iv1 = independent_variable(sys1)
+    iv2 = independent_variable(sys2)
+    isequal(iv1, iv2) &&
+    _eq_unordered(get_eqs(sys1), get_eqs(sys2)) &&
+    _eq_unordered(get_states(sys1), get_states(sys2)) &&
+    _eq_unordered(get_ps(sys1), get_ps(sys2)) &&
+    all(s1 == s2 for (s1, s2) in zip(get_systems(sys1), get_systems(sys2)))
+end
 
 function flatten(sys::ODESystem)
     systems = get_systems(sys)
diff --git a/src/systems/diffeqs/sdesystem.jl b/src/systems/diffeqs/sdesystem.jl
@@ -101,10 +101,10 @@ function SDESystem(deqs::AbstractVector{<:Equation}, neqs, iv, dvs, ps;
 end
 
 function generate_diffusion_function(sys::SDESystem, dvs = states(sys), ps = parameters(sys); kwargs...)
-    return build_function(sys.noiseeqs,
+    return build_function(get_noiseeqs(sys),
                           map(x->time_varying_as_func(value(x), sys), dvs),
                           map(x->time_varying_as_func(value(x), sys), ps),
-                          sys.iv; kwargs...)
+                          independent_variable(sys); kwargs...)
 end
 
 """
@@ -114,36 +114,34 @@ Choose correction_factor=-1//2 (1//2) to converte Ito -> Stratonovich (Stratonov
 """
 function stochastic_integral_transform(sys::SDESystem, correction_factor)
     # use the general interface
-    if typeof(sys.noiseeqs) <: Vector
-        eqs = vcat([sys.eqs[i].lhs ~ sys.noiseeqs[i] for i in eachindex(sys.states)]...)
-        de = ODESystem(eqs,sys.iv,sys.states,sys.ps)
+    if typeof(get_noiseeqs(sys)) <: Vector
+        eqs = vcat([equations(sys)[i].lhs ~ get_noiseeqs(sys)[i] for i in eachindex(states(sys))]...)
+        de = ODESystem(eqs,get_iv(sys),states(sys),parameters(sys))
 
         jac = calculate_jacobian(de, sparse=false, simplify=false)
-        ∇σσ′ = simplify.(jac*sys.noiseeqs)
+        ∇σσ′ = simplify.(jac*get_noiseeqs(sys))
 
-        deqs = vcat([sys.eqs[i].lhs ~ sys.eqs[i].rhs+ correction_factor*∇σσ′[i] for i in eachindex(sys.states)]...)
+        deqs = vcat([equations(sys)[i].lhs ~ equations(sys)[i].rhs+ correction_factor*∇σσ′[i] for i in eachindex(states(sys))]...)
     else
-        dimstate, m = size(sys.noiseeqs)
-        eqs = vcat([sys.eqs[i].lhs ~ sys.noiseeqs[i] for i in eachindex(sys.states)]...)
-        de = ODESystem(eqs,sys.iv,sys.states,sys.ps)
+        dimstate, m = size(get_noiseeqs(sys))
+        eqs = vcat([equations(sys)[i].lhs ~ get_noiseeqs(sys)[i] for i in eachindex(states(sys))]...)
+        de = ODESystem(eqs,get_iv(sys),states(sys),parameters(sys))
 
         jac = calculate_jacobian(de, sparse=false, simplify=false)
-        ∇σσ′ = simplify.(jac*sys.noiseeqs[:,1])
+        ∇σσ′ = simplify.(jac*get_noiseeqs(sys)[:,1])
         for k = 2:m
-            eqs = vcat([sys.eqs[i].lhs ~ sys.noiseeqs[Int(i+(k-1)*dimstate)] for i in eachindex(sys.states)]...)
-            de = ODESystem(eqs,sys.iv,sys.states,sys.ps)
+            eqs = vcat([equations(sys)[i].lhs ~ get_noiseeqs(sys)[Int(i+(k-1)*dimstate)] for i in eachindex(states(sys))]...)
+            de = ODESystem(eqs,get_iv(sys),states(sys),parameters(sys))
 
             jac = calculate_jacobian(de, sparse=false, simplify=false)
-            ∇σσ′ = ∇σσ′ + simplify.(jac*sys.noiseeqs[:,k])
+            ∇σσ′ = ∇σσ′ + simplify.(jac*get_noiseeqs(sys)[:,k])
         end
 
-        deqs = vcat([sys.eqs[i].lhs ~ sys.eqs[i].rhs + correction_factor*∇σσ′[i] for i in eachindex(sys.states)]...)
+        deqs = vcat([equations(sys)[i].lhs ~ equations(sys)[i].rhs + correction_factor*∇σσ′[i] for i in eachindex(states(sys))]...)
     end
 
 
-    de = SDESystem(deqs,sys.noiseeqs,sys.iv,sys.states,sys.ps)
-
-    de
+    SDESystem(deqs,get_noiseeqs(sys),get_iv(sys),states(sys),parameters(sys))
 end
 
 
@@ -161,7 +159,7 @@ Create an `SDEFunction` from the [`SDESystem`](@ref). The arguments `dvs` and `p
 are used to set the order of the dependent variable and parameter vectors,
 respectively.
 """
-function DiffEqBase.SDEFunction{iip}(sys::SDESystem, dvs = sys.states, ps = sys.ps,
+function DiffEqBase.SDEFunction{iip}(sys::SDESystem, dvs = states(sys), ps = parameters(sys),
                                      u0 = nothing;
                                      version = nothing, tgrad=false, sparse = false,
                                      jac = false, Wfact = false, eval_expression = true, kwargs...) where {iip}
@@ -215,7 +213,7 @@ function DiffEqBase.SDEFunction{iip}(sys::SDESystem, dvs = sys.states, ps = sys.
                      Wfact = _Wfact === nothing ? nothing : _Wfact,
                      Wfact_t = _Wfact_t === nothing ? nothing : _Wfact_t,
                      mass_matrix = _M,
-                     syms = Symbol.(sys.states))
+                     syms = Symbol.(states(sys)))
 end
 
 function DiffEqBase.SDEFunction(sys::SDESystem, args...; kwargs...)
@@ -287,7 +285,7 @@ function SDEFunctionExpr{iip}(sys::SDESystem, dvs = states(sys),
                          Wfact = Wfact,
                          Wfact_t = Wfact_t,
                          mass_matrix = M,
-                         syms = $(Symbol.(states(sys))),kwargs...)
+                         syms = $(Symbol.(states(sys))))
     end
     !linenumbers ? striplines(ex) : ex
 end
@@ -322,13 +320,14 @@ function DiffEqBase.SDEProblem{iip}(sys::SDESystem,u0map,tspan,parammap=DiffEqBa
     f, u0, p = process_DEProblem(SDEFunction{iip}, sys, u0map, parammap; kwargs...)
     sparsenoise === nothing && (sparsenoise = get(kwargs, :sparse, false))
 
-    if typeof(sys.noiseeqs) <: AbstractVector
+    noiseeqs = get_noiseeqs(sys)
+    if noiseeqs isa AbstractVector
         noise_rate_prototype = nothing
     elseif sparsenoise
-        I,J,V = findnz(SparseArrays.sparse(sys.noiseeqs))
+        I,J,V = findnz(SparseArrays.sparse(noiseeqs))
         noise_rate_prototype = SparseArrays.sparse(I,J,zero(eltype(u0)))
     else
-        noise_rate_prototype = zeros(eltype(u0),size(sys.noiseeqs))
+        noise_rate_prototype = zeros(eltype(u0),size(noiseeqs))
     end
 
     SDEProblem{iip}(f,f.g,u0,tspan,p;noise_rate_prototype=noise_rate_prototype,kwargs...)
@@ -363,13 +362,15 @@ function SDEProblemExpr{iip}(sys::SDESystem,u0map,tspan,
     linenumbers = get(kwargs, :linenumbers, true)
     sparsenoise === nothing && (sparsenoise = get(kwargs, :sparse, false))
 
-    if typeof(sys.noiseeqs) <: AbstractVector
+    noiseeqs = get_noiseeqs(sys)
+    if noiseeqs isa AbstractVector
         noise_rate_prototype = nothing
     elseif sparsenoise
-        I,J,V = findnz(SparseArrays.sparse(sys.noiseeqs))
+        I,J,V = findnz(SparseArrays.sparse(noiseeqs))
         noise_rate_prototype = SparseArrays.sparse(I,J,zero(eltype(u0)))
     else
-        noise_rate_prototype = zeros(eltype(u0),size(sys.noiseeqs))
+        T = u0 === nothing ? Float64 : eltype(u0)
+        noise_rate_prototype = zeros(T,size(get_noiseeqs(sys)))
     end
     ex = quote
         f = $f
diff --git a/src/systems/optimization/optimizationsystem.jl b/src/systems/optimization/optimizationsystem.jl
@@ -98,11 +98,10 @@ function generate_function(sys::OptimizationSystem, vs = states(sys), ps = param
                           conv = AbstractSysToExpr(sys),kwargs...)
 end
 
-equations(sys::OptimizationSystem) = isempty(sys.systems) ? sys.op : sys.op + reduce(+,namespace_expr.(sys.systems))
-namespace_expr(sys::OptimizationSystem) = namespace_expr(sys.op,sys.name,nothing)
+equations(sys::OptimizationSystem) = isempty(get_systems(sys)) ? get_op(sys) : get_op(sys) + reduce(+,namespace_expr.(get_systems(sys)))
+namespace_expr(sys::OptimizationSystem) = namespace_expr(get_op(sys),nameof(sys),nothing)
 
-hessian_sparsity(sys::OptimizationSystem) =
-    hessian_sparsity(sys.op, states(sys))
+hessian_sparsity(sys::OptimizationSystem) = hessian_sparsity(get_op(sys), states(sys))
 
 struct AutoModelingToolkit <: DiffEqBase.AbstractADType end
 
diff --git a/src/systems/reaction/reactionsystem.jl b/src/systems/reaction/reactionsystem.jl
@@ -446,7 +446,7 @@ function Base.convert(::Type{<:SDESystem},rs::ReactionSystem, combinatoric_ratel
               (noise_scaling===nothing) ?
                     get_ps(rs) :
                     union(get_ps(rs),toparam.(noise_scaling)),
-                    name=rs.name,systems=convert.(SDESystem,get_systems(rs)))
+                    name=nameof(rs),systems=convert.(SDESystem,get_systems(rs)))
 end
 
 """
diff --git a/test/build_targets.jl b/test/build_targets.jl
@@ -33,25 +33,25 @@ sys = ODESystem(eqs,t,[x,y],[a])
 @test ModelingToolkit.build_function(eqs,[x,y],[a],t,target = ModelingToolkit.CTarget(),
                                      lhsname=:internal_var___du,
                                      rhsnames=[:internal_var___u,:internal_var___p,:t]) ==
-      ModelingToolkit.build_function(sys.eqs,[x,y],[a],t,target = ModelingToolkit.CTarget(),
-                                           lhsname=:internal_var___du,
-                                           rhsnames=[:internal_var___u,:internal_var___p,:t])
+    ModelingToolkit.build_function(equations(sys),[x,y],[a],t,target = ModelingToolkit.CTarget(),
+                                   lhsname=:internal_var___du,
+                                   rhsnames=[:internal_var___u,:internal_var___p,:t])
 
 @test ModelingToolkit.build_function(eqs,[x,y],[a],t,target = ModelingToolkit.StanTarget()) ==
-      ModelingToolkit.build_function(sys.eqs,[x,y],[a],t,target = ModelingToolkit.StanTarget())
+      ModelingToolkit.build_function(equations(sys),[x,y],[a],t,target = ModelingToolkit.StanTarget())
 
 @test ModelingToolkit.build_function(eqs,[x,y],[a],t,target = ModelingToolkit.MATLABTarget()) ==
-      ModelingToolkit.build_function(sys.eqs,[x,y],[a],t,target = ModelingToolkit.MATLABTarget())
+      ModelingToolkit.build_function(equations(sys),[x,y],[a],t,target = ModelingToolkit.MATLABTarget())
 
 @test ModelingToolkit.build_function(eqs,[x,y],[a],t,target = ModelingToolkit.CTarget(),
                                      lhsname=:internal_var___du,
                                      rhsnames=[:internal_var___u,:internal_var___p,:t]) ==
-      ModelingToolkit.build_function(sys.eqs,sys.states,sys.ps,sys.iv,target = ModelingToolkit.CTarget(),
+      ModelingToolkit.build_function(equations(sys),states(sys),parameters(sys),independent_variable(sys),target = ModelingToolkit.CTarget(),
                                      lhsname=:internal_var___du,
                                      rhsnames=[:internal_var___u,:internal_var___p,:t])
 
 @test ModelingToolkit.build_function(eqs,[x,y],[a],t,target = ModelingToolkit.StanTarget()) ==
-      ModelingToolkit.build_function(sys.eqs,sys.states,sys.ps,sys.iv,target = ModelingToolkit.StanTarget())
+      ModelingToolkit.build_function(equations(sys),states(sys),parameters(sys),independent_variable(sys),target = ModelingToolkit.StanTarget())
 
 @test ModelingToolkit.build_function(eqs,[x,y],[a],t,target = ModelingToolkit.MATLABTarget()) ==
-      ModelingToolkit.build_function(sys.eqs,sys.states,sys.ps,sys.iv,target = ModelingToolkit.MATLABTarget())
+      ModelingToolkit.build_function(equations(sys),states(sys),parameters(sys),independent_variable(sys),target = ModelingToolkit.MATLABTarget())
diff --git a/test/lowering_solving.jl b/test/lowering_solving.jl
diff --git a/test/odesystem.jl b/test/odesystem.jl
diff --git a/test/reactionsystem.jl b/test/reactionsystem.jl
diff --git a/test/reduction.jl b/test/reduction.jl
diff --git a/test/sdesystem.jl b/test/sdesystem.jl
