Merge branch 'master' into new_default_prob_update_tests

Author: TorkelE

docs/src/assets/long_ploting_times/model_creation/mm_kinetics.svg

@@ -191,11 +191,12 @@ u₀map = [:m₁ => 0, :m₂ => 0, :m₃ => 0, :P₁ => 20, :P₂ => 0, :P₃ =>
 dprob = DiscreteProblem(repressilator, u₀map, tspan, pmap)
 
 # now, we create a JumpProblem, and specify Gillespie's Direct Method as the solver:
-jprob = JumpProblem(repressilator, dprob, Direct(), save_positions=(false,false))
+jprob = JumpProblem(repressilator, dprob, Direct())
 
 # now, let's solve and plot the jump process:
-sol = solve(jprob, SSAStepper(), saveat=10.)
+sol = solve(jprob, SSAStepper())
 plot(sol)
+plot(sol, plotdensity = 1000, fmt = :png) # hide
 ```
 
 We see that oscillations remain, but become much noisier. Note, in constructing

docs/src/assets/long_ploting_times/model_creation/sir_outbreaks.svg

@@ -115,12 +115,13 @@ using Plots
 oplt = plot(osol; title = "Reaction rate equation (ODE)")
 splt = plot(ssol; title = "Chemical Langevin equation (SDE)")
 jplt = plot(jsol; title = "Stochastic chemical kinetics (Jump)")
-plot(oplt, splt, jplt; lw = 2, size=(800,800), layout = (3,1)) 
+plot(oplt, splt, jplt; lw = 2, size=(800,800), layout = (3,1))
+oplt = plot(osol; title = "Reaction rate equation (ODE)", plotdensity = 1000, fmt = :png) # hide
+splt = plot(ssol; title = "Chemical Langevin equation (SDE)", plotdensity = 1000, fmt = :png) # hide
+jplt = plot(jsol; title = "Stochastic chemical kinetics (Jump)", plotdensity = 1000, fmt = :png) # hide
+plot(oplt, splt, jplt; lw = 2, size=(800,800), layout = (3,1), plotdensity = 1000, fmt = :png) # hide
 plot!(bottom_margin = 3Plots.Measures.mm) # hide
-nothing # hide
 ```
-![MM Kinetics](../../assets/long_ploting_times/model_creation/mm_kinetics.svg)
-Note that, due to the large amounts of the species involved, the stochastic trajectories are very similar to the deterministic one.
 
 ## [SIR infection model](@id basic_CRN_library_sir)
 The [SIR model](https://en.wikipedia.org/wiki/Compartmental_models_in_epidemiology#The_SIR_model) is the simplest model of the spread of an infectious disease. While the real system is very different from the chemical and cellular processes typically modelled with CRNs, it (and several other epidemiological systems) can be modelled using the same CRN formalism. The SIR model consists of three species: susceptible ($S$), infected ($I$), and removed ($R$) individuals, and two reaction events: infection and recovery.
@@ -160,9 +161,11 @@ jplt1 = plot(jsol1; title = "Outbreak")
 jplt2 = plot(jsol2; title = "Outbreak")
 jplt3 = plot(jsol3; title = "No outbreak")
 plot(jplt1, jplt2, jplt3; lw = 3, size=(800,700), layout = (3,1))
-nothing # hide
+jplt1 = plot(jsol1; title = "Outbreak", plotdensity = 1000, fmt = :png) # hide
+jplt2 = plot(jsol2; title = "Outbreak", plotdensity = 1000, fmt = :png) # hide
+jplt3 = plot(jsol3; title = "No outbreak", plotdensity = 1000, fmt = :png) # hide
+plot(jplt1, jplt2, jplt3; lw = 3, size=(800,700), layout = (3,1), plotdensity = 1000, fmt = :png) # hide
 ```
-![SIR Outbreak](../../assets/long_ploting_times/model_creation/sir_outbreaks.svg)
 
 ## [Chemical cross-coupling](@id basic_CRN_library_cc)
 In chemistry, [cross-coupling](https://en.wikipedia.org/wiki/Cross-coupling_reaction) is when a catalyst combines two substrates to form a product. In this example, the catalyst ($C$) first binds one substrate ($S₁$) to form an intermediary complex ($S₁C$). Next, the complex binds the second substrate ($S₂$) to form another complex ($CP$). Finally, the catalyst releases the now-formed product ($P$). This system is an extended version of the [Michaelis-Menten system presented earlier](@ref basic_CRN_library_mm).

docs/src/assets/long_ploting_times/model_simulation/incomplete_brusselator_simulation.svg

@@ -31,9 +31,8 @@ oprob = ODEProblem(brusselator, u0, tspan, ps)
 
 sol1 = solve(oprob, Tsit5())
 plot(sol1)
-nothing # hide
+plot(sol1, plotdensity = 1000, fmt = :png) # hide
 ```
-![Incomplete Brusselator Simulation](../assets/long_ploting_times/model_simulation/incomplete_brusselator_simulation.svg)
 
 We get a warning, indicating that the simulation was terminated. Furthermore, the resulting plot ends at $t ≈ 12$, meaning that the simulation was not completed (as the simulation's endpoint is $t = 20$). Indeed, we can confirm this by checking the *return code* of the solution object:
 ```@example ode_simulation_performance_1

docs/src/introduction_to_catalyst/introduction_to_catalyst.md

@@ -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).

docs/src/model_creation/examples/basic_CRN_library.md

@@ -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);

docs/src/model_simulation/ode_simulation_performance.md

@@ -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.

ext/CatalystBifurcationKitExtension/bifurcation_kit_extension.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