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 conjecture
+  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

@@ -61,9 +61,6 @@ 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₂]

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.

src/reactionsystem_conversions.jl

@@ -449,21 +449,6 @@ diff_2_zero(expr) = (Symbolics.is_derivative(expr) ? 0 : expr)
 
 COMPLETENESS_ERROR = "A ReactionSystem must be complete before it can be converted to other system types. A ReactionSystem can be marked as complete using the `complete` function."
 
-# Used to, when required, display a warning about conservation law removal and remake.
-function check_cons_warning(remove_conserved, remove_conserved_warn)
-    (remove_conserved && remove_conserved_warn) || return
-    @warn "You are creating a system or problem while eliminating conserved quantities. Please note,
-        due to limitations / design choices in ModelingToolkit if you use the created system to
-        create a problem (e.g. an `ODEProblem`), or are directly creating a problem, you *should not*
-        modify that problem's initial conditions for species (e.g. using `remake`). Changing initial
-        conditions must be done by creating a new Problem from your reaction system or the
-        ModelingToolkit system you converted it into with the new initial condition map.
-        Modification of parameter values is still possible, *except* for the modification of any
-        conservation law constants ($CONSERVED_CONSTANT_SYMBOL), which is not possible. You might
-        get this warning when creating a problem directly.
-
-        You can remove this warning by setting `remove_conserved_warn = false`."
-end
 
 ### System Conversions ###
 
@@ -482,19 +467,16 @@ Keyword args and default values:
 - `remove_conserved=false`, if set to `true` will calculate conservation laws of the
   underlying set of reactions (ignoring constraint equations), and then apply them to reduce
   the number of equations.
-- `remove_conserved_warn = true`: If `true`, if also `remove_conserved = true`, there will be
-  a warning regarding limitations of modifying problems generated from the created system.
 """
 function Base.convert(::Type{<:ODESystem}, rs::ReactionSystem; name = nameof(rs),
         combinatoric_ratelaws = get_combinatoric_ratelaws(rs),
-        include_zero_odes = true, remove_conserved = false, remove_conserved_warn = true,
+        include_zero_odes = true, remove_conserved = false,
         checks = false, default_u0 = Dict(), default_p = Dict(),
         defaults = _merge(Dict(default_u0), Dict(default_p)),
         kwargs...)
     # Error checks.
     iscomplete(rs) || error(COMPLETENESS_ERROR)
     spatial_convert_err(rs::ReactionSystem, ODESystem)
-    check_cons_warning(remove_conserved, remove_conserved_warn)
 
     fullrs = Catalyst.flatten(rs)
     remove_conserved && conservationlaws(fullrs)
@@ -529,19 +511,15 @@ Keyword args and default values:
 - `remove_conserved=false`, if set to `true` will calculate conservation laws of the
   underlying set of reactions (ignoring constraint equations), and then apply them to reduce
   the number of equations.
-- `remove_conserved_warn = true`: If `true`, if also `remove_conserved = true`, there will be
-  a warning regarding limitations of modifying problems generated from the created system.
 """
 function Base.convert(::Type{<:NonlinearSystem}, rs::ReactionSystem; name = nameof(rs),
         combinatoric_ratelaws = get_combinatoric_ratelaws(rs),
         include_zero_odes = true, remove_conserved = false, checks = false,
-        remove_conserved_warn = true, default_u0 = Dict(), default_p = Dict(),
-        defaults = _merge(Dict(default_u0), Dict(default_p)),
+        default_u0 = Dict(), default_p = Dict(), defaults = _merge(Dict(default_u0), Dict(default_p)),
         all_differentials_permitted = false, kwargs...)
     # Error checks.
     iscomplete(rs) || error(COMPLETENESS_ERROR)
     spatial_convert_err(rs::ReactionSystem, NonlinearSystem)
-    check_cons_warning(remove_conserved, remove_conserved_warn)
     if !isautonomous(rs)
         error("Attempting to convert a non-autonomous `ReactionSystem` (e.g. where some rate depend on $(get_iv(rs))) to a `NonlinearSystem`. This is not possible. if you are intending to compute system steady states, consider creating and solving a `SteadyStateProblem.")
     end
@@ -612,19 +590,15 @@ Notes:
 - `remove_conserved=false`, if set to `true` will calculate conservation laws of the
   underlying set of reactions (ignoring constraint equations), and then apply them to reduce
   the number of equations.
-- `remove_conserved_warn = true`: If `true`, if also `remove_conserved = true`, there will be
-  a warning regarding limitations of modifying problems generated from the created system.
 """
 function Base.convert(::Type{<:SDESystem}, rs::ReactionSystem;
         name = nameof(rs), combinatoric_ratelaws = get_combinatoric_ratelaws(rs),
         include_zero_odes = true, checks = false, remove_conserved = false,
-        remove_conserved_warn = true, default_u0 = Dict(), default_p = Dict(),
-        defaults = _merge(Dict(default_u0), Dict(default_p)),
+        default_u0 = Dict(), default_p = Dict(), defaults = _merge(Dict(default_u0), Dict(default_p)),
         kwargs...)
     # Error checks.
     iscomplete(rs) || error(COMPLETENESS_ERROR)
     spatial_convert_err(rs::ReactionSystem, SDESystem)
-    check_cons_warning(remove_conserved, remove_conserved_warn)
 
     flatrs = Catalyst.flatten(rs)
 
@@ -708,12 +682,12 @@ function DiffEqBase.ODEProblem(rs::ReactionSystem, u0, tspan,
         p = DiffEqBase.NullParameters(), args...;
         check_length = false, name = nameof(rs),
         combinatoric_ratelaws = get_combinatoric_ratelaws(rs),
-        include_zero_odes = true, remove_conserved = false, remove_conserved_warn = true,
-        checks = false, structural_simplify = false, kwargs...)
+        include_zero_odes = true, remove_conserved = false, checks = false, 
+        structural_simplify = false, kwargs...)
     u0map = symmap_to_varmap(rs, u0)
     pmap = symmap_to_varmap(rs, p)
     osys = convert(ODESystem, rs; name, combinatoric_ratelaws, include_zero_odes, checks,
-        remove_conserved, remove_conserved_warn)
+        remove_conserved)
 
     # Handles potential differential algebraic equations (which requires `structural_simplify`).
     if structural_simplify
@@ -732,12 +706,12 @@ function DiffEqBase.NonlinearProblem(rs::ReactionSystem, u0,
         p = DiffEqBase.NullParameters(), args...;
         name = nameof(rs), include_zero_odes = true,
         combinatoric_ratelaws = get_combinatoric_ratelaws(rs),
-        remove_conserved = false, remove_conserved_warn = true, checks = false,
-        check_length = false, all_differentials_permitted = false, kwargs...)
+        remove_conserved = false, checks = false, check_length = false, 
+        all_differentials_permitted = false, kwargs...)
     u0map = symmap_to_varmap(rs, u0)
     pmap = symmap_to_varmap(rs, p)
     nlsys = convert(NonlinearSystem, rs; name, combinatoric_ratelaws, include_zero_odes,
-        checks, all_differentials_permitted, remove_conserved, remove_conserved_warn)
+        checks, all_differentials_permitted, remove_conserved)
     nlsys = complete(nlsys)
     return NonlinearProblem(nlsys, u0map, pmap, args...; check_length,
         kwargs...)
@@ -748,11 +722,11 @@ function DiffEqBase.SDEProblem(rs::ReactionSystem, u0, tspan,
         p = DiffEqBase.NullParameters(), args...;
         name = nameof(rs), combinatoric_ratelaws = get_combinatoric_ratelaws(rs),
         include_zero_odes = true, checks = false, check_length = false, remove_conserved = false,
-        remove_conserved_warn = true, structural_simplify = false, kwargs...)
+        structural_simplify = false, kwargs...)
     u0map = symmap_to_varmap(rs, u0)
     pmap = symmap_to_varmap(rs, p)
     sde_sys = convert(SDESystem, rs; name, combinatoric_ratelaws,
-        include_zero_odes, checks, remove_conserved, remove_conserved_warn)
+        include_zero_odes, checks, remove_conserved)
 
     # Handles potential differential algebraic equations (which requires `structural_simplify`).
     if structural_simplify
@@ -876,12 +850,12 @@ function DiffEqBase.SteadyStateProblem(rs::ReactionSystem, u0,
         p = DiffEqBase.NullParameters(), args...;
         check_length = false, name = nameof(rs),
         combinatoric_ratelaws = get_combinatoric_ratelaws(rs),
-        remove_conserved = false, remove_conserved_warn = true, include_zero_odes = true,
-        checks = false, structural_simplify = false, kwargs...)
+        remove_conserved = false, include_zero_odes = true, checks = false, 
+        structural_simplify = false, kwargs...)
     u0map = symmap_to_varmap(rs, u0)
     pmap = symmap_to_varmap(rs, p)
     osys = convert(ODESystem, rs; name, combinatoric_ratelaws, include_zero_odes, checks,
-        remove_conserved, remove_conserved_warn)
+        remove_conserved)
 
     # Handles potential differential algebraic equations (which requires `structural_simplify`).
     if structural_simplify

src/steady_state_stability.jl

@@ -118,7 +118,7 @@ function steady_state_jac(rs::ReactionSystem; u0 = [sp => 0.0 for sp in unknowns
     # Creates an `ODEProblem` with a Jacobian. Dummy values for `u0` and `ps` must be provided.
     ps = [p => 0.0 for p in parameters(rs)]
     return ODEProblem(rs, u0, 0, ps; jac = true, combinatoric_ratelaws,
-        remove_conserved = true, remove_conserved_warn = false)
+        remove_conserved = true)
 end
 
 # Converts a `u` vector from a vector of values to a map.