refactor: rename states to unknowns, and related changes

- also some bug fixes

docs/src/basics/AbstractSystem.md

@@ -38,18 +38,18 @@ same keyword arguments, which are:
 
 Each `AbstractSystem` has lists of variables in context, such as distinguishing
 parameters vs states. In addition, an `AbstractSystem` can also hold other
-`AbstractSystem` types. Direct accessing of the values, such as `sys.states`,
-gives the immediate list, while the accessor functions `states(sys)` gives the
+`AbstractSystem` types. Direct accessing of the values, such as `sys.unknowns`,
+gives the immediate list, while the accessor functions `unknowns(sys)` gives the
 total set, which includes that of all systems held inside.
 
 The values which are common to all `AbstractSystem`s are:
 
   - `equations(sys)`: All equations that define the system and its subsystems.
-  - `states(sys)`: All the states in the system and its subsystems.
+  - `unknowns(sys)`: All the unknowns in the system and its subsystems.
   - `parameters(sys)`: All parameters of the system and its subsystems.
   - `nameof(sys)`: The name of the current-level system.
   - `get_eqs(sys)`: Equations that define the current-level system.
-  - `get_states(sys)`: States that are in the current-level system.
+  - `get_unknowns(sys)`: States that are in the current-level system.
   - `get_ps(sys)`: Parameters that are in the current-level system.
   - `get_systems(sys)`: Subsystems of the current-level system.
 

docs/src/basics/Linearization.md

@@ -38,7 +38,7 @@ The named tuple `matrices` contains the matrices of the linear statespace repres
 
 ```@example LINEARIZE
 using ModelingToolkit: inputs, outputs
-[states(simplified_sys); inputs(simplified_sys); outputs(simplified_sys)]
+[unknowns(simplified_sys); inputs(simplified_sys); outputs(simplified_sys)]
 ```
 
 ## Operating point

docs/src/examples/modelingtoolkitize_index_reduction.md

@@ -32,7 +32,7 @@ traced_sys = modelingtoolkitize(pendulum_prob)
 pendulum_sys = structural_simplify(dae_index_lowering(traced_sys))
 prob = ODAEProblem(pendulum_sys, [], tspan)
 sol = solve(prob, Tsit5(), abstol = 1e-8, reltol = 1e-8)
-plot(sol, idxs = states(traced_sys))
+plot(sol, idxs = unknowns(traced_sys))
 ```
 
 ## Explanation
@@ -154,10 +154,10 @@ prob = ODEProblem(pendulum_sys, Pair[], tspan)
 sol = solve(prob, Rodas4())
 
 using Plots
-plot(sol, idxs = states(traced_sys))
+plot(sol, idxs = unknowns(traced_sys))
 ```
 
-Note that plotting using `states(traced_sys)` is done so that any
+Note that plotting using `unknowns(traced_sys)` is done so that any
 variables which are symbolically eliminated, or any variable reordering
 done for enhanced parallelism/performance, still show up in the resulting
 plot and the plot is shown in the same order as the original numerical
@@ -173,7 +173,7 @@ traced_sys = modelingtoolkitize(pendulum_prob)
 pendulum_sys = structural_simplify(dae_index_lowering(traced_sys))
 prob = ODAEProblem(pendulum_sys, Pair[], tspan)
 sol = solve(prob, Tsit5(), abstol = 1e-8, reltol = 1e-8)
-plot(sol, idxs = states(traced_sys))
+plot(sol, idxs = unknowns(traced_sys))
 ```
 
 And there you go: this has transformed the model from being too hard to

docs/src/systems/JumpSystem.md

@@ -9,7 +9,7 @@ JumpSystem
 ## Composition and Accessor Functions
 
   - `get_eqs(sys)` or `equations(sys)`: The equations that define the jump system.
-  - `get_states(sys)` or `states(sys)`: The set of states in the jump system.
+  - `get_unknowns(sys)` or `unknowns(sys)`: The set of unknowns in the jump system.
   - `get_ps(sys)` or `parameters(sys)`: The parameters of the jump system.
   - `get_iv(sys)`: The independent variable of the jump system.
 

docs/src/systems/NonlinearSystem.md

@@ -9,7 +9,7 @@ NonlinearSystem
 ## Composition and Accessor Functions
 
   - `get_eqs(sys)` or `equations(sys)`: The equations that define the nonlinear system.
-  - `get_states(sys)` or `states(sys)`: The set of states in the nonlinear system.
+  - `get_unknowns(sys)` or `unknowns(sys)`: The set of unknowns in the nonlinear system.
   - `get_ps(sys)` or `parameters(sys)`: The parameters of the nonlinear system.
   - `get_u0_p(sys, u0map, parammap)` Numeric arrays for the initial condition and parameters given `var => value` maps.
 

docs/src/systems/ODESystem.md

@@ -9,7 +9,7 @@ ODESystem
 ## Composition and Accessor Functions
 
   - `get_eqs(sys)` or `equations(sys)`: The equations that define the ODE.
-  - `get_states(sys)` or `states(sys)`: The set of states in the ODE.
+  - `get_unknowns(sys)` or `unknowns(sys)`: The set of unknowns in the ODE.
   - `get_ps(sys)` or `parameters(sys)`: The parameters of the ODE.
   - `get_iv(sys)`: The independent variable of the ODE.
   - `get_u0_p(sys, u0map, parammap)` Numeric arrays for the initial condition and parameters given `var => value` maps.

docs/src/systems/OptimizationSystem.md

@@ -9,7 +9,7 @@ OptimizationSystem
 ## Composition and Accessor Functions
 
   - `get_op(sys)`: The objective to be minimized.
-  - `get_states(sys)` or `states(sys)`: The set of states for the optimization.
+  - `get_unknowns(sys)` or `unknowns(sys)`: The set of unknowns for the optimization.
   - `get_ps(sys)` or `parameters(sys)`: The parameters for the optimization.
   - `get_constraints(sys)` or `constraints(sys)`: The constraints for the optimization.
 

docs/src/systems/SDESystem.md

@@ -16,7 +16,7 @@ sde = SDESystem(ode, noiseeqs)
 ## Composition and Accessor Functions
 
   - `get_eqs(sys)` or `equations(sys)`: The equations that define the SDE.
-  - `get_states(sys)` or `states(sys)`: The set of states in the SDE.
+  - `get_unknowns(sys)` or `unknowns(sys)`: The set of unknowns in the SDE.
   - `get_ps(sys)` or `parameters(sys)`: The parameters of the SDE.
   - `get_iv(sys)`: The independent variable of the SDE.
 

ext/MTKBifurcationKitExt.jl

@@ -30,12 +30,12 @@ struct ObservableRecordFromSolution{S, T}
             p_vals) where {S, T}
         obs_eqs = observed(nsys)
         target_obs_idx = findfirst(isequal(plot_var, eq.lhs) for eq in observed(nsys))
-        state_end_idxs = length(states(nsys))
+        state_end_idxs = length(unknowns(nsys))
         param_end_idxs = state_end_idxs + length(parameters(nsys))
 
         bif_par_idx = state_end_idxs + bif_idx
         # Gets the (base) substitution values for states.
-        subs_vals_states = Pair.(states(nsys), u0_vals)
+        subs_vals_states = Pair.(unknowns(nsys), u0_vals)
         # Gets the (base) substitution values for parameters.
         subs_vals_params = Pair.(parameters(nsys), p_vals)
         # Gets the (base) substitution values for observables.
@@ -95,16 +95,16 @@ function BifurcationKit.BifurcationProblem(nsys::NonlinearSystem,
 
     # Converts the input state guess.
     u0_bif_vals = ModelingToolkit.varmap_to_vars(u0_bif,
-        states(nsys);
+        unknowns(nsys);
         defaults = nsys.defaults)
     p_vals = ModelingToolkit.varmap_to_vars(ps, parameters(nsys); defaults = nsys.defaults)
 
     # Computes bifurcation parameter and the plotting function.
     bif_idx = findfirst(isequal(bif_par), parameters(nsys))
     if !isnothing(plot_var)
         # If the plot var is a normal state.
-        if any(isequal(plot_var, var) for var in states(nsys))
-            plot_idx = findfirst(isequal(plot_var), states(nsys))
+        if any(isequal(plot_var, var) for var in unknowns(nsys))
+            plot_idx = findfirst(isequal(plot_var), unknowns(nsys))
             record_from_solution = (x, p) -> x[plot_idx]
 
             # If the plot var is an observed state.
@@ -134,7 +134,7 @@ end
 # When input is a ODESystem.
 function BifurcationKit.BifurcationProblem(osys::ODESystem, args...; kwargs...)
     nsys = NonlinearSystem([0 ~ eq.rhs for eq in equations(osys)],
-        states(osys),
+        unknowns(osys),
         parameters(osys);
         name = nameof(osys))
     return BifurcationKit.BifurcationProblem(nsys, args...; kwargs...)

src/ModelingToolkit.jl

@@ -36,7 +36,7 @@ using PrecompileTools, Reexport
     using RecursiveArrayTools
 
     using SymbolicIndexingInterface
-    export independent_variables, states, parameters
+    export independent_variables, unknowns, parameters
     import SymbolicUtils
     import SymbolicUtils: istree, arguments, operation, similarterm, promote_symtype,
         Symbolic, isadd, ismul, ispow, issym, FnType,
@@ -196,7 +196,6 @@ export ODEProblem, SDEProblem
 export NonlinearFunction, NonlinearFunctionExpr
 export NonlinearProblem, BlockNonlinearProblem, NonlinearProblemExpr
 export OptimizationProblem, OptimizationProblemExpr, constraints
-export AutoModelingToolkit
 export SteadyStateProblem, SteadyStateProblemExpr
 export JumpProblem, DiscreteProblem
 export NonlinearSystem, OptimizationSystem, ConstraintsSystem

src/inputoutput.jl

@@ -5,7 +5,7 @@ using Symbolics: get_variables
 Return all variables that mare marked as inputs. See also [`unbound_inputs`](@ref)
 See also [`bound_inputs`](@ref), [`unbound_inputs`](@ref)
 """
-inputs(sys) = [filter(isinput, states(sys)); filter(isinput, parameters(sys))]
+inputs(sys) = [filter(isinput, unknowns(sys)); filter(isinput, parameters(sys))]
 
 """
     outputs(sys)
@@ -17,7 +17,7 @@ function outputs(sys)
     o = observed(sys)
     rhss = [eq.rhs for eq in o]
     lhss = [eq.lhs for eq in o]
-    unique([filter(isoutput, states(sys))
+    unique([filter(isoutput, unknowns(sys))
         filter(isoutput, parameters(sys))
         filter(x -> istree(x) && isoutput(x), rhss) # observed can return equations with complicated expressions, we are only looking for single Terms
         filter(x -> istree(x) && isoutput(x), lhss)])
@@ -205,7 +205,7 @@ function generate_control_function(sys::AbstractODESystem, inputs = unbound_inpu
 
     sys, _ = io_preprocessing(sys, inputs, []; simplify, kwargs...)
 
-    dvs = states(sys)
+    dvs = unknowns(sys)
     ps = parameters(sys)
     ps = setdiff(ps, inputs)
     if disturbance_inputs !== nothing
@@ -300,7 +300,7 @@ function inputs_to_parameters!(state::TransformationState, io)
 
     @set! sys.eqs = isempty(input_to_parameters) ? equations(sys) :
                     fast_substitute(equations(sys), input_to_parameters)
-    @set! sys.states = setdiff(states(sys), keys(input_to_parameters))
+    @set! sys.unknowns = setdiff(unknowns(sys), keys(input_to_parameters))
     ps = parameters(sys)
 
     if io !== nothing

src/structural_transformation/StructuralTransformations.jl

@@ -11,7 +11,7 @@ using SymbolicUtils: similarterm, istree
 
 using ModelingToolkit
 using ModelingToolkit: ODESystem, AbstractSystem, var_from_nested_derivative, Differential,
-    states, equations, vars, Symbolic, diff2term, value,
+    unknowns, equations, vars, Symbolic, diff2term, value,
     operation, arguments, Sym, Term, simplify, solve_for,
     isdiffeq, isdifferential, isirreducible,
     empty_substitutions, get_substitutions,

src/structural_transformation/codegen.jl

@@ -298,7 +298,7 @@ function build_torn_function(sys;
         rhss)
 
     states = Any[fullvars[i] for i in states_idxs]
-    @set! sys.unknown_states = states
+    @set! sys.solved_unknowns = states
     syms = map(Symbol, states)
 
     pre = get_postprocess_fbody(sys)
@@ -410,9 +410,9 @@ function build_observed_function(state, ts, var_eq_matching, var_sccs,
     required_algvars = Set(intersect(algvars, vars))
     obs = observed(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)
-    namespaced_to_sts = Dict(states(sys, x) => x for x in states(sys))
-    sts = Set(states(sys))
+    namespaced_to_obs = Dict(unknowns(sys, x.lhs) => x.lhs for x in obs)
+    namespaced_to_sts = Dict(unknowns(sys, x) => x for x in unknowns(sys))
+    sts = Set(unknowns(sys))
 
     # FIXME: This is a rather rough estimate of dependencies. We assume
     # the expression depends on everything before the `maxidx`.

src/structural_transformation/pantelides.jl

@@ -65,7 +65,7 @@ function pantelides_reassemble(state::TearingState, var_eq_matching)
             out_eqs[sort(filter(x -> x !== unassigned, var_eq_matching))]))
 
     @set! sys.eqs = final_eqs
-    @set! sys.states = final_vars
+    @set! sys.unknowns = final_vars
     return sys
 end
 

src/structural_transformation/symbolics_tearing.jl

@@ -540,7 +540,7 @@ function tearing_reassemble(state::TearingState, var_eq_matching;
 
     sys = state.sys
     @set! sys.eqs = neweqs
-    @set! sys.states = Any[v
+    @set! sys.unknowns = Any[v
                            for (i, v) in enumerate(fullvars)
                                if diff_to_var[i] === nothing && ispresent(i)]
     @set! sys.substitutions = Substitutions(subeqs, deps)