Merge pull request #1206 from SciML/fix_nlsys_conslaws

Re-enable conservation laws for `NonlinearSystem`s

Author: isaacsas

HISTORY.md

@@ -7,29 +7,47 @@
   (at the time the release is made). If you need a dependency version increased,
   please open an issue and we can update it and make a new Catalyst release once
   testing against the newer dependency version is complete.
-- It is no longer recommended to install and use the full OrdinaryDiffEq library to access specific ODE solvers.
-  Instead, only install the specific OrdinaryDiffEq sub-libraries that contain the desired
-  solver. This significantly reduces installation and package loading times. I.e. to use the default
-  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.
-    (3) Every symbol not declared in (1) or (2) that occurs in an expression provided after `@equations` is inferred as a variable.
-    (4) Every symbol not declared in (1), (2), or (3) that occurs either as a reaction rate or stoichiometric coefficient is inferred to be a parameter.
-  E.g. in
-  ```julia
-  @reaction_network begin
-      @equations V1 + S ~ V2^2
-      (p + S + V1), S --> 0
-  end
-  ```
-  `S` is inferred as a species, `V1` and `V2` as variables, and `p` as a parameter. The previous special cases for the `@observables`, `@compounds`, and `@differentials` options still hold. Finally, the `@require_declaration` options (described in more detail below) can now be used to require everything to be explicitly declared.
-- New formula for determining whether the default differentials have been used within an `@equations` option. Now, if any expression `D(...)` is encountered (where `...` can be anything), this is inferred as usage of the default differential D. E.g. in the following equations `D` is inferred as a differential with respect to the default independent variable:
+- It is no longer recommended to install and use the full OrdinaryDiffEq library
+  to access specific ODE solvers. Instead, only install the specific
+  OrdinaryDiffEq sub-libraries that contain the desired solver. This
+  significantly reduces installation and package loading times. I.e. to use the
+  default 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 with the exception of `NonlinearProblem`s or `NonlinearSystem`s.
+  For NonlinearProblems it is safe to use `remake` if only updating `u0` values,
+  but it is not safe to update the value of the conserved constant, `Γ`. See
+  [the FAQ](https://docs.sciml.ai/Catalyst/stable/faqs/#faq_remake_nonlinprob)
+  for details.
+- 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.
+    3. Every symbol not declared in (1) or (2) that occurs in an expression
+       provided after `@equations` is inferred as a variable.
+    4. Every symbol not declared in (1), (2), or (3) that occurs either as a
+       reaction rate or stoichiometric coefficient is inferred to be a
+       parameter. E.g. in
+       ```julia
+       @reaction_network begin
+           @equations V1 + S ~ V2^2
+           (p + S + V1), S --> 0
+       end
+       ```
+       `S` is inferred as a species, `V1` and `V2` as variables, and `p` as a
+       parameter. The previous special cases for the `@observables`, `@compounds`,
+       and `@differentials` options still hold. Finally, the `@require_declaration`
+       options (described in more detail below) can now be used to require everything
+       to be explicitly declared.
+- New formula for determining whether the default differentials have been used
+  within an `@equations` option. Now, if any expression `D(...)` is encountered
+  (where `...` can be anything), this is inferred as usage of the default
+  differential D. E.g. in the following equations `D` is inferred as a
+  differential with respect to the default independent variable:
   ```julia
   @reaction_network begin
     @equations D(V) + V ~ 1
@@ -38,7 +56,9 @@
     @equations D(D(V)) ~ 1
   end
   ```
-  Please note that this cannot be used at the same time as `D` is used to represent a species, variable, or parameter (including if these are implicitly designated as such by e.g. appearing as a reaction reactant). 
+  Please note that this cannot be used at the same time as `D` is used to
+  represent a species, variable, or parameter (including if these are implicitly
+  designated as such by e.g. appearing as a reaction reactant). 
 - Array symbolics support is more consistent with ModelingToolkit v9. Parameter
   arrays are no longer scalarized by Catalyst, while species and variables
   arrays still are (as in ModelingToolkit). As such, parameter arrays should now
@@ -51,7 +71,11 @@
   *not* to do this as it has signifcant performance costs with ModelingToolkit
   v9. Note, scalarized parameter arrays passed to the two-argument
   `ReactionSystem` constructor may become unscalarized.
-- Functional (e.g. time-dependent) parameters can now be used in Catalyst models. These can e.g. be used to incorporate arbitrary time-dependent functions (as a parameter) in a model. For more details on how to use these, please read: https://docs.sciml.ai/Catalyst/stable/model_creation/functional_parameters/.
+- Functional (e.g. time-dependent) parameters can now be used in Catalyst
+  models. These can e.g. be used to incorporate arbitrary time-dependent
+  functions (as a parameter) in a model. For more details on how to use these,
+  please read:
+  https://docs.sciml.ai/Catalyst/stable/model_creation/functional_parameters/.
 - Scoped species/variables/parameters are now treated similar to the latest MTK
   releases (≥ 9.49).
 - A tutorial on making interactive plot displays using Makie has been added.

docs/src/faqs.md

@@ -271,7 +271,7 @@ so the reaction
 ```@example faq7
 using Catalyst
 rn = @reaction_network begin
-  k, X --> ∅
+    k, X --> ∅
 end
 convert(ODESystem, rn)
 ```
@@ -280,7 +280,7 @@ using any of the following non-filled arrows when declaring the reaction: `<=`,
 `⇐`, `⟽`, `=>`, `⇒`, `⟾`, `⇔`, `⟺`. This means that the reaction
 ```@example faq7
 rn = @reaction_network begin
-  k, X => ∅
+    k, X => ∅
 end
 convert(ODESystem, rn)
 ```
@@ -290,7 +290,7 @@ will be degraded at a constant rate even when very small or equal to 0).
 Note, stoichiometric coefficients are still included, i.e. the reaction
 ```@example faq7
 rn = @reaction_network begin
-  k, 2*X ⇒ ∅
+    k, 2*X ⇒ ∅
 end
 convert(ODESystem, rn)
 ```
@@ -303,7 +303,7 @@ ModelingToolkit. e.g., this is should work
 using Catalyst
 myHill(x) = 2*x^3/(x^3+1.5^3)
 rn = @reaction_network begin
-  myHill(X), ∅ --> X
+    myHill(X), ∅ --> X
 end
 ```
 In some cases, it may be necessary or desirable to register functions with
@@ -322,3 +322,137 @@ rn = @reaction_network begin
     (k1, k2), A <--> B
 end
 ```
+
+## [What to be aware of when using `remake` with conservation law elimination and NonlinearProblems?](@id faq_remake_nonlinprob)
+
+When constructing `NonlinearSystem`s or `NonlinearProblem`s with `remove_conserved = true`, i.e.
+```julia
+# for rn a ReactionSystem
+nsys = convert(NonlinearSystem, rn; remove_conserved = true)
+
+# or 
+nprob = NonlinearProblem(rn, u0, p; remove_conserved = true)
+```
+`remake` is currently unable to correctly update all `u0` values when the
+conserved constant(s), `Γ`, are updated. As an example consider the following
+```@example faq_remake
+using Catalyst, NonlinearSolve
+rn = @reaction_network begin
+    (k₁,k₂), X₁ <--> X₂
+    (k₃,k₄), X₁ + X₂ --> 2X₃
+end
+u0 = [:X₁ => 1.0, :X₂ => 2.0, :X₃ => 3.0]
+ps = [:k₁ => 0.1, :k₂ => 0.2, :k₃ => 0.3, :k₄ => 0.4]
+nlsys = convert(NonlinearSystem, rn; remove_conserved = true, conseqs_remake_warn = false)
+nlsys = complete(nlsys)
+equations(nlsys)
+```
+If we generate a `NonlinearProblem` from this system the conservation constant,
+`Γ[1]`, is automatically set to `X₁ + X₂ + X₃ = 6` and the initial values are
+those in `u0`. i.e if
+```@example faq_remake
+nlprob1 = NonlinearProblem(nlsys, u0, ps)
+```
+then
+```@example faq_remake
+nlprob1[(:X₁, :X₂, :X₃)] == (1.0, 2.0, 3.0)
+```
+and
+```@example faq_remake
+nlprob1.ps[:Γ][1] == 6.0
+```
+If we now try to change a value of `X₁`, `X₂`, or `X₃` using `remake`, the
+conserved constant will be recalculated. i.e. if
+```@example faq_remake
+nlprob2 = remake(nlprob1; u0 = [:X₂ => 3.0])
+```
+compare  
+```@example faq_remake
+println("Correct u0 is: ", (1.0, 3.0, 3.0), "\n", "remade value is: ", nlprob2[(:X₁, :X₂, :X₃)]) 
+```
+and
+```@example faq_remake
+println("Correct Γ is: ", 7.0, "\n", "remade value is: ", nlprob2.ps[:Γ][1])
+```
+However, if we try to directly change the value of `Γ` it is not always the case
+that a `u0` value will correctly update so that the conservation law is
+conserved. Consider
+```@example faq_remake
+nlprob3 = remake(nlprob1; u0 = [:X₂ => nothing], p = [:Γ => [4.0]])
+```
+Setting `[:X₂ => nothing]` for other problem types communicates that the
+`u0` value for `X₂` should be solved for. However, if we examine the values we
+find
+```@example faq_remake
+println("Correct u0 is: ", (1.0, 0.0, 3.0), "\n", "remade value is: ", nlprob3[(:X₁, :X₂, :X₃)]) 
+```
+and
+```@example faq_remake
+println("Correct Γ is: ", 4.0, "\n", "remade value is: ", nlprob3.ps[:Γ][1])
+```
+As such, the `u0` value for `X₂` has not updated, and the conservation law is
+now violated by the `u0` values, i.e,
+```@example faq_remake
+(nlprob3[:X₁] + nlprob3[:X₂] + nlprob3[:X₃]) == nlprob3.ps[:Γ][1] 
+```
+Currently, the only way to avoid this issue is to manually specify updated
+values for the `u0` components, which will ensure that `Γ` updates appropriately
+as in the first example. i.e. we manually set `X₂` to the value it should be and
+`Γ` will be updated accordingly:
+```@example faq_remake
+nlprob4 = remake(nlprob1; u0 = [:X₂ => 0.0])
+```
+so that
+```@example faq_remake
+println("Correct u0 is: ", (1.0, 0.0, 3.0), "\n", "remade value is: ", nlprob4[(:X₁, :X₂, :X₃)]) 
+```
+and
+```@example faq_remake
+println("Correct Γ is: ", 4.0, "\n", "remade value is: ", nlprob4.ps[:Γ][1])
+```
+
+Finally, we note there is one extra consideration to take into account if using
+`structural_simplify`. In this case one of `X₁`, `X₂`, or `X₃` will be moved to
+being an observed. It will then always correspond to the updated value if one
+tries to manually change `Γ`. Let's see what happens here directly
+```@example faq_remake
+nlsys = convert(NonlinearSystem, rn; remove_conserved = true, conseqs_remake_warn = false)
+nlsys = structural_simplify(nlsys)
+nlprob1 = NonlinearProblem(nlsys, u0, ps)
+```
+We can now try to change just `Γ` and implicitly the observed variable that was
+removed will be assumed to have changed its initial value to compensate for it.
+Let's confirm this. First we find the observed variable that was elminated.
+```@example faq_remake
+obs_unknown = only(observed(nlsys)).lhs
+```
+We can figure out its index in `u0` via
+```@example faq_remake
+obs_symbol = ModelingToolkit.getname(obs_unknown)
+obsidx = findfirst(p -> p[1] == obs_symbol, u0)
+```
+Let's now remake 
+```@example faq_remake
+nlprob2 = remake(nlprob1; u0 = [obs_unknown => nothing], p = [:Γ => [8.0]])
+```
+Here we indicate that the observed variable should be treated as unspecified
+during initialization. Since the observed variable is not considered an unknown,
+everything now works, with the observed variable's assumed initial value
+adjusted to allow `Γ = 8`:
+```@example faq_remake
+correct_u0 = last.(u0)
+correct_u0[obsidx] = 8 - sum(correct_u0) + correct_u0[obsidx]
+println("Correct u0 is: ", (1.0, 2.0, 5.0), "\n", "remade value is: ", nlprob2[(:X₁, :X₂, :X₃)]) 
+```
+and `Γ` becomes
+```@example faq_remake
+println("Correct Γ is: ", 8.0, "\n", "remade value is: ", nlprob2.ps[:Γ][1])
+```
+Unfortunately, as with our first example, trying to enforce that a
+non-eliminated species should have its initial value updated instead of the
+observed species will not work.
+
+*Summary:* it is not recommended to directly update `Γ` via `remake`, but to
+instead update values of the initial guesses in `u0` to obtain a desired `Γ`. At
+this time the behavior when updating `Γ` can result in `u0` values that do not
+satisfy the conservation law defined by `Γ` as illustrated above. 

docs/src/model_creation/conservation_laws.md

@@ -133,7 +133,7 @@ Generally, for a conservation law where we have:
 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`).
+    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`). Moreover, care is needed when using `remake` with `NonlinearProblem`s for which conserved equations have been removed. See [the FAQ](https://docs.sciml.ai/Catalyst/stable/faqs/#faq_remake_nonlinprob) for details on what is supported and what may lead to `u0` values that do not satisfy the conservation law in the special case of `NonlinearProblem`s.
 
 ### [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:

ext/CatalystBifurcationKitExtension/bifurcation_kit_extension.jl

@@ -37,6 +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)
+    nsys = convert(NonlinearSystem, rs; defaults, remove_conserved = true, conseqs_remake_warn = false)
     return complete(nsys)
 end

ext/CatalystHomotopyContinuationExtension/homotopy_continuation_extension.jl

@@ -1,16 +1,16 @@
 ### Homotopy Continuation Based Steady State Finding ###
 
 """
-    hc_steady_states(rs::ReactionSystem, ps; filter_negative=true, neg_thres=-1e-20, u0=typeof(ps)(), kwargs...)
+    hc_steady_states(rs::ReactionSystem, ps; filter_negative = true, neg_thres = -1e-16, u0 = typeof(ps)(), kwargs...)
 
 Uses homotopy continuation via HomotopyContinuation.jl to find the steady states of the ODE system corresponding to the provided reaction system.
 
 Arguments:
 - `rs::ReactionSystem`: The reaction system for which we want to find the steady states.
 - `ps`: The parameter values for which we want to find the steady states.
-- `filter_negative=true`: If set to true, solutions with any species concentration <neg_thres is removed from the output.
-- `neg_thres=-1e-20`: Determine the minimum values for which a species concentration is to be considered non-negative. Species concentrations ``> neg_thres`` but `< 0.0` are set to `0.0`.
-- `u0=nothing`: Initial conditions for which we want to find the steady states. For systems with conservation laws this are required to compute conserved quantities. Initial conditions are not required for all species, only those involved in conserved quantities (if this set is unknown, it is recommended to provide initial conditions for all species).
+- `filter_negative = true`: If set to true, solutions with any species concentration <neg_thres is removed from the output.
+- `neg_thres = -1e-15`: Determine the minimum values for which a species concentration is to be considered non-negative. Species concentrations ``> neg_thres`` but `< 0.0` are set to `0.0`.
+- `u0 = nothing`: Initial conditions for which we want to find the steady states. For systems with conservation laws this are required to compute conserved quantities. Initial conditions are not required for all species, only those involved in conserved quantities (if this set is unknown, it is recommended to provide initial conditions for all species).
 - `kwargs...`: any additional arguments (like `show_progress= true`) are passed into HomotopyContinuation.jl's `solve` call.
 
 Examples
@@ -36,7 +36,7 @@ Notes:
 ```
   """
 function Catalyst.hc_steady_states(rs::ReactionSystem, ps; filter_negative = true,
-        neg_thres = -1e-20, u0 = [], kwargs...)
+        neg_thres = -1e-15, u0 = [], kwargs...)
     if !isautonomous(rs)
         error("Attempting to compute steady state for a non-autonomous system (e.g. where some rate depend on $(get_iv(rs))). This is not possible.")
     end
@@ -51,7 +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))
+    ns = complete(convert(NonlinearSystem, rs; remove_conserved = true, conseqs_remake_warn = false))
     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)
@@ -164,7 +164,7 @@ function reorder_sols!(sols, ss_poly, rs::ReactionSystem)
 end
 
 # Filters away solutions with negative species concentrations (and for neg_thres < val < 0.0, sets val=0.0).
-function filter_negative_f(sols; neg_thres = -1e-20)
+function filter_negative_f(sols; neg_thres = -1e-15)
     for sol in sols, idx in 1:length(sol)
         (neg_thres < sol[idx] < 0) && (sol[idx] = 0)
     end