Merge pull request #785 from SciML/states_to_unknowns

Change states to unknowns

Author: TorkelE

Project.toml

@@ -44,7 +44,7 @@ HomotopyContinuation = "2.9"
 LaTeXStrings = "1.3.0"
 Latexify = "0.14, 0.15, 0.16"
 MacroTools = "0.5.5"
-ModelingToolkit = "8.73"
+ModelingToolkit = "9.5"
 Parameters = "0.12"
 Reexport = "0.2, 1.0"
 Requires = "1.0"

docs/Project.toml

@@ -39,7 +39,7 @@ Distributions = "0.25"
 Documenter = "0.27"
 HomotopyContinuation = "2.6"
 Latexify = "0.15, 0.16"
-ModelingToolkit = "8.47"
+ModelingToolkit = "9.5"
 NonlinearSolve = "3.4.0"
 Optim = "1"
 Optimization = "3.19"

docs/src/api/catalyst_api.md

@@ -95,11 +95,11 @@ retrieve info from just a base [`ReactionSystem`](@ref) `rn`, ignoring
 sub-systems of `rn`, one can use the ModelingToolkit accessors (these provide
 direct access to the corresponding internal fields of the `ReactionSystem`)
 
-* `ModelingToolkit.get_states(rn)` is a vector that collects all the species
+* `ModelingToolkit.get_unknowns(rn)` is a vector that collects all the species
   defined within `rn`, ordered by species and then non-species variables.
 * `Catalyst.get_species(rn)` is a vector of all the species variables in the system. The
   entries in `get_species(rn)` correspond to the first `length(get_species(rn))`
-  components in `get_states(rn)`.
+  components in `get_unknowns(rn)`.
 * `ModelingToolkit.get_ps(rn)` is a vector that collects all the parameters
   defined *within* reactions in `rn`.
 * `ModelingToolkit.get_eqs(rn)` is a vector that collects all the
@@ -120,10 +120,10 @@ To retrieve information from the full reaction network represented by a system
 `rn`, which corresponds to information within both `rn` and all sub-systems, one
 can call:
 
-* `ModelingToolkit.states(rn)` returns all species *and variables* across the
+* `ModelingToolkit.unknowns(rn)` returns all species *and variables* across the
   system, *all sub-systems*, and all constraint systems. Species are ordered
-  before non-species variables in `states(rn)`, with the first `numspecies(rn)`
-  entires in `states(rn)` being the same as `species(rn)`.
+  before non-species variables in `unknowns(rn)`, with the first `numspecies(rn)`
+  entires in `unknowns(rn)` being the same as `species(rn)`.
 * [`species(rn)`](@ref) is a vector collecting all the chemical species within
   the system and any sub-systems that are also `ReactionSystems`.
 * `ModelingToolkit.parameters(rn)` returns all parameters across the
@@ -196,7 +196,6 @@ ModelingToolkit.extend
 ModelingToolkit.compose
 Catalyst.flatten
 merge!(network1::ReactionSystem, network2::ReactionSystem)
-reorder_states!
 ```
 
 ## Network analysis and representations

docs/src/catalyst_applications/advanced_simulations.md

@@ -173,7 +173,7 @@ want the callback to trigger. The `affect!` function determines what happens to
 the simulation when the callback is triggered. It takes a single object, an
 `integrator` and makes some modification to it (please read more about
 integrators [here](https://docs.sciml.ai/DiffEqDocs/stable/basics/integrator/)).
-Here, we access the system's state `X` through the `integrator`, and add
+Here, we access the system's unknown `X` through the `integrator`, and add
 `5.0` to its amount. We can now simulate our system using the
 callback:
 ```@example ex2

docs/src/catalyst_applications/bifurcation_diagrams.md

@@ -86,7 +86,7 @@ Here, in the bistable region, we only see a single branch. The reason is that th
 
 
 ## Systems with conservation laws
-Some systems are under-determined at steady-state, so that for a given parameter set they have an infinite number of possible steady state solutions, preventing bifurcation diagrams from being computed. Similar to when we [compute steady states for fixed parameter values](@ref homotopy_continuation_conservation_laws), we can utilise Catalyst's ability to detect and eliminate conservation laws to resolve this issue. This requires us to provide information of the species concentrations at which we wish to compute the bifurcation diagram (to determine the values of conserved quantities). These are provided to the `BifurcationProblem` using the `u0` argument.
+Some systems are under-determined at steady state, so that for a given parameter set they have an infinite number of possible steady state solutions, preventing bifurcation diagrams from being computed. Similar to when we [compute steady states for fixed parameter values](@ref homotopy_continuation_conservation_laws), we can utilise Catalyst's ability to detect and eliminate conservation laws to resolve this issue. This requires us to provide information of the species concentrations at which we wish to compute the bifurcation diagram (to determine the values of conserved quantities). These are provided to the `BifurcationProblem` using the `u0` argument.
 
 To illustrate this, we will create a simple model of a kinase that is produced and degraded (at rates $p$ and $d$). The kinase facilitates the phosphorylation of a protein ($X$), which is dephosphorylated at a constant rate. For this system, we will compute a bifurcation diagram, showing how the concentration of the phosphorylated protein ($Xp$) depends on the degradation rate of the kinase ($d$). We will set the total amount of protein ($X+Xp$) to $1.0$.
 ```@example ex2

docs/src/catalyst_applications/simulation_structure_interfacing.md

@@ -1,7 +1,7 @@
 # [Interfacing problems, integrators, and solutions](@id simulation_structure_interfacing)
 When simulating a model, one begins with creating a [problem](https://docs.sciml.ai/DiffEqDocs/stable/basics/problem/). Next, a simulation is performed on a problem, during which the state of the simulation is recorded through an [integrator](https://docs.sciml.ai/DiffEqDocs/stable/basics/integrator/). Finally, the simulation output is returned as a [solution](https://docs.sciml.ai/DiffEqDocs/stable/basics/solution/). This tutorial describes how to access, or modify the state, or parameter, values of problems, integrators, and solutions structures.
 
-Generally, when we have a structure `simulation_struct` and want to interface with the state (or parameter) `G`, we use `simulation_struct[:G]` to access the value, and `simulation_struct[:G] = 5.0` to set it to a new value. However, see the following examples for full details.
+Generally, when we have a structure `simulation_struct` and want to interface with the unknown (or parameter) `G`, we use `simulation_struct[:G]` to access the value, and `simulation_struct[:G] = 5.0` to set it to a new value. However, see the following examples for full details.
 
 ## Interfacing problem objects
 
@@ -27,7 +27,7 @@ with the notation being identical for parameters:
 oprob[:k1]
 ```
 
-If we want to change a state's initial condition value, we use the following notation
+If we want to change an unknown's initial condition value, we use the following notation
 ```@example ex1
 oprob[:X1] = 10.0
 ```
@@ -68,7 +68,7 @@ During a simulation, the solution is stored in an integrator object, we will her
 ```@example ex1
 integrator = init(oprob)
 ```
-Using a similar syntax to problems, we can get the current values of a state within the integrator:
+Using a similar syntax to problems, we can get the current values of an unknown within the integrator:
 ```@example ex1
 integrator[:X1]
 ```
@@ -89,15 +89,15 @@ Finally, we consider solution objects. First, we simulate our problem:
 ```@example ex1
 sol = solve(oprob)
 ```
-For solutions, when we access a state, we get its whole simulation vector:
+For solutions, when we access an unknown, we get its whole simulation vector:
 ```@example ex1
 sol[:X1]
 ```
 while when we access a parameter we only get a single value:
 ```@example ex1
 sol[:k1]
 ```
-Finally, we note that we cannot change the values of solution states or parameters (i.e. both `sol[:X1] = 0.0` and `sol[:k1] = 0.0` generate errors).
+Finally, we note that we cannot change the values of solution unknowns or parameters (i.e. both `sol[:X1] = 0.0` and `sol[:k1] = 0.0` generate errors).
 
 ## [Interfacing using symbolic representation](@id simulation_structure_interfacing_symbolic_representation)
 
@@ -127,7 +127,7 @@ oprob[k1]
 ```
 To interface with integrators and solutions we use a similar syntax.
 
-Finally, instead of using `@unpack` to access a symbolic variable or parameter, we can access it directly using `rn.X1`, and thus access a state of our `ODEProblem` using
+Finally, instead of using `@unpack` to access a symbolic variable or parameter, we can access it directly using `rn.X1`, and thus access an unknown of our `ODEProblem` using
 ```@example ex2
 oprob[rn.X1]
 ```

docs/src/catalyst_functionality/compositional_modeling.md

@@ -72,7 +72,7 @@ Catalyst.get_rxs(rn)
 ```
 To see all the species, parameters and reactions in the tree we can use
 ```@example ex1
-species(rn)   # or states(rn)
+species(rn)   # or unknowns(rn)
 ```
 ```@example ex1
 parameters(rn)  # or reactionparameters(rn)

docs/src/catalyst_functionality/constraint_equations.md

@@ -47,7 +47,7 @@ rn = @reaction_network begin
 end
 ```
 Notice, here we interpolated the variable `V` with `$V` to ensure we use the
-same symbolic state variable in the `rn` as we used in building `osys`. See the
+same symbolic unknown variable in the `rn` as we used in building `osys`. See the
 doc section on [interpolation of variables](@ref
 dsl_description_interpolation_of_variables) for more information.
 
@@ -99,7 +99,7 @@ lower-level approach for creating events via the DifferentialEquations.jl
 callback interface is illustrated in the [Advanced Simulation Options](@ref
 advanced_simulations) tutorial.
 
-Let's first create our equations and states/species again
+Let's first create our equations and unknowns/species again
 ```@example ceq3
 using Catalyst, DifferentialEquations, Plots
 
@@ -124,7 +124,7 @@ oprob = ODEProblem(rs, [], (0.0, 10.0))
 sol = solve(oprob, Tsit5())
 plot(sol; plotdensity = 1000)
 ```
-We can also model discrete events. Similar to our example with continuous events, we start by creating reaction equations, parameters, variables, and states. 
+We can also model discrete events. Similar to our example with continuous events, we start by creating reaction equations, parameters, variables, and unknowns. 
 ```@example ceq3
 @parameters k_on switch_time k_off
 @variables t

docs/src/catalyst_functionality/dsl_description.md

@@ -546,7 +546,7 @@ end
 Please see the API [Rate Laws](@ref api_rate_laws) section for more details.
 
 ## Including non-species variables
-Non-species state variables can be specified in the DSL using the `@variables`
+Non-species unknown variables can be specified in the DSL using the `@variables`
 macro. These are declared similarly to species. For example,
 ```@example tut2
 rn_with_volume = @reaction_network begin
@@ -562,10 +562,10 @@ and one non-species
 ```@example tut2
 nonspecies(rn_with_volume)
 ```
-giving two state variables, always internally ordered by species and then
+giving two unknown variables, always internally ordered by species and then
 nonspecies:
 ```@example tut2
-states(rn_with_volume)
+unknowns(rn_with_volume)
 ```
 
 `rn_with_volume` could then be extended with constraint equations for how `V(t)`
@@ -584,9 +584,9 @@ rn_with_s = @reaction_network begin
 end
 show(stdout, MIME"text/plain"(), rn_with_s)  # hide
 ```
-where we see all states are now functions of `s`.
+where we see all unknowns are now functions of `s`.
 
-Similarly, if one wants states to be functions of more than one independent
+Similarly, if one wants unknowns to be functions of more than one independent
 variable, for example to encode a spatial problem, one can list more than one
 variable, i.e. `@ivs t x y`. Here the first listed independent variable is
 always chosen to represent time. For example,

docs/src/catalyst_functionality/example_networks/smoluchowski_coagulation_equation.md

@@ -51,7 +51,7 @@ end
 ```
 We'll store the reaction rates in `pars` as `Pair`s, and set the initial condition that only monomers are present at ``t=0`` in `u₀map`.
 ```julia
-# state variables are X, pars stores rate parameters for each rx
+# unknown variables are X, pars stores rate parameters for each rx
 @variables t
 @species k[1:nr] (X(t))[1:N]
 pars = Pair.(collect(k), kv)