rebase

Author: TorkelE

docs/pages.jl

@@ -37,7 +37,7 @@ pages = Any[
     "Steady state analysis" => Any[
         "steady_state_functionality/homotopy_continuation.md",
         "steady_state_functionality/nonlinear_solve.md",
-        # Stability analysis.
+        "steady_state_functionality/steady_state_stability_computation.md",        
         "steady_state_functionality/bifurcation_diagrams.md"
         # Dynamical systems analysis.
     ],

docs/src/catalyst_applications/steady_state_stability_computation.md renamed to docs/src/steady_state_functionality/steady_state_stability_computation.md

@@ -1,6 +1,9 @@
 # Steady state stability computation
 After system steady states have been found using [HomotopyContinuation.jl](@ref homotopy_continuation), [NonlinearSolve.jl](@ref nonlinear_solve), or other means, their stability can be computed using Catalyst's `steady_state_stability` function. Systems with conservation laws will automatically have these removed, permitting stability computation on systems with singular Jacobian.
 
+!!! warn 
+    Catalyst currently computes steady state stabilities using the naive approach of checking whether a system's largest eigenvalue is negative. While more advanced stability computation methods exist (and would be a welcome addition to Catalyst), there is no direct plans to implement these.
+
 ## Basic examples
 Let us consider the following basic example:
 ```@example stability_1
@@ -51,4 +54,10 @@ ps_2 = [:v => 4.0, :K => 1.5, :n => 2, :d => 1.0]
 steady_states_2 = hc_steady_states(sa_loop, ps)
 stability_2 = steady_state_stability(steady_states_2, sa_loop, ps_2; ss_jac=ss_jac)
 nothing # hide
+```
+
+It is possible to designate that a [sparse Jacobian](@ref ref) should be used using the `sparse = true` option (either to `steady_state_jac` or directly to `steady_state_stability`):
+```@example stability_1
+ss_jac = steady_state_jac(sa_loop; sparse = true)
+nothing # hide
 ```

ext/CatalystBifurcationKitExtension/bifurcation_kit_extension.jl

@@ -3,7 +3,9 @@
 # Creates a BifurcationProblem, using a ReactionSystem as an input.
 function BK.BifurcationProblem(rs::ReactionSystem, u0_bif, ps, bif_par, args...;
                                plot_var=nothing, record_from_solution=BK.record_sol_default, jac=true, u0=[], kwargs...)
-    is_autonomous(rs) || error("Attempting to compute steady state for non-autonomous system (e.g. does some rate depend on $(rs.iv)?), this is not possible.")
+    if !is_autonomous(rs) 
+        error("Attempting to create a `BifurcationProblem` for a non-autonomous system (e.g. where some rate depend on $(rs.iv)). This is not possible.")
+    end
 
     # Converts symbols to symbolics.
     (bif_par isa Symbol) && (bif_par = ModelingToolkit.get_var_to_name(rs)[bif_par])

ext/CatalystHomotopyContinuationExtension/homotopy_continuation_extension.jl

@@ -36,7 +36,9 @@ Notes:
 ```
   """
 function Catalyst.hc_steady_states(rs::ReactionSystem, ps; filter_negative=true, neg_thres=-1e-20, u0=[], kwargs...)
-    is_autonomous(rs) || error("Attempting to compute steady state for non-autonomous system (e.g. does some rate depend on $(rs.iv)?), this is not possible.")
+    if !is_autonomous(rs) 
+        error("Attempting to compute steady state for a non-autonomous system (e.g. where some rate depend on $(rs.iv)). This is not possible.")
+    end
     ss_poly = steady_state_polynomial(rs, ps, u0)
     sols = HC.real_solutions(HC.solve(ss_poly; kwargs...))
     reorder_sols!(sols, ss_poly, rs)

src/Catalyst.jl

@@ -103,6 +103,7 @@ export species, nonspecies, reactionparams, reactions, nonreactions, speciesmap,
 export numspecies, numreactions, numreactionparams, setdefaults!
 export make_empty_network, reactionparamsmap
 export dependants, dependents, substoichmat, prodstoichmat, netstoichmat
+export is_autonomous
 export reactionrates
 export isequivalent
 export set_default_noise_scaling
@@ -150,7 +151,7 @@ export @compound, @compounds
 export iscompound, components, coefficients, component_coefficients
 export balance_reaction
 
-# for functions I am unsure where to best place them.
+# Functionality for computing the stability of system steady states.
 include("steady_state_stability.jl")
 export steady_state_stability, steady_state_jac
 

src/reactionsystem.jl

@@ -1203,6 +1203,35 @@ function dependants(rx, network)
     dependents(rx, network)
 end
 
+"""
+    is_autonomous(rs::ReactionSystem)
+
+Checks if a system is autonomous (i.e. no rate or equation depend on the independent variable(s)).
+Example:
+```julia
+rs1 = @reaction_system
+    (p,d), 0 <--> X
+end
+is_autonomous(rs1) # Returns `true`.
+
+rs2 = @reaction_system
+    (p/t,d), 0 <--> X
+end
+is_autonomous(rs2) # Returns `false`.
+```
+"""
+function is_autonomous(rs::ReactionSystem)
+    # Get all variables occuring in reactions and then other equations.
+    dep_var_param_rxs = [ModelingToolkit.get_variables(rate) for rate in reactionrates(rs)]
+    dep_var_param_eqs = [ModelingToolkit.get_variables(eq) for eq in filter(eq -> !(eq isa Reaction), equations(rs))]
+    dep_var_param = reduce(vcat,[dep_var_param_rxs; dep_var_param_eqs])
+
+    # Checks for iv and spatial ivs
+    any(isequal(get_iv(rs), var) for var in dep_var_param) && (return false)
+    any(isequal(siv, var) for siv in get_sivs(rs) for var in dep_var_param) && (return false)
+    return true
+end
+
 
 ### `ReactionSystem` Remaking ###
 

src/reactionsystem_conversions.jl

@@ -518,8 +518,14 @@ function Base.convert(::Type{<:NonlinearSystem}, rs::ReactionSystem; name = name
                       include_zero_odes = true, remove_conserved = false, checks = false,
                       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)
+    if !is_autonomous(rs) 
+        error("Attempting to convert a non-autonomous `ReactionSystem` (e.g. where some rate depend on $(rs.iv)) to a `NonlinearSystem`. This is not possible. if you are intending to compute system steady states, consider creating and solving a `SteadyStateProblem.")
+    end
+
+    # Generates system equations.
     fullrs = Catalyst.flatten(rs)
     remove_conserved && conservationlaws(fullrs)
     ists, ispcs = get_indep_sts(fullrs, remove_conserved)

src/steady_state_stability.jl

@@ -1,18 +1,19 @@
 ### Stability Analysis ###
 
 """
-    stability(u::Vector{T}, rs::ReactionSystem, p; sparse=false, ss_jac = steady_state_jac(u, rs, p; sparse=sparse), t = Inf, non_autonomous_warn = true)
+    stability(u::Vector{T}, rs::ReactionSystem, p; sparse = false, 
+              ss_jac = steady_state_jac(u, rs, p; sparse=sparse))
 
 Compute the stability of a steady state (Returned as a `Bool`, with `true` indicating stability).
 
 Arguments:
     - `u`: The steady state for which we want to compute stability.
-    - `rs`: The reaction system model for which we want to compute stability.
+    - `rs`: The `ReactionSystem` model for which we want to compute stability.
     - `p`: The parameter set for which we want to compute stability.
-    - `sparse=false`: If we wish to create a sparse Jacobian for the stability computation.
-    - `ss_jac = steady_state_jac(u, rs; sparse=sparse)`: It is possible to pre-compute the Jacobian used for stability computation using `steady_state_jac`. If stability is computed for many states, precomputing the Jacobian may speed up evaluation.
-    - `t = Inf`: The time point at which stability is computed. 
-    - `non_autonomous_warn = true`: If the system is non-autonomous (e.g. a rate depends on t), a warning will be given. Set this to false to remove that. Alternatively, specify a nonInf value for `t`.
+    - `sparse = false`: If we wish to create a sparse Jacobian for the stability computation.
+    - `ss_jac = steady_state_jac(u, rs; sparse = sparse)`: It is possible to pre-compute the 
+      Jacobian used for stability computation using `steady_state_jac`. If stability is computed 
+      for many states, precomputing the Jacobian may speed up evaluation.
 
 Example:
 ```julia
@@ -22,45 +23,49 @@ rn = @reaction_network begin
 end
 p = [:p => 1.0, :d => 0.5]
 
-# Finds (the trivial) steady state, and computes stability.
+# Finds (the trivial) steady state, and computes its stability.
 steady_state = [2.0]
 steady_state_stability(steady_state, rn, p)
-```
 
 Notes:
-    y states (each being a `Vector`) is provided as `u`, stability is computed for each state separately.
+    - The input `u` can (currently) be both a vector of values (like `[1.0, 2.0]`) and a vector map 
+    (like `[X => 1.0, Y => 2.0]`). The reason is that most methods for computing stability only
+    produces vectors of values. However, if possible, it is recommended to work with values on the
+    form of maps.
+```
 """
-function steady_state_stability(u::Vector{T}, rs::ReactionSystem, p; 
-                                sparse=false, ss_jac = steady_state_jac(rs; u0=u, sparse=sparse), 
-                                t = Inf, non_autonomous_warn =true) where T
+function steady_state_stability(u::Vector, rs::ReactionSystem, ps; sparse = false, 
+                                ss_jac = steady_state_jac(rs; u0 = u, sparse = sparse))
     # Warning checks.
-    !is_autonomous(rs) && non_autonomous_warn && @warn "Attempting to compute stability for a non-autonomous system. Set `non_autonomous_warn = false` to disable this warning."
-
-    # Because Jacobian currently requires u and p to be a normal vector.
-    # Can be removed once this get fixed in MTK.
-    if (u isa Vector{<:Pair}) || (u isa Dict) 
-        u_dict = Dict(symmap_to_varmap(rs, u))
-        u = [u_dict[var] for var in states(rs)]        
+    if !is_autonomous(rs) 
+        error("Attempting to compute stability for a non-autonomous system (e.g. where some rate depend on $(rs.iv)). This is not possible.")
     end
-    if (p isa Vector{<:Pair}) || (p isa Dict)
-        p_dict = Dict(symmap_to_varmap(rs, p))
-        p = [p_dict[var] for var in parameters(rs)]
+
+    # If `u` is a vector of values, we convert it to a map. Also, if there are conservation laws,
+    # corresponding values must be eliminated from `u`.
+    u = steady_state_u_conversion(u, rs)
+    if length(u) > length(unknowns(ss_jac.f.sys))
+        u = filter(u_v -> any(isequal(u_v[1]), unknowns(ss_jac.f.sys)), u)
     end
 
     # Computes stability (by checking that the real part of all eigenvalues are negative).
-    jac = ss_jac(u, p, Inf)
-    return maximum(real.(eigvals(jac))) < 0
-end
-# Computes the stability for a vector of steady states.
-function steady_state_stability(us::Vector{Vector{T}}, rs::ReactionSystem, p; 
-                                sparse=false, ss_jac = steady_state_jac(rs; u0=us[1], sparse=sparse)) where T
-    return [steady_state_stability(u, rs, p) for u in us]
+    # Here, `ss_jac` is a `ODEProblem` with dummy values for `u0` and `p`.
+    J = zeros(length(u), length(u))
+    ss_jac = remake(ss_jac; u0 = u, p = ps)
+    ss_jac.f.jac(J, ss_jac.u0, ss_jac.p, Inf)
+    return maximum(real(ev) for ev in (eigvals(J))) < 0
 end
 
 """
-    steady_state_jac(rs::ReactionSystem; u0=[], sparse=false)
+    steady_state_jac(rs::ReactionSystem; u0 = [], sparse = false)
+
+Creates the Jacobian function which can be used as input to `steady_state_stability`. Useful when 
+a large number of stability computation has to be carried out in a performant manner.
 
-Creates the Jacobian function which can be used as input to `steady_state_stability`. Useful when a large number of stability computation has to be carried out in a performant manner.
+Arguments:
+    - `rs`: The reaction system model for which we want to compute stability.
+    - `u0 = []`: For systems with conservation laws, a `u` is required to compute the conserved quantities.
+    - `sparse = false`: If we wish to create a sparse Jacobian for the stability computation.
  
 Example:
 ```julia
@@ -76,28 +81,34 @@ ss_jac = steady_state_jacobian(rn)
 # Finds (the trivial) steady state, and computes stability.
 steady_state = [2.0]
 steady_state_stability(steady_state, rn, p; ss_jac)
-```
 
-Arguments:
-    - `rs`: The reaction system model for which we want to compute stability.
-    - `u0=[]`: For systems with conservation laws, a `u` is required to compute the conserved quantities.
-    - `sparse=false`: If we wish to create a sparse Jacobian for the stability computation.
+Notes:
+    - In practise, this function returns an `ODEProblem` (with a computed Jacobian) set up in 
+    such a way that it can be used by the `steady_state_stability` function.
+```
 """
-function steady_state_jac(rs::ReactionSystem; u0=[], sparse=false) where T
+function steady_state_jac(rs::ReactionSystem; u0 = [sp => 0.0 for sp in unknowns(rs)], 
+                          sparse = false, combinatoric_ratelaws = get_combinatoric_ratelaws(rs))
     # If u0 is a vector of values, must be converted to something MTK understands.
-    if (u0 isa Vector{<:Number})
-        if length(u0) == length(states(rs))
-            u0 = Pair.(states(rs), u0)
-        elseif !isempty(u0)
-            error("You are trying to compute a stability matrix, providing u0 to compute conservation laws. Your provided u0 vector has length < the number of system states. If you provide a u0 vector, these have to be identical.")
-        end
-    end
 
     # Converts u0 to values MTK understands, and checks that potential conservation laws are accounted for.
-    u0 = symmap_to_varmap(rs, u0)
+    u0 = steady_state_u_conversion(u0, rs)
     conservationlaw_errorcheck(rs, u0)
 
-    # Creates a ODESystem and jacobian.
-    osys = convert(ODESystem, rs; remove_conserved = true, defaults=Dict(u0))
-    return ODEFunction(osys; jac=true, sparse=sparse).jac
+    # 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 = get_combinatoric_ratelaws(rs))
+end
+
+# Converts a `u` vector from a vector of values to a map.
+function steady_state_u_conversion(u, rs::ReactionSystem)
+    if (u isa Vector{<:Number})
+        if length(u) == length(unknowns(rs))
+            u = [sp => v for (sp,v) in zip(unknowns(rs), u)]
+        else
+            error("You are trying to generate a stability Jacobian, providing u0 to compute conservation laws. Your provided u0 vector has length < the number of system states. If you provide a u0 vector, these have to be identical.")
+        end
+    end
+    return symmap_to_varmap(rs, u)
 end

test/miscellaneous_tests/api.jl

@@ -462,7 +462,6 @@ let
         (p + X*(p1+p2),d), 0 <--> X
         (kB,kD), 2X <--> X
     end
-
     @test !is_autonomous(rn1)
     @test !is_autonomous(rn2)
     @test is_autonomous(rn3)
@@ -488,9 +487,15 @@ let
         (p + X*(p1+p2),d), 0 <--> X
         (kB,kD), 2X <--> X
     end
-
     @test !is_autonomous(rn4)
     @test !is_autonomous(rn5)
     @test !is_autonomous(rn6)
     @test is_autonomous(rn7)
+
+    # Using a coupled CRN/equation model.
+    rn7 = @reaction_network begin
+        @equations D(V) ~ X/(1+t) - V
+        (p,d), 0 <--> X
+    end
+    @test !is_autonomous(rn7)
 end