Remake for conservation laws update

HISTORY.md

@@ -13,6 +13,9 @@
   solver that auto-switches between explicit and implicit methods, install `OrdinaryDiffEqDefault`. To 
   use `Tsit5` install `OrdinaryDiffEqTsit5`, etc. The possible sub-libraries, each containing different solvers,
   can be viewed [here](https://github.com/SciML/OrdinaryDiffEq.jl/tree/master/lib).
+- It should now be safe to use `remake` on problems which have had conservation laws removed.
+  The warning regarding this when eliminating conservation laws have been removed (in conjunction
+  with the `remove_conserved_warn` argument for disabling this warning).
 - New formula for inferring variables from equations (declared using the `@equations` options) in the DSL. The order of inference of species/variables/parameters is now:
     (1) Every symbol explicitly declared using `@species`, `@variables`, and `@parameters` are assigned to the correct category.
     (2) Every symbol used as a reaction reactant is inferred as a species.

docs/src/model_creation/conservation_laws.md

@@ -1,5 +1,5 @@
 # [Working with Conservation Laws](@id conservation_laws)
-Catalyst can detect, and eliminate for differential-equation based models, *conserved quantities*, i.e. linear combinations of species which are conserved via the chemistry. This functionality is both automatically utilised by Catalyst (e.g. to [remove singular Jacobians during steady state computations](@ref homotopy_continuation_conservation_laws)), but is also available for users to utilise directly (e.g. to potentially [improve simulation performance](@ref ode_simulation_performance_conservation_laws)).
+Catalyst can detect, and eliminate for differential-equation based models, *conserved quantities*, i.e. linear combinations of species which are conserved via their chemistry. This functionality is both automatically utilised by Catalyst (e.g. to [remove singular Jacobians during steady state computations](@ref homotopy_continuation_conservation_laws)), but is also available for users to utilise directly (e.g. to potentially [improve simulation performance](@ref ode_simulation_performance_conservation_laws)).
 
 To illustrate conserved quantities, let us consider the following [two-state](@ref basic_CRN_library_two_states) model:
 ```@example conservation_laws
@@ -39,7 +39,7 @@ Using the `unknowns` function we can confirm that the ODE only has a single unkn
 ```@example conservation_laws
 unknowns(osys)
 ```
-Next, using `parameters` we note that an additional parameter, `Γ[1]` has been added to the system:
+Next, using `parameters` we note that an additional (vector) parameter, `Γ` has been added to the system:
 ```@example conservation_laws
 parameters(osys)
 ```
@@ -53,27 +53,27 @@ ps = [:k₁ => 10.0, :k₂ => 2.0]
 oprob = ODEProblem(rs, u0, (0.0, 1.0), ps; remove_conserved = true)
 nothing # hide
 ```
-Here, while `Γ[1]` becomes a parameter of the new system, it has a [default value](@ref dsl_advanced_options_default_vals) equal to the corresponding conservation law. Hence, its value is computed from the initial condition `[:X₁ => 80.0, :X₂ => 20.0]`, and does not need to be provided in the parameter vector. Next, we can simulate and plot our model using normal syntax:
+Here, while `Γ` becomes a parameter of the new system, it has a [default value](@ref dsl_advanced_options_default_vals) equal to the corresponding conservation law. Hence, its value is computed from the initial condition `[:X₁ => 80.0, :X₂ => 20.0]`, and does not need to be provided in the parameter vector. Next, we can simulate and plot our model using normal syntax:
 ```@example conservation_laws
 sol = solve(oprob)
 plot(sol)
 ```
 !!! note
     Any species eliminated using `remove_conserved = true` will not be plotted by default. However, it can be added to the plot using [the `idxs` plotting option](@ref simulation_plotting_options). E.g. here we would use `plot(sol; idxs = [:X₁, :X₂])` to plot both species.
 
-!!! danger
-    Currently, there is a bug in MTK where the parameters, `Γ`, associated with conservation laws are not updated properly in response to [`remake`](@ref simulation_structure_interfacing_problems_remake) (or [other problem-updating functions, such as `getu`](@ref simulation_structure_interfacing_functions)). Hence, problems created using `remove_conserved = true` *should not* be modified, i.e. by changing initial conditions or the values of the `Γ`'s directly. Instead, to change these values new problems must be generated with the new initial condition map.
-
 While `X₂` is an observable (and not unknown) of the ODE, we can [access it](@ref simulation_structure_interfacing_problems) just like if `remove_conserved = true` had not been used:
 ```@example conservation_laws
 sol[:X₂]
 ```
 !!! note
-    Generally, `remove_conserved = true` should not change any model workflows. I.e. anything that works without this option should also work when an `ODEProblem` is created using `remove_conserved = true`.
+    Generally, `remove_conserved = true` should not change any modelling workflows. I.e. anything that works without this option should also work when an `ODEProblem` is created using `remove_conserved = true`.
 
 !!! note
     The `remove_conserved = true` option is available when creating `SDEProblem`s, `NonlinearProblem`s, and `SteadyStateProblem`s (and their corresponding systems). However, it cannot be used when creating `JumpProblem`s.
 
+!!! warn
+    Users of the [ModelingToolkit.jl](https://github.com/SciML/ModelingToolkit.jl) package might be familiar with the `structural_simplify` function. While it can be applied to Catalyst models as well, generally, this should be avoided (as `remove_conserved` performs a similar role, but is better adapted to these models). Furthermore, applying `structural_simplify` will interfere with conservation law removal, preventing users from accessing eliminated quantities.
+
 ## [Conservation law accessor functions](@id conservation_laws_accessors)
 
 For any given `ReactionSystem` model, we can use `conservationlaw_constants` to compute all of a system's conserved quantities:
@@ -89,3 +89,54 @@ Finally, the `conservationlaws` function yields a $m \times n$ matrix, where $n$
 conservationlaws(rs)
 ```
 I.e. in this case we have a single conserved quantity, which contains a single copy each of the system's two species.
+
+## [Updating conservation law values directly](@id conservation_laws_prob_updating)
+Previously we noted that conservation law elimination adds the conservation law as a system parameter (named $Γ$). E.g. here we find it among the parameters of the ODE model generated by the previously used two-state system:
+```@example conservation_laws
+parameters(convert(ODESystem, rs; remove_conserved = true))
+```
+An advantage of this is that we can set conserved quantities's values directly in simulations. Here we simulate the model while specifying that $X₁ + X₂ = 10.0$ (and while also specifying an initial condition for $X₁$, but none for $X₂$).
+```@example conservation_laws
+u0 = [:X₁ => 6.0]
+ps = [:k₁ => 1.0, :k₂ => 2.0, :Γ => [10.0]]
+oprob = ODEProblem(rs, u0, 1.0, ps; remove_conserved = true)
+sol = solve(oprob)
+nothing # hide
+```
+Generally, for each conservation law, one can omit specifying either the conservation law constant, or one initial condition (whichever quantity is missing will be computed from the other ones). 
+!!! note
+    As previously mentioned, the conservation law parameter $Γ$ is a [*vector-valued* parameter](@ref dsl_advanced_options_vector_variables) with one value for each conservation law. That is why we above provide its value as a vector (`:Γ => [5.0]`). If there had been multiple conservation laws, we would provide the `:Γ` vector with one value for each of them (e.g. `:Γ => [10.0, 15.0]`).
+!!! warn
+    If you specify the value of a conservation law parameter, you *must not* specify the value of all species of that conservation law (this can result in an error). Instead, the value of exactly one species must be left unspecified.
+
+Just like when we create a problem, if we [update the species (or conservation law parameter) values of `oprob`](@ref simulation_structure_interfacing_problems), the remaining ones will be recomputed to generate an accurate conservation law. E.g. here we create an `ODEProblem`, check the value of the conservation law, and then confirm that its value is updated with $X₁$.
+```@example conservation_laws
+u0 = [:X₁ => 6.0, :X₂ => 4.0]
+ps = [:k₁ => 1.0, :k₂ => 2.0]
+oprob = ODEProblem(rs, u0, 10.0, ps; remove_conserved = true)
+oprob.ps[:Γ][1]
+```
+```@example conservation_laws
+oprob = remake(oprob; u0 = [:X₁ => 16.0])
+oprob.ps[:Γ][1]
+```
+It is also possible to update the value of $Γ$. Here, as $X₂$ is the species eliminated by the conservation law (which can be checked using `conservedequations`), $X₂$'s value will be modified to ensure that $Γ$'s new value is correct:
+```@example conservation_laws
+oprob = remake(oprob; p = [:Γ => [30.0]])
+oprob[:X₂]
+```
+
+Generally, for a conservation law where we have:
+- One (or more) non-eliminated species.
+- One eliminated species.
+- A conservation law parameter.
+If the value of the conservation law parameter $Γ$'s value *has never been specified*, its value will be updated to accommodate changes in species values. If the conservation law parameter ($Γ$)'s value *has been specified* (either when the `ODEProblem` was created, or using `remake`), then the eliminated species's value will be updated to accommodate changes in the conservation law parameter or non-eliminated species's values. Furthermore, in this case, the value of the eliminated species *cannot be updated*. 
+
+!!! warn
+    When updating the values of problems with conservation laws it is additionally important to use `remake` (as opposed to direct indexing, e.g. setting `oprob[:X₁] = 16.0`).
+
+### [Extracting the conservation law parameter's symbolic variable](@id conservation_laws_prob_updating_symvar)
+Catalyst represents its models using the [Symbolics.jl](https://github.com/JuliaSymbolics/Symbolics.jl) computer algebraic system, something which allows the user to [form symbolic expressions of model components](@ref simulation_structure_interfacing_symbolic_representation). If you wish to extract and use the symbolic variable corresponding to a model's conserved quantities, you can use the following syntax:
+```@example conservation_laws
+Γ = Catalyst.get_networkproperties(rs).conservedconst
+```

ext/CatalystBifurcationKitExtension/bifurcation_kit_extension.jl

@@ -37,7 +37,6 @@ function bkext_make_nsys(rs, u0)
     cons_default = [cons_eq.rhs for cons_eq in cons_eqs]
     cons_default = Catalyst.get_networkproperties(rs).conservedconst => cons_default
     defaults = Dict([u0; cons_default])
-    nsys = convert(NonlinearSystem, rs; defaults,
-        remove_conserved = true, remove_conserved_warn = false)
+    nsys = convert(NonlinearSystem, rs; defaults, remove_conserved = true)
     return complete(nsys)
 end

ext/CatalystHomotopyContinuationExtension/homotopy_continuation_extension.jl

@@ -51,8 +51,7 @@ function steady_state_polynomial(rs::ReactionSystem, ps, u0)
     # Creates the appropriate nonlinear system, and converts parameters to a form that can
     # be substituted in later.
     rs = Catalyst.expand_registered_functions(rs)
-    ns = complete(convert(NonlinearSystem, rs;
-        remove_conserved = true, remove_conserved_warn = false))
+    ns = complete(convert(NonlinearSystem, rs; remove_conserved = true))
     pre_varmap = [symmap_to_varmap(rs, u0)..., symmap_to_varmap(rs, ps)...]
     Catalyst.conservationlaw_errorcheck(rs, pre_varmap)
     p_dict = make_p_val_dict(pre_varmap, rs, ns)

ext/CatalystStructuralIdentifiabilityExtension/structural_identifiability_extension.jl

@@ -169,7 +169,7 @@ function make_osys(rs::ReactionSystem; remove_conserved = true)
         error("Identifiability should only be computed for complete systems. A ReactionSystem can be marked as complete using the `complete` function.")
     end
     rs = complete(Catalyst.expand_registered_functions(flatten(rs)))
-    osys = complete(convert(ODESystem, rs; remove_conserved, remove_conserved_warn = false))
+    osys = complete(convert(ODESystem, rs; remove_conserved))
     vars = [unknowns(rs); parameters(rs)]
 
     # Computes equations for system conservation laws.