Instead of collecting constants up front and adding them to the fields of a system, collect the constants just prior to building the functions.

Author: lamorton

src/systems/abstractsystem.jl

@@ -167,7 +167,6 @@ for prop in [:eqs
              :iv
              :states
              :ps
-             :cs
              :var_to_name
              :ctrls
              :defaults
@@ -377,7 +376,6 @@ end
 namespace_variables(sys::AbstractSystem) = states(sys, states(sys))
 namespace_parameters(sys::AbstractSystem) = parameters(sys, parameters(sys))
 namespace_controls(sys::AbstractSystem) = controls(sys, controls(sys))
-namespace_constants(sys::AbstractSystem) = constants(sys, constants(sys))
 
 function namespace_defaults(sys)
     defs = defaults(sys)
@@ -439,12 +437,6 @@ function parameters(sys::AbstractSystem)
     unique(isempty(systems) ? ps : [ps; reduce(vcat, namespace_parameters.(systems))])
 end
 
-function constants(sys::AbstractSystem)
-    cs = get_cs(sys)
-    systems = get_systems(sys)
-    unique(isempty(systems) ? cs : [cs; reduce(vcat, namespace_constants.(systems))])
-end
-
 function controls(sys::AbstractSystem)
     ctrls = get_ctrls(sys)
     systems = get_systems(sys)

src/systems/diffeqs/abstractodesystem.jl

@@ -128,7 +128,7 @@ function generate_function(sys::AbstractODESystem, dvs = states(sys), ps = param
            [eq.rhs for eq in eqs]
 
     # Swap constants for their values
-    cs = constants(sys)
+    cs = collect_constants(eqs)
     if !isempty(cs) > 0
         cmap = map(x -> x => getdefault(x), cs)
         rhss = map(x -> substitute(x, cmap), rhss)
@@ -142,6 +142,7 @@ function generate_function(sys::AbstractODESystem, dvs = states(sys), ps = param
     pre, sol_states = get_substitutions_and_solved_states(sys,
                                                           no_postprocess = has_difference)
 
+
     if implicit_dae
         build_function(rhss, ddvs, u, p, t; postprocess_fbody = pre, states = sol_states,
                        kwargs...)

src/systems/diffeqs/odesystem.jl

@@ -36,8 +36,6 @@ struct ODESystem <: AbstractODESystem
     states::Vector
     """Parameter variables. Must not contain the independent variable."""
     ps::Vector
-    """Symbolic constants."""
-    cs::Vector
     """Array variables."""
     var_to_name::Any
     """Control parameters (some subset of `ps`)."""
@@ -122,7 +120,7 @@ struct ODESystem <: AbstractODESystem
     """
     metadata::Any
 
-    function ODESystem(deqs, iv, dvs, ps, cs, var_to_name, ctrls, observed, tgrad,
+    function ODESystem(deqs, iv, dvs, ps, var_to_name, ctrls, observed, tgrad,
                        jac, ctrl_jac, Wfact, Wfact_t, name, systems, defaults,
                        torn_matching, connector_type, connections, preface, cevents,
                        devents, tearing_state = nothing, substitutions = nothing,
@@ -135,16 +133,16 @@ struct ODESystem <: AbstractODESystem
             check_equations(equations(cevents), iv)
         end
         if checks == true || (checks & CheckUnits) > 0
-            all_dimensionless([dvs; ps; iv; cs]) || check_units(deqs)
+            all_dimensionless([dvs; ps; iv]) || check_units(deqs)
         end
-        new(deqs, iv, dvs, ps, cs, var_to_name, ctrls, observed, tgrad, jac,
+        new(deqs, iv, dvs, ps, var_to_name, ctrls, observed, tgrad, jac,
             ctrl_jac, Wfact, Wfact_t, name, systems, defaults, torn_matching,
             connector_type, connections, preface, cevents, devents, tearing_state,
             substitutions, metadata)
     end
 end
 
-function ODESystem(deqs::AbstractVector{<:Equation}, iv, dvs, ps, cs;
+function ODESystem(deqs::AbstractVector{<:Equation}, iv, dvs, ps;
                    controls = Num[],
                    observed = Equation[],
                    systems = ODESystem[],
@@ -166,7 +164,6 @@ function ODESystem(deqs::AbstractVector{<:Equation}, iv, dvs, ps, cs;
     iv′ = value(scalarize(iv))
     dvs′ = value.(scalarize(dvs))
     ps′ = value.(scalarize(ps))
-    cs′ = value.(scalarize(cs))
     ctrl′ = value.(scalarize(controls))
 
     if !(isempty(default_u0) && isempty(default_p))
@@ -179,7 +176,6 @@ function ODESystem(deqs::AbstractVector{<:Equation}, iv, dvs, ps, cs;
     var_to_name = Dict()
     process_variables!(var_to_name, defaults, dvs′)
     process_variables!(var_to_name, defaults, ps′)
-    process_variables!(var_to_name, defaults, cs′)
     isempty(observed) || collect_var_to_name!(var_to_name, (eq.lhs for eq in observed))
 
     tgrad = RefValue(EMPTY_TGRAD)
@@ -193,7 +189,7 @@ function ODESystem(deqs::AbstractVector{<:Equation}, iv, dvs, ps, cs;
     end
     cont_callbacks = SymbolicContinuousCallbacks(continuous_events)
     disc_callbacks = SymbolicDiscreteCallbacks(discrete_events)
-    ODESystem(deqs, iv′, dvs′, ps′, cs′, var_to_name, ctrl′, observed, tgrad, jac,
+    ODESystem(deqs, iv′, dvs′, ps′, var_to_name, ctrl′, observed, tgrad, jac,
               ctrl_jac, Wfact, Wfact_t, name, systems, defaults, nothing,
               connector_type, nothing, preface, cont_callbacks, disc_callbacks,
               metadata, checks = checks)
@@ -205,7 +201,6 @@ function ODESystem(eqs, iv = nothing; kwargs...)
     diffvars = OrderedSet()
     allstates = OrderedSet()
     ps = OrderedSet()
-    cs = OrderedSet() #Constants
     # reorder equations such that it is in the form of `diffeq, algeeq`
     diffeq = Equation[]
     algeeq = Equation[]
@@ -223,8 +218,8 @@ function ODESystem(eqs, iv = nothing; kwargs...)
     compressed_eqs = Equation[] # equations that need to be expanded later, like `connect(a, b)`
     for eq in eqs
         eq.lhs isa Union{Symbolic, Number} || (push!(compressed_eqs, eq); continue)
-        collect_vars!(allstates, ps, cs, eq.lhs, iv)
-        collect_vars!(allstates, ps, cs, eq.rhs, iv)
+        collect_vars!(allstates, ps, eq.lhs, iv)
+        collect_vars!(allstates, ps, eq.rhs, iv)
         if isdiffeq(eq)
             diffvar, _ = var_from_nested_derivative(eq.lhs)
             isequal(iv, iv_from_nested_derivative(eq.lhs)) ||
@@ -240,7 +235,16 @@ function ODESystem(eqs, iv = nothing; kwargs...)
     algevars = setdiff(allstates, diffvars)
     # the orders here are very important!
     return ODESystem(Equation[diffeq; algeeq; compressed_eqs], iv,
-                     collect(Iterators.flatten((diffvars, algevars))), ps, cs; kwargs...)
+                     collect(Iterators.flatten((diffvars, algevars))), ps; kwargs...)
+end
+
+function collect_constants(eqs) #Does this need to be different for other system types?
+    constants = Set()
+    for eq in eqs
+        collect_constants!(constants,eq.lhs)
+        collect_constants!(constants,eq.rhs)
+    end
+    return collect(constants)
 end
 
 # NOTE: equality does not check cached Jacobian
@@ -265,7 +269,6 @@ function flatten(sys::ODESystem, noeqs = false)
                          get_iv(sys),
                          states(sys),
                          parameters(sys),
-                         constants(sys),
                          observed = observed(sys),
                          continuous_events = continuous_events(sys),
                          discrete_events = discrete_events(sys),
@@ -298,6 +301,13 @@ function build_explicit_observed_function(sys, ts;
     dep_vars = scalarize(setdiff(vars, ivs))
 
     obs = observed(sys)
+
+    cs = collect_constants(obs)
+    if !isempty(cs) > 0
+        cmap = map(x -> x => getdefault(x), cs)
+        obs = map(x -> x.lhs ~ substitute(x.rhs, cmap), obs)
+    end
+
     sts = Set(states(sys))
     observed_idx = Dict(x.lhs => i for (i, x) in enumerate(obs))
     namespaced_to_obs = Dict(states(sys, x.lhs) => x.lhs for x in obs)

src/utils.jl

@@ -466,46 +466,61 @@ function find_derivatives!(vars, expr, f)
     return vars
 end
 
-function collect_vars!(states, parameters, constants, expr, iv)
+function collect_vars!(states, parameters, expr, iv)
     if expr isa Sym
-        collect_var!(states, parameters, constants, expr, iv)
+        collect_var!(states, parameters, expr, iv)
     else
         for var in vars(expr)
             if istree(var) && operation(var) isa Differential
                 var, _ = var_from_nested_derivative(var)
             end
-            collect_var!(states, parameters, constants, var, iv)
+            collect_var!(states, parameters, var, iv)
         end
     end
     return nothing
 end
 
-function collect_vars_difference!(states, parameters, constants, expr, iv)
+function collect_constants!(constants, expr)
     if expr isa Sym
-        collect_var!(states, parameters, constants, expr, iv)
+        collect_constant!(constants, expr)
+    else
+        for var in vars(expr)
+            collect_constant!(constants, var)
+        end
+    end
+    return nothing
+end
+
+function collect_vars_difference!(states, parameters, expr, iv)
+    if expr isa Sym
+        collect_var!(states, parameters, expr, iv)
     else
         for var in vars(expr)
             if istree(var) && operation(var) isa Difference
                 var, _ = var_from_nested_difference(var)
             end
-            collect_var!(states, parameters, constants, var, iv)
+            collect_var!(states, parameters, var, iv)
         end
     end
     return nothing
 end
 
-function collect_var!(states, parameters, constants, var, iv)
+function collect_var!(states, parameters, var, iv)
     isequal(var, iv) && return nothing
     if isparameter(var) || (istree(var) && isparameter(operation(var)))
         push!(parameters, var)
-    elseif isconstant(var)
-        push!(constants,var)
-    else
+    elseif !isconstant(var)
         push!(states, var)
     end
     return nothing
 end
 
+function collect_constant!(constants, var)
+    if isconstant(var)
+        push!(constants,var)
+    end
+end
+
 function get_postprocess_fbody(sys)
     if has_preface(sys) && (pre = preface(sys); pre !== nothing)
         pre_ = let pre = pre
@@ -549,6 +564,13 @@ function get_substitutions_and_solved_states(sys; no_postprocess = false)
         pre = no_postprocess ? (ex -> ex) : get_postprocess_fbody(sys)
     else
         @unpack subs = get_substitutions(sys)
+        # Swap constants for their values
+        cs = collect_constants(subs)
+        if !isempty(cs) > 0
+            cmap = map(x -> x => getdefault(x), cs)
+            subs = map(x -> x.lhs ~ substitute(x.rhs, cmap), subs)
+        end
+
         sol_states = Code.NameState(Dict(eq.lhs => Symbol(eq.lhs) for eq in subs))
         if no_postprocess
             pre = ex -> Let(Assignment[Assignment(eq.lhs, eq.rhs) for eq in subs], ex,

test/constants.jl

@@ -9,13 +9,17 @@ MT = ModelingToolkit
 D = Differential(t)
 eqs = [D(x) ~ a]
 @named sys = ODESystem(eqs)
-prob = ODEProblem(sys, [0, ], [0.0, 1.0],[])
-sol = solve(prob,Tsit5())
+prob = ODEProblem(sys, [0, ], [0.0, 1.0], [])
+sol = solve(prob, Tsit5())
 
-# Test structural_simplify handling
-eqs = [D(x) ~ t,
+# Test structural_simplify substitutions & observed values
+eqs = [D(x) ~ 1,
     w ~ a]
 @named sys = ODESystem(eqs)
 simp = structural_simplify(sys);
-@test isequal(simp.substitutions.subs[1], w~a)
+@test isequal(simp.substitutions.subs[1], eqs[2])
+@test isequal(equations(simp)[1], eqs[1])
+prob = ODEProblem(simp, [0, ], [0.0, 1.0], [])
+sol = solve(prob, Tsit5())
+@test sol[w][1] == 1 
 

test/runtests.jl

@@ -48,4 +48,4 @@ println("Last test requires gcc available in the path!")
 @safetestset "FuncAffect Test" begin include("funcaffect.jl") end
 
 # Reference tests go Last
-#@safetestset "Latexify recipes Test" begin include("latexify.jl") end
+@safetestset "Latexify recipes Test" begin include("latexify.jl") end