Merge pull request #955 from SciML/conservation_laws_warning

[v14] Conservation law warnings

Author: TorkelE

ext/CatalystBifurcationKitExtension/bifurcation_kit_extension.jl

@@ -22,7 +22,8 @@ function BK.BifurcationProblem(rs::ReactionSystem, u0_bif, ps, bif_par, args...;
 
     # Creates NonlinearSystem.
     Catalyst.conservationlaw_errorcheck(rs, vcat(ps, u0))
-    nsys = convert(NonlinearSystem, rs; remove_conserved = true, defaults = Dict(u0))
+    nsys = convert(NonlinearSystem, rs; defaults = Dict(u0),
+        remove_conserved = true, remove_conserved_warn = false)
     nsys = complete(nsys)
 
     # Makes BifurcationProblem (this call goes through the ModelingToolkit-based BifurcationKit extension).

ext/CatalystHomotopyContinuationExtension/homotopy_continuation_extension.jl

@@ -49,7 +49,8 @@ end
 # For a given reaction system, parameter values, and initial conditions, find the polynomial that HC solves to find steady states.
 function steady_state_polynomial(rs::ReactionSystem, ps, u0)
     rs = Catalyst.expand_registered_functions(rs)
-    ns = complete(convert(NonlinearSystem, rs; remove_conserved = true))
+    ns = complete(convert(NonlinearSystem, rs;
+        remove_conserved = true, remove_conserved_warn = false))
     pre_varmap = [symmap_to_varmap(rs, u0)..., symmap_to_varmap(rs, ps)...]
     Catalyst.conservationlaw_errorcheck(rs, pre_varmap)
     p_vals = ModelingToolkit.varmap_to_vars(pre_varmap, parameters(ns);

ext/CatalystStructuralIdentifiabilityExtension/structural_identifiability_extension.jl

@@ -167,7 +167,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))
+    osys = complete(convert(ODESystem, rs; remove_conserved, remove_conserved_warn = false))
     vars = [unknowns(rs); parameters(rs)]
 
     # Computes equations for system conservation laws.

src/reactionsystem_conversions.jl

@@ -450,6 +450,22 @@ 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 removeal 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 ###
 
 """
@@ -467,15 +483,20 @@ 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, checks = false,
-        default_u0 = Dict(), default_p = Dict(),
+        include_zero_odes = true, remove_conserved = false, remove_conserved_warn = true,
+        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)
     ists, ispcs = get_indep_sts(fullrs, remove_conserved)
@@ -509,16 +530,19 @@ 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,
-        default_u0 = Dict(), default_p = Dict(),
+        remove_conserved_warn = true, 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
@@ -589,17 +613,19 @@ 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.
-- Does not currently support `ReactionSystem`s that include coupled algebraic or
-  differential 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,
-        default_u0 = Dict(), default_p = Dict(),
+        remove_conserved_warn = true, 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)
 
@@ -651,7 +677,6 @@ function Base.convert(::Type{<:JumpSystem}, rs::ReactionSystem; name = nameof(rs
         kwargs...)
     iscomplete(rs) || error(COMPLETENESS_ERROR)
     spatial_convert_err(rs::ReactionSystem, JumpSystem)
-
     (remove_conserved !== nothing) &&
         throw(ArgumentError("Catalyst does not support removing conserved species when converting to JumpSystems."))
 
@@ -684,12 +709,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,
+        include_zero_odes = true, remove_conserved = false, remove_conserved_warn = 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, remove_conserved_warn)
 
     # Handles potential differential algebraic equations (which requires `structural_simplify`).
     if structural_simplify
@@ -708,12 +733,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, checks = false,
+        remove_conserved = false, remove_conserved_warn = true, 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)
+        checks, all_differentials_permitted, remove_conserved, remove_conserved_warn)
     nlsys = complete(nlsys)
     return NonlinearProblem(nlsys, u0map, pmap, args...; check_length,
         kwargs...)
@@ -723,12 +748,12 @@ end
 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, structural_simplify = false, kwargs...)
+        include_zero_odes = true, checks = false, check_length = false, remove_conserved = false,
+        remove_conserved_warn = true, 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)
+        include_zero_odes, checks, remove_conserved, remove_conserved_warn)
 
     # Handles potential differential algebraic equations (which requires `structural_simplify`).
     if structural_simplify
@@ -772,12 +797,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, include_zero_odes = true,
+        remove_conserved = false, remove_conserved_warn = true, 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, remove_conserved_warn)
 
     # Handles potential differential algebraic equations (which requires `structural_simplify`).
     if structural_simplify

src/steady_state_stability.jl

@@ -117,8 +117,8 @@ 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, remove_conserved = true,
-        combinatoric_ratelaws = combinatoric_ratelaws)
+    return ODEProblem(rs, u0, 0, ps; jac = true, combinatoric_ratelaws,
+        remove_conserved = true, remove_conserved_warn = false)
 end
 
 # Converts a `u` vector from a vector of values to a map.

test/network_analysis/conservation_laws.jl

@@ -11,6 +11,9 @@ seed = rand(rng, 1:100)
 # Fetch test networks.
 include("../test_networks.jl")
 
+# Except where we test the warnings, we do not want to print this warning.
+remove_conserved_warn = false
+
 ### Basic Tests ###
 
 # Tests basic functionality on system with known conservation laws.
@@ -114,23 +117,23 @@ let
     tspan = (0.0, 20.0)
 
     # Simulates model using ODEs and checks that simulations are identical.
-    osys = complete(convert(ODESystem, rn; remove_conserved = true))
+    osys = complete(convert(ODESystem, rn; remove_conserved = true, remove_conserved_warn))
     oprob1 = ODEProblem(osys, u0, tspan, p)
     oprob2 = ODEProblem(rn, u0, tspan, p)
-    oprob3 = ODEProblem(rn, u0, tspan, p; remove_conserved = true)
+    oprob3 = ODEProblem(rn, u0, tspan, p; remove_conserved = true, remove_conserved_warn)
     osol1 = solve(oprob1, Tsit5(); abstol = 1e-8, reltol = 1e-8, saveat= 0.2)
     osol2 = solve(oprob2, Tsit5(); abstol = 1e-8, reltol = 1e-8, saveat= 0.2)
     osol3 = solve(oprob3, Tsit5(); abstol = 1e-8, reltol = 1e-8, saveat= 0.2)
     @test osol1[sps] ≈ osol2[sps] ≈ osol3[sps]
 
     # Checks that steady states found using nonlinear solving and steady state simulations are identical.
-    nsys = complete(convert(NonlinearSystem, rn; remove_conserved = true))
+    nsys = complete(convert(NonlinearSystem, rn; remove_conserved = true, remove_conserved_warn))
     nprob1 = NonlinearProblem{true}(nsys, u0, p)
     nprob2 = NonlinearProblem(rn, u0, p)
-    nprob3 = NonlinearProblem(rn, u0, p; remove_conserved = true)
+    nprob3 = NonlinearProblem(rn, u0, p; remove_conserved = true, remove_conserved_warn)
     ssprob1 = SteadyStateProblem{true}(osys, u0, p)
     ssprob2 = SteadyStateProblem(rn, u0, p)
-    ssprob3 = SteadyStateProblem(rn, u0, p; remove_conserved = true)
+    ssprob3 = SteadyStateProblem(rn, u0, p; remove_conserved = true, remove_conserved_warn)
     nsol1 = solve(nprob1, NewtonRaphson(); abstol = 1e-8)
     # Nonlinear problems cannot find steady states properly without removing conserved species.
     nsol3 = solve(nprob3, NewtonRaphson(); abstol = 1e-8)
@@ -142,10 +145,10 @@ let
     # Creates SDEProblems using various approaches.
     u0_sde = [A => 100.0, B => 20.0, C => 5.0, D => 10.0, E => 3.0, F1 => 8.0, F2 => 2.0,
         F3 => 20.0]
-    ssys = complete(convert(SDESystem, rn; remove_conserved = true))
+    ssys = complete(convert(SDESystem, rn; remove_conserved = true, remove_conserved_warn))
     sprob1 = SDEProblem(ssys, u0_sde, tspan, p)
     sprob2 = SDEProblem(rn, u0_sde, tspan, p)
-    sprob3 = SDEProblem(rn, u0_sde, tspan, p; remove_conserved = true)
+    sprob3 = SDEProblem(rn, u0_sde, tspan, p; remove_conserved = true, remove_conserved_warn)
 
     # Checks that the SDEs f and g function evaluates to the same thing.
     ind_us = ModelingToolkit.get_unknowns(ssys)
@@ -186,9 +189,9 @@ let
     u0_2 = [rn.X => 2.0, rn.Y => 3.0, rn.XY => 4.0]
     u0_3 = [:X => 2.0, :Y => 3.0, :XY => 4.0]
     ps = (kB => 2, kD => 1.5)
-    oprob1 = ODEProblem(rn, u0_1, 10.0, ps; remove_conserved = true)
-    oprob2 = ODEProblem(rn, u0_2, 10.0, ps; remove_conserved = true)
-    oprob3 = ODEProblem(rn, u0_3, 10.0, ps; remove_conserved = true)
+    oprob1 = ODEProblem(rn, u0_1, 10.0, ps; remove_conserved = true, remove_conserved_warn)
+    oprob2 = ODEProblem(rn, u0_2, 10.0, ps; remove_conserved = true, remove_conserved_warn)
+    oprob3 = ODEProblem(rn, u0_3, 10.0, ps; remove_conserved = true, remove_conserved_warn)
     @test solve(oprob1)[sps] ≈ solve(oprob2)[sps] ≈ solve(oprob3)[sps]
 end
 
@@ -200,7 +203,7 @@ let
     end
     u0 = Dict([:X1 => 100.0, :X2 => 120.0])
     ps = [:k1 => 0.2, :k2 => 0.15]
-    sprob = SDEProblem(rn, u0, 10.0, ps; remove_conserved = true)
+    sprob = SDEProblem(rn, u0, 10.0, ps; remove_conserved = true, remove_conserved_warn)
 
     # Checks that conservation laws hold in all simulations.
     sol = solve(sprob, ImplicitEM(); seed)
@@ -213,7 +216,7 @@ let
     rn = @reaction_network begin
         (k1,k2), X1 <--> X2
     end
-    osys = complete(convert(ODESystem, rn; remove_conserved = true))
+    osys = complete(convert(ODESystem, rn; remove_conserved = true, remove_conserved_warn))
     u0 = [osys.X1 => 1.0, osys.X2 => 1.0]
     ps_1 = [osys.k1 => 2.0, osys.k2 => 3.0]
     ps_2 = [osys.k1 => 2.0, osys.k2 => 3.0, osys.Γ[1] => 4.0]
@@ -234,7 +237,7 @@ let
     rn = @reaction_network begin
         (k1,k2), X1 <--> X2
     end
-    @test_throws ArgumentError convert(JumpSystem, rn; remove_conserved = true)
+    @test_throws ArgumentError convert(JumpSystem, rn; remove_conserved = true, remove_conserved_warn)
 end
 
 # Checks that `conserved` metadata is added correctly to parameters.
@@ -245,11 +248,44 @@ let
         (k1,k2), X1 <--> X2
         (k1,k2), Y1 <--> Y2
     end
-    osys = convert(ODESystem, rs; remove_conserved = true)
+    osys = convert(ODESystem, rs; remove_conserved = true, remove_conserved_warn)
 
     # Checks that the correct parameters have the `conserved` metadata.
     @test Catalyst.isconserved(osys.Γ[1])
     @test Catalyst.isconserved(osys.Γ[2])
     @test !Catalyst.isconserved(osys.k1)
     @test !Catalyst.isconserved(osys.k2)
 end
+
+# Checks that conservation law elimination warnings are generated in the correct cases.
+let
+    # Prepare model.
+    rn = @reaction_network begin 
+        (k1,k2), X1 <--> X2
+    end
+    u0 = [:X1 => 1.0, :X2 => 2.0]
+    tspan = (0.0, 1.0)
+    ps = [:k1 => 3.0, :k2 => 4.0]
+
+    # Check warnings in system conversion.
+    for XSystem in [ODESystem, SDESystem, NonlinearSystem]
+        @test_nowarn convert(XSystem, rn)
+        @test_logs (:warn, r"You are creating a system or problem while eliminating conserved quantities. Please *") convert(XSystem, rn; remove_conserved = true)
+        @test_nowarn convert(XSystem, rn; remove_conserved_warn = false)
+        @test_nowarn convert(XSystem, rn; remove_conserved = true, remove_conserved_warn = false)
+    end
+
+    # Checks during problem creation (separate depending on whether they have a time span or not).
+    for XProblem in [ODEProblem, SDEProblem]
+        @test_nowarn XProblem(rn, u0, tspan, ps)
+        @test_logs (:warn, r"You are creating a system or problem while eliminating conserved quantities. Please *") XProblem(rn, u0, tspan, ps; remove_conserved = true)
+        @test_nowarn XProblem(rn, u0, tspan, ps; remove_conserved_warn = false)
+        @test_nowarn XProblem(rn, u0, tspan, ps; remove_conserved = true, remove_conserved_warn = false)
+    end
+    for XProblem in [NonlinearProblem, SteadyStateProblem]
+        @test_nowarn XProblem(rn, u0, ps)
+        @test_logs (:warn, r"You are creating a system or problem while eliminating conserved quantities. Please *") XProblem(rn, u0, ps; remove_conserved = true)
+        @test_nowarn XProblem(rn, u0, ps; remove_conserved_warn = false)
+        @test_nowarn XProblem(rn, u0, ps; remove_conserved = true, remove_conserved_warn = false)
+    end
+end