From 9cc5f9cd63fb8994f2783b1ac05f9163928cb3f5 Mon Sep 17 00:00:00 2001 From: Torkel Date: Tue, 21 May 2024 13:29:38 -0400 Subject: [PATCH] doc changes --- .../steady_state_stability_computation.md | 18 +++++++++--------- src/steady_state_stability.jl | 4 ++-- 2 files changed, 11 insertions(+), 11 deletions(-) diff --git a/docs/src/steady_state_functionality/steady_state_stability_computation.md b/docs/src/steady_state_functionality/steady_state_stability_computation.md index 7f84c78899..fdd12530e8 100644 --- a/docs/src/steady_state_functionality/steady_state_stability_computation.md +++ b/docs/src/steady_state_functionality/steady_state_stability_computation.md @@ -2,7 +2,7 @@ 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 real part 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. Furthermore, Catalyst uses a arbitrary tolerance $tol ~ 1.5*10^-7$ to determine whether a computed eigenvalue is far away enough from 0 to be reliably used. This selected threshold, however, have not been subject to further analysis (and can be changed through the `tol` argument). + Catalyst currently computes steady state stabilities using the naive approach of checking whether a system's largest eigenvalue real part 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. Furthermore, Catalyst uses a arbitrary tolerance $tol ≈ 1.5*10^-7$ to determine whether a computed eigenvalue is far away enough from 0 to be reliably used. This selected threshold, however, have not been subject to further analysis (and can be changed through the `tol` argument). ## Basic examples Let us consider the following basic example: @@ -19,7 +19,7 @@ steady_state = [:X => 4.0] steady_state_stability(steady_state, rn, ps) ``` -Next, let us consider the following self-activation loop: +Next, let us consider the following [self-activation loop](@ref ref): ```@example stability_1 sa_loop = @reaction_network begin (hill(X,v,K,n),d), 0 <--> X @@ -31,18 +31,15 @@ import HomotopyContinuation ps = [:v => 2.0, :K => 0.5, :n => 3, :d => 1.0] steady_states = hc_steady_states(sa_loop, ps) ``` -Next, we can apply `steady_state_stability` directly to this steady state vector, receiving the stability for each: +Next, we can apply `steady_state_stability` to each steady state yielding a vector of their stabilities: ```@example stability_1 -steady_state_stability(steady_states, sa_loop, ps) +[steady_state_stability(sstate, sa_loop, ps) for sstate in steady_states] ``` -!!! note - For systems with [conservation laws](@ref homotopy_continuation_conservation_laws), `steady_state_jac` must be supplied a `u0` vector (indicating species concentrations for conservation law computation). This is required to eliminate the conserved quantities, preventing a singular Jacobian. These are supplied using the `u0` optional argument. - ## Pre-computing the Jacobian to increase performance when computing stability for many steady states Catalyst uses the system Jacobian to compute steady state stability, and the Jacobian is computed once for each call to `steady_state_stability`. If you repeatedly compute stability for steady states of the same system, pre-computing the Jacobian and supplying it to the `steady_state_stability` function can improve performance. -In this example we use the self-activation loop from previously, pre-computes the Jacobian, and uses it to multiple `steady_state_stability` calls: +In this example we use the self-activation loop from previously, pre-computes its Jacobian, and uses it to multiple `steady_state_stability` calls: ```@example stability_1 ss_jac = steady_state_jac(sa_loop) @@ -60,4 +57,7 @@ It is possible to designate that a [sparse Jacobian](@ref ref) should be used us ```@example stability_1 ss_jac = steady_state_jac(sa_loop; sparse = true) nothing # hide -``` \ No newline at end of file +``` + +!!! warn + For systems with [conservation laws](@ref homotopy_continuation_conservation_laws), `steady_state_jac` must be supplied a `u0` vector (indicating species concentrations for conservation law computation). This is required to eliminate the conserved quantities, preventing a singular Jacobian. These are supplied using the `u0` optional argument. \ No newline at end of file diff --git a/src/steady_state_stability.jl b/src/steady_state_stability.jl index e4b5c4a5cb..dc901f577b 100644 --- a/src/steady_state_stability.jl +++ b/src/steady_state_stability.jl @@ -119,8 +119,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 = get_combinatoric_ratelaws(rs)) + return ODEProblem(rs, u0, 0, ps; jac = true, remove_conserved = true, sparse = sparse, + combinatoric_ratelaws = combinatoric_ratelaws) end # Converts a `u` vector from a vector of values to a map.