From 3d4afd71d807366fa5e3acfe0f5e28e2ee039b63 Mon Sep 17 00:00:00 2001 From: Torkel Date: Mon, 4 Dec 2023 10:56:21 -0500 Subject: [PATCH 01/28] init --- docs/src/api.md | 1 + .../steady_state_stability_computation.md | 54 ++++++++++ src/Catalyst.jl | 5 + src/steady_state_stability.jl | 100 ++++++++++++++++++ test/miscellaneous_tests/api.jl | 48 +++++++++ .../stability_computation.jl | 79 ++++++++++++++ 6 files changed, 287 insertions(+) create mode 100644 docs/src/catalyst_functionality/steady_state_stability_computation.md create mode 100644 src/steady_state_stability.jl create mode 100644 test/miscellaneous_tests/stability_computation.jl diff --git a/docs/src/api.md b/docs/src/api.md index 51a403c1a5..5e267791a3 100644 --- a/docs/src/api.md +++ b/docs/src/api.md @@ -166,6 +166,7 @@ speciesmap paramsmap reactionparamsmap isspecies +is_autonomous Catalyst.isconstant Catalyst.isbc ``` diff --git a/docs/src/catalyst_functionality/steady_state_stability_computation.md b/docs/src/catalyst_functionality/steady_state_stability_computation.md new file mode 100644 index 0000000000..4f812ce27d --- /dev/null +++ b/docs/src/catalyst_functionality/steady_state_stability_computation.md @@ -0,0 +1,54 @@ +# 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. + +## Basic examples +Let us consider the following basic example: +```@example stability_1 +using Catalyst +rn = @reaction_network begin + (p,d), 0 <--> X +end +``` +It has a single (stable) steady state at $X = p/d$. We can confirm stability using the `steady_state_stability` function, to which we provide the steady state, the reaction system, and the parameter values: +```@example stability_1 +ps = [:p => 2.0, :d => 0.5] +steady_state = [:X => 4.0] +steady_state_stability(steady_state, rn, ps) +``` + +Next, let us consider the following self-activation loop: +```@example stability_1 +sa_loop = @reaction_network begin + (hill(X,v,K,n),d), 0 <--> X +end +``` +For certain parameter choices, this system exhibits multi-stability. Here, we can find the steady states [using homotopy continuation](@ref homotopy_continuation): +```@example stability_1 +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: +```@example stability_1 +steady_state_stability(steady_states, sa_loop, ps) +``` + +## 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: +```@example stability_1 +ss_jac = steady_state_jac(sa_loop) + +ps_1 = [:v => 2.0, :K => 0.5, :n => 3, :d => 1.0] +steady_states_1 = hc_steady_states(sa_loop, ps) +stability_1 = steady_state_stability(steady_states_1, sa_loop, ps_1; ss_jac=ss_jac) + +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 +``` + +!!! 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. \ No newline at end of file diff --git a/src/Catalyst.jl b/src/Catalyst.jl index b347f0ab31..f87548b47f 100644 --- a/src/Catalyst.jl +++ b/src/Catalyst.jl @@ -42,6 +42,7 @@ import Base: (==), hash, size, getindex, setindex, isless, Sort.defalg, length, import MacroTools, Graphs import Graphs: DiGraph, SimpleGraph, SimpleDiGraph, vertices, edges, add_vertices!, nv, ne import DataStructures: OrderedDict, OrderedSet +import LinearAlgebra.eigvals import Parameters: @with_kw_noshow import Symbolics: occursin, wrap @@ -149,6 +150,10 @@ export @compound, @compounds export iscompound, components, coefficients, component_coefficients export balance_reaction +# for functions I am unsure where to best place them. +include("steady_state_stability.jl") +export steady_state_stability, steady_state_jac + ### Extensions ### # HomotopyContinuation diff --git a/src/steady_state_stability.jl b/src/steady_state_stability.jl new file mode 100644 index 0000000000..33cf55529b --- /dev/null +++ b/src/steady_state_stability.jl @@ -0,0 +1,100 @@ +### 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_war=true) + +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. + - `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_war=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`. + +Example: +```julia +# Creates model. +rn = @reaction_network begin + (p,d), 0 <--> X +end +p = [:p => 1.0, :d => 0.5] + +# Finds (the trivial) steady state, and computes 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. +""" +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_war=true) where T + # Warning checks. + !is_autonomous(rs) && non_autonomous_war && @warn "Attempting to compute stability for a non-autonomous system. Set `non_autonomous_war=false` to disable this warning." + + # Because Jacobian currently requires ps 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)] + 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)] + end + + # Computes stability + 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] +end + +""" + 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. + +Example: +```julia +# Creates model. +rn = @reaction_network begin + (p,d), 0 <--> X +end +p = [:p => 1.0, :d => 0.5] + +# Creates the steady state jacobian. +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. +""" +function steady_state_jac(rs::ReactionSystem; u0=[], sparse=false) where T + # 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) + 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 +end \ No newline at end of file diff --git a/test/miscellaneous_tests/api.jl b/test/miscellaneous_tests/api.jl index 2c42475628..1340f5d85e 100644 --- a/test/miscellaneous_tests/api.jl +++ b/test/miscellaneous_tests/api.jl @@ -446,3 +446,51 @@ let neweqs = getfield.(equations(ns), :rhs) @test_throws MethodError Catalyst.to_multivariate_poly(neweqs) end + +# Tests `is_autonomous` function. +let + # Using default iv. + rn1 = @reaction_network begin + (p + X*(p1/(t+p3)),d), 0 <--> X + (kB,kD), 2X <--> X + end + rn2 = @reaction_network begin + (hill(X, v/t, K, n),d), 0 <--> X + (kB,kD), 2X <--> X + end + rn3 = @reaction_network begin + (p + X*(p1+p2),d), 0 <--> X + (kB,kD), 2X <--> X + end + + @test !is_autonomous(rn1) + @test !is_autonomous(rn2) + @test is_autonomous(rn3) + + # Using non-default iv. + rn4 = @reaction_network begin + @ivs i1 i2 + (p + X*(p1/(1+i1)),d), 0 <--> X + (kB,kD), 2X <--> X + end + rn5 = @reaction_network begin + @ivs i1 i2 + (p + X*(i2+p2),d), 0 <--> X + (kB,kD), 2X <--> X + end + rn6 = @reaction_network begin + @ivs i1 i2 + (hill(X, v/i1, i2, n),d), 0 <--> X + (kB,kD), 2X <--> X + end + rn7 = @reaction_network begin + @ivs i1 i2 + (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) +end \ No newline at end of file diff --git a/test/miscellaneous_tests/stability_computation.jl b/test/miscellaneous_tests/stability_computation.jl new file mode 100644 index 0000000000..93c22784d8 --- /dev/null +++ b/test/miscellaneous_tests/stability_computation.jl @@ -0,0 +1,79 @@ +### Fetch Packages ### + +# Fetch packages. +using Catalyst, OrdinaryDiffEq +using Random, Test + +# Sets rnd number. +using StableRNGs +rng = StableRNG(12345) + +### Run Tests ### + +# Tests that stability is correctly assessed (using simulation) in multi stable system. +# Tests that `steady_state_jac` function works. +# Tests with and without sparsity. +# tests using symbolic input. +let + # System which may have between 1 and 7 fixed points. + rn = @reaction_network begin + v/20.0 + hillar(X,Y,v,K,n), 0 --> X + v/20.0 + hillar(Y,X,v,K,n), 0 --> Y + d, (X,Y) --> 0 + end + ss_jac = steady_state_jac(rn) + ss_jac_sparse = steady_state_jac(rn; sparse=true) + + # Repeats several times, most cases should be encountered several times. + for i = 1:50 + # Generates random parameter values (which can generate all steady states cases). + p = [:v => 1.0 + 3*rand(rng), :K => 0.5 + 2*rand(rng), :n => rand(rng,[1,2,3,4]), :d => 0.5 + rand(rng)] + + sss = hc_steady_states(rn, p) + stabs_1 = steady_state_stability(sss, rn, p) + stabs_2 = steady_state_stability(sss, rn, p; sparse=true) + stabs_3 = steady_state_stability(sss, rn, p; ss_jac=ss_jac) + stabs_4 = steady_state_stability(sss, rn, p; ss_jac=ss_jac_sparse) + + for (idx,ss) in enumerate(sss) + oprob = ODEProblem(rn, [1.001, 0.999] .* ss, (0.0,1000.0), p) + sol_end = solve(oprob, Rosenbrock23())[end] + stabs_5 = ss ≈ sol_end + @test stabs_1[idx] == stabs_2[idx] == stabs_3[idx] == stabs_4[idx] == stabs_5 + end + end +end + +# Checks stability for system with known stability structure. +# Tests for system with conservation laws. +# Tests for various input forms of u0 and ps. +let + rn = complete(@reaction_network begin + k1+Z, Y --> 2X + k2, 2X --> X + Y + k3, X + Y --> Y + k4, X --> 0 + (kD1+X, kD2), 2Z <--> Z2 + end) + @unpack k1, k2, k3, k4, kD1, kD2, X, Y, Z, Z2 = rn + u0_1 = [X => 1.0, Y => 1.0, Z => 1.0, Z2 => 1.0] + u0_2 = [:X => 1.0, :Y => 1.0, :Z => 1.0, :Z2 => 1.0] + u0_3 = [rn.X => 1.0, rn.Y => 1.0, rn.Z => 1.0, rn.Z2 => 1.0] + u0_4 = [1.0, 1.0, 1.0, 1.0] + ps_1 = [k1 => 8.0, k2 => 2.0, k3 => 1.0, k4 => 1.5, kD1 => 0.5, kD2 => 2.0] + ps_2 = [:k1 => 8.0, :k2 => 2.0, :k3 => 1.0, :k4 => 1.5, :kD1 => 0.5, :kD2 => 2.0] + ps_3 = [rn.k1 => 8.0, rn.k2 => 2.0, rn.k3 => 1.0, rn.k4 => 1.5, rn.kD1 => 0.5, rn.kD2 => 4.0] + ps_4 = [8.0, 2.0, 1.0, 1.5, 0.5, 4.0] + + sss = hc_steady_states(rn, ps_1; u0=u0_1) + for u0 in [u0_1, u0_2, u0_3, u0_4], ps in [ps_1, ps_2, ps_3, ps_4] + stab_1 = steady_state_stability(sss, rn, ps) + @test length(stab_1) == 3 + @test count(stab_1) == 2 + + ss_jac = steady_state_jac(rn; u0=u0) + stab_2 = steady_state_stability(sss, rn, ps; ss_jac=ss_jac) + @test length(stab_2) == 3 + @test count(stab_2) == 2 + end +end \ No newline at end of file From f2c460c05b54cbc1aa163293116f58a78b3639f7 Mon Sep 17 00:00:00 2001 From: Torkel Date: Mon, 4 Dec 2023 11:02:18 -0500 Subject: [PATCH 02/28] add non-autonomous errors to bifurcation computations and homotopy continuation computation. --- .../bifurcation_kit_extension.jl | 1 + .../homotopy_continuation_extension.jl | 1 + test/extensions/homotopy_continuation.jl | 9 +++++++++ 3 files changed, 11 insertions(+) diff --git a/ext/CatalystBifurcationKitExtension/bifurcation_kit_extension.jl b/ext/CatalystBifurcationKitExtension/bifurcation_kit_extension.jl index f99389c6df..7c940190df 100644 --- a/ext/CatalystBifurcationKitExtension/bifurcation_kit_extension.jl +++ b/ext/CatalystBifurcationKitExtension/bifurcation_kit_extension.jl @@ -3,6 +3,7 @@ # 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.") # Converts symbols to symbolics. (bif_par isa Symbol) && (bif_par = ModelingToolkit.get_var_to_name(rs)[bif_par]) diff --git a/ext/CatalystHomotopyContinuationExtension/homotopy_continuation_extension.jl b/ext/CatalystHomotopyContinuationExtension/homotopy_continuation_extension.jl index 6d0a762c27..ca4e6408c4 100644 --- a/ext/CatalystHomotopyContinuationExtension/homotopy_continuation_extension.jl +++ b/ext/CatalystHomotopyContinuationExtension/homotopy_continuation_extension.jl @@ -36,6 +36,7 @@ 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.") ss_poly = steady_state_polynomial(rs, ps, u0) sols = HC.real_solutions(HC.solve(ss_poly; kwargs...)) reorder_sols!(sols, ss_poly, rs) diff --git a/test/extensions/homotopy_continuation.jl b/test/extensions/homotopy_continuation.jl index bc0f30e3a4..d669be5460 100644 --- a/test/extensions/homotopy_continuation.jl +++ b/test/extensions/homotopy_continuation.jl @@ -138,4 +138,13 @@ let # Computes bifurcation diagram. @test_throws Exception hc_steady_states(incomplete_network, p_start) +end + +# Tests that non-autonomous system throws an error +let + rs = @reaction_network begin + (k,t), 0 <--> X + end + ps = [k => 1.0] + @test_throws Exception hc_steady_states(rs, ps) end \ No newline at end of file From f35157a95b92ab0290bba74a0fac0b938e6f447d Mon Sep 17 00:00:00 2001 From: Torkel Date: Mon, 4 Dec 2023 11:03:57 -0500 Subject: [PATCH 03/28] update history file --- HISTORY.md | 13 +++++++++++++ 1 file changed, 13 insertions(+) diff --git a/HISTORY.md b/HISTORY.md index 3b51630da3..41b23916d5 100644 --- a/HISTORY.md +++ b/HISTORY.md @@ -85,6 +85,19 @@ plot(bif_dia; xguide="k1", yguide="X") ``` - Automatically handles elimination of conservation laws for computing bifurcation diagrams. - Updated Bifurcation documentation with respect to this new feature. +- Added function `is_autonomous` to check if a `ReactionSystem` is autonomous. +- Added function `steady_state_stability` to compute stability for steady states. Example: +```julia +# Creates model. +rn = @reaction_network begin + (p,d), 0 <--> X +end +p = [:p => 1.0, :d => 0.5] + +# Finds (the trivial) steady state, and computes stability. +steady_state = [2.0] +steady_state_stability(steady_state, rn, p) +``` ## Catalyst 13.5 - Added a CatalystHomotopyContinuationExtension extension, which exports the `hc_steady_state` function if HomotopyContinuation is exported. `hc_steady_state` finds the steady states of a reaction system using the homotopy continuation method. This feature is only available for julia versions 1.9+. Example: From ab02348c7e33efb542c7bc4d2d651ec47d2ad3ac Mon Sep 17 00:00:00 2001 From: Torkel Date: Mon, 4 Dec 2023 11:05:28 -0500 Subject: [PATCH 04/28] Code readability improvements. --- src/steady_state_stability.jl | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/src/steady_state_stability.jl b/src/steady_state_stability.jl index 33cf55529b..dd197091e9 100644 --- a/src/steady_state_stability.jl +++ b/src/steady_state_stability.jl @@ -35,7 +35,8 @@ function steady_state_stability(u::Vector{T}, rs::ReactionSystem, p; # Warning checks. !is_autonomous(rs) && non_autonomous_war && @warn "Attempting to compute stability for a non-autonomous system. Set `non_autonomous_war=false` to disable this warning." - # Because Jacobian currently requires ps to be a normal vector, can be removed once this get fixed in MTK. + # 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)] @@ -45,12 +46,13 @@ function steady_state_stability(u::Vector{T}, rs::ReactionSystem, p; p = [p_dict[var] for var in parameters(rs)] end - # Computes stability + # 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 +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] end From 0aa46a0b9d99ba8e2243d07b2c3a9f09b7b63ac5 Mon Sep 17 00:00:00 2001 From: Torkel Date: Mon, 4 Dec 2023 12:06:40 -0500 Subject: [PATCH 05/28] up --- test/miscellaneous_tests/stability_computation.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/test/miscellaneous_tests/stability_computation.jl b/test/miscellaneous_tests/stability_computation.jl index 93c22784d8..ac8054e762 100644 --- a/test/miscellaneous_tests/stability_computation.jl +++ b/test/miscellaneous_tests/stability_computation.jl @@ -3,6 +3,7 @@ # Fetch packages. using Catalyst, OrdinaryDiffEq using Random, Test +import HomotopyContinuation # Sets rnd number. using StableRNGs From 502f264f96f62d26463d0bf47befd72421a8f85c Mon Sep 17 00:00:00 2001 From: Torkel Date: Mon, 4 Dec 2023 12:33:45 -0500 Subject: [PATCH 06/28] up --- test/extensions/homotopy_continuation.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/extensions/homotopy_continuation.jl b/test/extensions/homotopy_continuation.jl index d669be5460..7a260e3b59 100644 --- a/test/extensions/homotopy_continuation.jl +++ b/test/extensions/homotopy_continuation.jl @@ -145,6 +145,6 @@ let rs = @reaction_network begin (k,t), 0 <--> X end - ps = [k => 1.0] + ps = [:k => 1.0] @test_throws Exception hc_steady_states(rs, ps) end \ No newline at end of file From ac3d9ccbf0409fe2e74da17ec8380cee958f72b3 Mon Sep 17 00:00:00 2001 From: Torkel Date: Mon, 4 Dec 2023 13:15:25 -0500 Subject: [PATCH 07/28] up --- .../steady_state_stability_computation.md | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename docs/src/{catalyst_functionality => catalyst_applications}/steady_state_stability_computation.md (100%) diff --git a/docs/src/catalyst_functionality/steady_state_stability_computation.md b/docs/src/catalyst_applications/steady_state_stability_computation.md similarity index 100% rename from docs/src/catalyst_functionality/steady_state_stability_computation.md rename to docs/src/catalyst_applications/steady_state_stability_computation.md From 9eeead55c36f1cf36e4bd36d145eb110bce8ff9c Mon Sep 17 00:00:00 2001 From: Torkel Date: Mon, 4 Dec 2023 17:24:44 -0500 Subject: [PATCH 08/28] up --- .../homotopy_continuation_extension.jl | 8 -------- 1 file changed, 8 deletions(-) diff --git a/ext/CatalystHomotopyContinuationExtension/homotopy_continuation_extension.jl b/ext/CatalystHomotopyContinuationExtension/homotopy_continuation_extension.jl index ca4e6408c4..a5355a6aaf 100644 --- a/ext/CatalystHomotopyContinuationExtension/homotopy_continuation_extension.jl +++ b/ext/CatalystHomotopyContinuationExtension/homotopy_continuation_extension.jl @@ -58,14 +58,6 @@ function steady_state_polynomial(rs::ReactionSystem, ps, u0) return poly_type_convert(ss_poly) end -# If u0s are not given while conservation laws are present, throws an error. -function conservationlaw_errorcheck(rs, pre_varmap) - vars_with_vals = Set(p[1] for p in pre_varmap) - any(s -> s in vars_with_vals, species(rs)) && return - isempty(conservedequations(rs)) || - error("The system has conservation laws but initial conditions were not provided for some species.") -end - # Parses and expression and return a version where any exponents that are Float64 (but an int, like 2.0) are turned into Int64s. make_int_exps(expr) = wrap(Rewriters.Postwalk(Rewriters.PassThrough(___make_int_exps))(unwrap(expr))).val function ___make_int_exps(expr) From 69f6f63c99ebd4693515d0c25dd93a39b4c491fc Mon Sep 17 00:00:00 2001 From: Torkel Date: Mon, 4 Dec 2023 18:26:38 -0500 Subject: [PATCH 09/28] add some additional tests --- test/miscellaneous_tests/stability_computation.jl | 13 +++++++++++++ 1 file changed, 13 insertions(+) diff --git a/test/miscellaneous_tests/stability_computation.jl b/test/miscellaneous_tests/stability_computation.jl index ac8054e762..d8ad4f9521 100644 --- a/test/miscellaneous_tests/stability_computation.jl +++ b/test/miscellaneous_tests/stability_computation.jl @@ -30,17 +30,23 @@ let # Generates random parameter values (which can generate all steady states cases). p = [:v => 1.0 + 3*rand(rng), :K => 0.5 + 2*rand(rng), :n => rand(rng,[1,2,3,4]), :d => 0.5 + rand(rng)] + # Computes stability using various jacobian options. sss = hc_steady_states(rn, p) stabs_1 = steady_state_stability(sss, rn, p) stabs_2 = steady_state_stability(sss, rn, p; sparse=true) stabs_3 = steady_state_stability(sss, rn, p; ss_jac=ss_jac) stabs_4 = steady_state_stability(sss, rn, p; ss_jac=ss_jac_sparse) + # Confirms stability using simulations. for (idx,ss) in enumerate(sss) oprob = ODEProblem(rn, [1.001, 0.999] .* ss, (0.0,1000.0), p) sol_end = solve(oprob, Rosenbrock23())[end] stabs_5 = ss ≈ sol_end @test stabs_1[idx] == stabs_2[idx] == stabs_3[idx] == stabs_4[idx] == stabs_5 + + # Checks stability when steady state is given on a pair form ([:X => x_val, :Y => y_val]). + stabs_6 = steady_state_stability(Pair.(states(rn),ss), rn, p) + @test stabs_5 == stabs_6 end end end @@ -49,6 +55,7 @@ end # Tests for system with conservation laws. # Tests for various input forms of u0 and ps. let + # Creates model. rn = complete(@reaction_network begin k1+Z, Y --> 2X k2, 2X --> X + Y @@ -56,6 +63,8 @@ let k4, X --> 0 (kD1+X, kD2), 2Z <--> Z2 end) + + # Creates various forms of input. @unpack k1, k2, k3, k4, kD1, kD2, X, Y, Z, Z2 = rn u0_1 = [X => 1.0, Y => 1.0, Z => 1.0, Z2 => 1.0] u0_2 = [:X => 1.0, :Y => 1.0, :Z => 1.0, :Z2 => 1.0] @@ -66,6 +75,7 @@ let ps_3 = [rn.k1 => 8.0, rn.k2 => 2.0, rn.k3 => 1.0, rn.k4 => 1.5, rn.kD1 => 0.5, rn.kD2 => 4.0] ps_4 = [8.0, 2.0, 1.0, 1.5, 0.5, 4.0] + # Computes stability using various input forms, and checks that the output is correct. sss = hc_steady_states(rn, ps_1; u0=u0_1) for u0 in [u0_1, u0_2, u0_3, u0_4], ps in [ps_1, ps_2, ps_3, ps_4] stab_1 = steady_state_stability(sss, rn, ps) @@ -77,4 +87,7 @@ let @test length(stab_2) == 3 @test count(stab_2) == 2 end + + # Confirms error when computing Jacobian with wrong length of u0. + @test_throws Exception steady_state_jac(rn; u0=[1.0, 1.0]) end \ No newline at end of file From 75bd7cf20c49cb68b6fa42c866a3484e8b1eb8d1 Mon Sep 17 00:00:00 2001 From: Torkel Date: Wed, 6 Dec 2023 16:51:57 -0500 Subject: [PATCH 10/28] doc up --- .../steady_state_stability_computation.md | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/docs/src/catalyst_applications/steady_state_stability_computation.md b/docs/src/catalyst_applications/steady_state_stability_computation.md index 4f812ce27d..dd650a1491 100644 --- a/docs/src/catalyst_applications/steady_state_stability_computation.md +++ b/docs/src/catalyst_applications/steady_state_stability_computation.md @@ -33,6 +33,9 @@ Next, we can apply `steady_state_stability` directly to this steady state vector steady_state_stability(steady_states, sa_loop, ps) ``` +!!! 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. @@ -48,7 +51,4 @@ 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 -``` - -!!! 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. \ No newline at end of file +``` \ No newline at end of file From 4a890aff1e4c2130455851ab81bab2de934f6f47 Mon Sep 17 00:00:00 2001 From: Torkel Date: Sat, 9 Dec 2023 16:56:17 -0500 Subject: [PATCH 11/28] doc modification --- docs/src/steady_state_functionality/nonlinear_solve.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/src/steady_state_functionality/nonlinear_solve.md b/docs/src/steady_state_functionality/nonlinear_solve.md index 0d4a903714..acf20d56c9 100644 --- a/docs/src/steady_state_functionality/nonlinear_solve.md +++ b/docs/src/steady_state_functionality/nonlinear_solve.md @@ -66,7 +66,7 @@ Like previously, the found steady state is unaffected by the initial `u` guess. ## [Systems with conservation laws](@id nonlinear_solve_conservation_laws) -As described in the section on homotopy continuation, when finding the steady states of systems with conservation laws, [additional considerations have to be taken](homotopy_continuation_conservation_laws). E.g. consider the following two-state system: +As described in the section on homotopy continuation, when finding the steady states of systems with conservation laws, [additional considerations have to be taken](@ref homotopy_continuation_conservation_laws). E.g. consider the following two-state system: ```@example nonlinear_solve2 using Catalyst, NonlinearSolve # hide two_state_model = @reaction_network begin From 0d50188f16a5fee966aaddaf1925bcb6e07ea428 Mon Sep 17 00:00:00 2001 From: Torkel Date: Fri, 2 Feb 2024 10:33:55 -0500 Subject: [PATCH 12/28] merge fix --- src/steady_state_stability.jl | 11 +- test/extensions/bifurcation_kit.jl | 15 ++- .../simulation_and_solving/solve_nonlinear.jl | 103 ------------------ 3 files changed, 20 insertions(+), 109 deletions(-) delete mode 100644 test/simulation_and_solving/solve_nonlinear.jl diff --git a/src/steady_state_stability.jl b/src/steady_state_stability.jl index dd197091e9..c3f3f213a8 100644 --- a/src/steady_state_stability.jl +++ b/src/steady_state_stability.jl @@ -1,7 +1,7 @@ ### 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_war=true) + 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) Compute the stability of a steady state (Returned as a `Bool`, with `true` indicating stability). @@ -11,8 +11,8 @@ Arguments: - `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_war=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`. + - `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`. Example: ```julia @@ -31,9 +31,10 @@ Notes: y states (each being a `Vector`) is provided as `u`, stability is computed for each state separately. """ 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_war=true) where T + sparse=false, ss_jac = steady_state_jac(rs; u0=u, sparse=sparse), + t = Inf, non_autonomous_warn =true) where T # Warning checks. - !is_autonomous(rs) && non_autonomous_war && @warn "Attempting to compute stability for a non-autonomous system. Set `non_autonomous_war=false` to disable this warning." + !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. diff --git a/test/extensions/bifurcation_kit.jl b/test/extensions/bifurcation_kit.jl index 5435775840..7900ed436b 100644 --- a/test/extensions/bifurcation_kit.jl +++ b/test/extensions/bifurcation_kit.jl @@ -174,7 +174,7 @@ let @test all(1 ./ k1s .* (1 .- xs) .≈ xs) # Checks that there is an error if information for conserved quantities computation is not provided. - @test_throws Exception bprob = BifurcationProblem(rn, u_guess, p_start, k1; plot_var = X1) + @test_throws Exception BifurcationProblem(rn, u_guess, p_start, k1; plot_var = X1) end @@ -214,4 +214,17 @@ let # Computes bifurcation diagram. @test_throws Exception BifurcationProblem(incomplete_network, u0_guess, p_start, :p) +end + +# Tests that computation for non-autonomous systems yields appropriate errors. +let + # Create t-dependant model. + rn = @reaction_network begin + (p/t,d), 0 <--> X + end + u0_guess = [:X => 1.0] + p_start = [:p => 1.0, :d => 0.2] + + # Attempts to build a BifurcationProblem. + @test_throws Exception BifurcationProblem(rn, u0_guess, p_start, :p) end \ No newline at end of file diff --git a/test/simulation_and_solving/solve_nonlinear.jl b/test/simulation_and_solving/solve_nonlinear.jl deleted file mode 100644 index 1dcb8cd5c1..0000000000 --- a/test/simulation_and_solving/solve_nonlinear.jl +++ /dev/null @@ -1,103 +0,0 @@ -### Prepares Tests ### - -# Fetch packages. -using Catalyst, NonlinearSolve, OrdinaryDiffEq, SteadyStateDiffEq -using Random, Test - -# Sets stable rng number. -using StableRNGs -rng = StableRNG(12345) - -# Fetch test functions. -include("../test_functions.jl") - -### Basic Tests ### - -# Creates a simple problem and find steady states just different approaches. -# Compares to analytic solution. -let - # Model with easily computable steady states. - steady_state_network_1 = @reaction_network begin - (k1, k2), ∅ ↔ X1 - (k3, k4), ∅ ↔ 3X2 - (k5, k6*X1), ∅ ↔ X3 - end - - # Creates NonlinearProblem. - u0 = rnd_u0(steady_state_network_1, rng) - ps = rnd_ps(steady_state_network_1, rng) - nlprob = NonlinearProblem(steady_state_network_1, u0, ps) - - # Solves it using standard algorithm and simulation based algorithm. - sol1 = solve(nlprob; abstol=1e-12, reltol=1e-12) - sol2 = solve(nlprob, DynamicSS(Rosenbrock23()); abstol=1e-12, reltol=1e-12) - - # Tests solutions are correct. - @test sol1[:X1] ≈ nlprob.ps[:k1] / nlprob.ps[:k2] atol=1e-10 - @test sol1[:X2]^3 / factorial(3) ≈ nlprob.ps[:k3] / nlprob.ps[:k4] atol=1e-10 - @test sol1[:X3] ≈ (nlprob.ps[:k5] * nlprob.ps[:k2]) / (nlprob.ps[:k6] * nlprob.ps[:k1]) atol=1e-10 - @test sol2[:X1] ≈ nlprob.ps[:k1] / nlprob.ps[:k2] atol=1e-10 - @test sol2[:X2]^3 / factorial(3) ≈ nlprob.ps[:k3] / nlprob.ps[:k4] atol=1e-10 - @test sol2[:X3] ≈ (nlprob.ps[:k5] * nlprob.ps[:k2]) / (nlprob.ps[:k6] * nlprob.ps[:k1]) atol=1e-10 -end - -# Creates a system with multiple steady states. -# Checks that corresponding ODEFunction return 0.0 in both cases. -# Checks for manually computed function as well. -let - # Creates steady state network. - steady_state_network_2 = @reaction_network begin - v / 10 + hill(X, v, K, n), ∅ → X - d, X → ∅ - end - - # Creates NonlinearProblem. - u0 = rnd_u0(steady_state_network_2, rng; min = 1.0) - ps = rnd_ps(steady_state_network_2, rng) - nlprob = NonlinearProblem(steady_state_network_2, u0, ps) - - # Solves it using standard algorithm and simulation based algorithm. - sol1 = solve(nlprob; abstol=1e-12, reltol=1e-12) - sol2 = solve(nlprob, DynamicSS(Rosenbrock23()); abstol=1e-12, reltol=1e-12) - - # Computes NonlinearFunction (manually and automatically). - nlprob_eval_1 = remake(nlprob; u0 = [:X => sol1[1]]) - nlprob_eval_2 = remake(nlprob; u0 = [:X => sol2[1]]) - function nf_manual(u,p) - X = u[:X] - v, K, n, d = p - return v/10 + v * X^n/(X^n + K^n) - d*X - end - - # Tests solutions are correct. - @test nlprob_eval_1.f(nlprob_eval_1.u0, nlprob_eval_1.p)[1] ≈ 0.0 atol=1e-10 - @test nlprob_eval_2.f(nlprob_eval_2.u0, nlprob_eval_2.p)[1] ≈ 0.0 atol=1e-10 - @test nf_manual(sol1, last.(ps)) ≈ 0.0 atol=1e-10 - @test nf_manual(sol2, last.(ps)) ≈ 0.0 atol=1e-10 -end - -# Checks for system with conservation laws. -# Checks using interfacing with output solution. -let - # Creates steady state network, unpack the parameter values. - steady_state_network_3 = @reaction_network begin - (p,d), 0 <--> X - (k1, k2), 2Y <--> Y2 - (k3, k4), X + Y2 <--> XY2 - end - @unpack X, Y, Y2, XY2 = steady_state_network_3 - - # Creates NonlinearProblem. - u0 = [steady_state_network_3.X => rand(), steady_state_network_3.Y => rand() + 1.0, steady_state_network_3.Y2 => rand() + 3.0, steady_state_network_3.XY2 => 0.0] - p = [:p => rand()+1.0, :d => 0.5, :k1 => 1.0, :k2 => 2.0, :k3 => 3.0, :k4 => 4.0] - nl_prob_1 = NonlinearProblem(steady_state_network_3, u0, p; remove_conserved = true) - nl_prob_2 = NonlinearProblem(steady_state_network_3, u0, p) - - # Solves it using standard algorithm and simulation based algorithm. - sol1 = solve(nl_prob_1; abstol=1e-12, reltol=1e-12) - sol2 = solve(nl_prob_2, DynamicSS(Rosenbrock23()); abstol=1e-12, reltol=1e-12) - - # Checks output using the ODE's drift function - @test f_eval(steady_state_network_3, [:X => sol1[X], :Y => sol1[Y], :Y2 => sol1[Y2], :XY2 => sol1[XY2]], p, 0.0) ≈ [0.0, 0.0, 0.0, 0.0] atol=1e-10 - @test f_eval(steady_state_network_3, [:X => sol2[X], :Y => sol2[Y], :Y2 => sol2[Y2], :XY2 => sol2[XY2]], p, 0.0) ≈ [0.0, 0.0, 0.0, 0.0] atol=1e-10 -end \ No newline at end of file From 1f72c7e713504d7bf45747b2a4c909f97b318e7f Mon Sep 17 00:00:00 2001 From: Torkel Date: Mon, 20 May 2024 14:59:57 -0400 Subject: [PATCH 13/28] rebase --- docs/pages.jl | 2 +- .../steady_state_stability_computation.md | 9 ++ .../bifurcation_kit_extension.jl | 4 +- .../homotopy_continuation_extension.jl | 4 +- src/Catalyst.jl | 3 +- src/reactionsystem.jl | 29 +++++ src/reactionsystem_conversions.jl | 6 + src/steady_state_stability.jl | 107 ++++++++++-------- test/miscellaneous_tests/api.jl | 9 +- .../stability_computation.jl | 60 +++++----- test/runtests.jl | 1 + .../simulation_and_solving/solve_nonlinear.jl | 103 +++++++++++++++++ 12 files changed, 251 insertions(+), 86 deletions(-) rename docs/src/{catalyst_applications => steady_state_functionality}/steady_state_stability_computation.md (83%) create mode 100644 test/simulation_and_solving/solve_nonlinear.jl diff --git a/docs/pages.jl b/docs/pages.jl index 04f255eede..c0358a9613 100644 --- a/docs/pages.jl +++ b/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. ], diff --git a/docs/src/catalyst_applications/steady_state_stability_computation.md b/docs/src/steady_state_functionality/steady_state_stability_computation.md similarity index 83% rename from docs/src/catalyst_applications/steady_state_stability_computation.md rename to docs/src/steady_state_functionality/steady_state_stability_computation.md index dd650a1491..18f2ad0c47 100644 --- a/docs/src/catalyst_applications/steady_state_stability_computation.md +++ b/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 ``` \ No newline at end of file diff --git a/ext/CatalystBifurcationKitExtension/bifurcation_kit_extension.jl b/ext/CatalystBifurcationKitExtension/bifurcation_kit_extension.jl index 7c940190df..3e2ea2effc 100644 --- a/ext/CatalystBifurcationKitExtension/bifurcation_kit_extension.jl +++ b/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]) diff --git a/ext/CatalystHomotopyContinuationExtension/homotopy_continuation_extension.jl b/ext/CatalystHomotopyContinuationExtension/homotopy_continuation_extension.jl index a5355a6aaf..926e7ff0eb 100644 --- a/ext/CatalystHomotopyContinuationExtension/homotopy_continuation_extension.jl +++ b/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) diff --git a/src/Catalyst.jl b/src/Catalyst.jl index f87548b47f..33ee2b87fb 100644 --- a/src/Catalyst.jl +++ b/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 diff --git a/src/reactionsystem.jl b/src/reactionsystem.jl index 47a0af56b4..7292042b7f 100644 --- a/src/reactionsystem.jl +++ b/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 ### diff --git a/src/reactionsystem_conversions.jl b/src/reactionsystem_conversions.jl index f122d3f62d..b7c7e656f3 100644 --- a/src/reactionsystem_conversions.jl +++ b/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) diff --git a/src/steady_state_stability.jl b/src/steady_state_stability.jl index c3f3f213a8..95add620f9 100644 --- a/src/steady_state_stability.jl +++ b/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 \ No newline at end of file diff --git a/test/miscellaneous_tests/api.jl b/test/miscellaneous_tests/api.jl index 1340f5d85e..ccf0ce2d26 100644 --- a/test/miscellaneous_tests/api.jl +++ b/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 \ No newline at end of file diff --git a/test/miscellaneous_tests/stability_computation.jl b/test/miscellaneous_tests/stability_computation.jl index d8ad4f9521..1cdc3304b7 100644 --- a/test/miscellaneous_tests/stability_computation.jl +++ b/test/miscellaneous_tests/stability_computation.jl @@ -1,8 +1,7 @@ ### Fetch Packages ### # Fetch packages. -using Catalyst, OrdinaryDiffEq -using Random, Test +using Catalyst, OrdinaryDiffEq, SteadyStateDiffEq, Test import HomotopyContinuation # Sets rnd number. @@ -14,7 +13,7 @@ rng = StableRNG(12345) # Tests that stability is correctly assessed (using simulation) in multi stable system. # Tests that `steady_state_jac` function works. # Tests with and without sparsity. -# tests using symbolic input. +# Tests using symbolic input. let # System which may have between 1 and 7 fixed points. rn = @reaction_network begin @@ -23,29 +22,30 @@ let d, (X,Y) --> 0 end ss_jac = steady_state_jac(rn) - ss_jac_sparse = steady_state_jac(rn; sparse=true) + ss_jac_sparse = steady_state_jac(rn; sparse = true) - # Repeats several times, most cases should be encountered several times. - for i = 1:50 + # Repeats several times, most steady state stability cases should be encountered several times. + for repeat = 1:50 # Generates random parameter values (which can generate all steady states cases). - p = [:v => 1.0 + 3*rand(rng), :K => 0.5 + 2*rand(rng), :n => rand(rng,[1,2,3,4]), :d => 0.5 + rand(rng)] + ps = (:v => 1.0 + 3*rand(rng), :K => 0.5 + 2*rand(rng), :n => rand(rng,[1,2,3,4]), + :d => 0.5 + rand(rng)) # Computes stability using various jacobian options. - sss = hc_steady_states(rn, p) - stabs_1 = steady_state_stability(sss, rn, p) - stabs_2 = steady_state_stability(sss, rn, p; sparse=true) - stabs_3 = steady_state_stability(sss, rn, p; ss_jac=ss_jac) - stabs_4 = steady_state_stability(sss, rn, p; ss_jac=ss_jac_sparse) + sss = hc_steady_states(rn, ps) + stabs_1 = [steady_state_stability(ss, rn, ps) for ss in sss] + stabs_2 = [steady_state_stability(ss, rn, ps; sparse = true) for ss in sss] + stabs_3 = [steady_state_stability(ss, rn, ps; ss_jac = ss_jac) for ss in sss] + stabs_4 = [steady_state_stability(ss, rn, ps; ss_jac = ss_jac_sparse) for ss in sss] # Confirms stability using simulations. - for (idx,ss) in enumerate(sss) - oprob = ODEProblem(rn, [1.001, 0.999] .* ss, (0.0,1000.0), p) - sol_end = solve(oprob, Rosenbrock23())[end] - stabs_5 = ss ≈ sol_end + for (idx, ss) in enumerate(sss) + ssprob = SteadyStateProblem(rn, [1.001, 0.999] .* ss, ps) + sol = solve(ssprob, DynamicSS(Vern7()); abstol = 1e-8, reltol = 1e-8) + stabs_5 = isapprox(ss, sol.u; atol = 1e-6) @test stabs_1[idx] == stabs_2[idx] == stabs_3[idx] == stabs_4[idx] == stabs_5 - # Checks stability when steady state is given on a pair form ([:X => x_val, :Y => y_val]). - stabs_6 = steady_state_stability(Pair.(states(rn),ss), rn, p) + # Checks stability when steady state is given on a pair form ([X => x_val, Y => y_val]). + stabs_6 = steady_state_stability(Pair.(unknowns(rn), ss), rn, ps) @test stabs_5 == stabs_6 end end @@ -56,13 +56,13 @@ end # Tests for various input forms of u0 and ps. let # Creates model. - rn = complete(@reaction_network begin + rn = @reaction_network begin k1+Z, Y --> 2X k2, 2X --> X + Y k3, X + Y --> Y k4, X --> 0 (kD1+X, kD2), 2Z <--> Z2 - end) + end # Creates various forms of input. @unpack k1, k2, k3, k4, kD1, kD2, X, Y, Z, Z2 = rn @@ -73,21 +73,17 @@ let ps_1 = [k1 => 8.0, k2 => 2.0, k3 => 1.0, k4 => 1.5, kD1 => 0.5, kD2 => 2.0] ps_2 = [:k1 => 8.0, :k2 => 2.0, :k3 => 1.0, :k4 => 1.5, :kD1 => 0.5, :kD2 => 2.0] ps_3 = [rn.k1 => 8.0, rn.k2 => 2.0, rn.k3 => 1.0, rn.k4 => 1.5, rn.kD1 => 0.5, rn.kD2 => 4.0] - ps_4 = [8.0, 2.0, 1.0, 1.5, 0.5, 4.0] # Computes stability using various input forms, and checks that the output is correct. - sss = hc_steady_states(rn, ps_1; u0=u0_1) - for u0 in [u0_1, u0_2, u0_3, u0_4], ps in [ps_1, ps_2, ps_3, ps_4] - stab_1 = steady_state_stability(sss, rn, ps) - @test length(stab_1) == 3 - @test count(stab_1) == 2 - - ss_jac = steady_state_jac(rn; u0=u0) - stab_2 = steady_state_stability(sss, rn, ps; ss_jac=ss_jac) - @test length(stab_2) == 3 - @test count(stab_2) == 2 + sss = hc_steady_states(rn, ps_1; u0 = u0_1) + for u0 in [u0_1, u0_2, u0_3, u0_4], ps in [ps_1, ps_2, ps_3] + stab_1 = [steady_state_stability(ss, rn, ps) for ss in sss] + ss_jac = steady_state_jac(rn; u0 = u0) + stab_2 = [steady_state_stability(ss, rn, ps; ss_jac = ss_jac) for ss in sss] + @test length(stab_1) == length(stab_2) == 3 + @test count(stab_1) == count(stab_2) == 2 end # Confirms error when computing Jacobian with wrong length of u0. - @test_throws Exception steady_state_jac(rn; u0=[1.0, 1.0]) + @test_throws Exception steady_state_jac(rn; u0 = [1.0, 1.0]) end \ No newline at end of file diff --git a/test/runtests.jl b/test/runtests.jl index 744ade8cdb..d7953daf40 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -33,6 +33,7 @@ using SafeTestsets, Test #if GROUP == "All" || GROUP == "Miscellaneous-NetworkAnalysis" # Tests various miscellaneous features. @time @safetestset "API" begin include("miscellaneous_tests/api.jl") end + @time @safetestset "Steady State Stability Computations" begin include("miscellaneous_tests/stability_computation.jl") end @time @safetestset "Compound Species" begin include("miscellaneous_tests/compound_macro.jl") end @time @safetestset "Reaction Balancing" begin include("miscellaneous_tests/reaction_balancing.jl") end @time @safetestset "Units" begin include("miscellaneous_tests/units.jl") end diff --git a/test/simulation_and_solving/solve_nonlinear.jl b/test/simulation_and_solving/solve_nonlinear.jl new file mode 100644 index 0000000000..1dcb8cd5c1 --- /dev/null +++ b/test/simulation_and_solving/solve_nonlinear.jl @@ -0,0 +1,103 @@ +### Prepares Tests ### + +# Fetch packages. +using Catalyst, NonlinearSolve, OrdinaryDiffEq, SteadyStateDiffEq +using Random, Test + +# Sets stable rng number. +using StableRNGs +rng = StableRNG(12345) + +# Fetch test functions. +include("../test_functions.jl") + +### Basic Tests ### + +# Creates a simple problem and find steady states just different approaches. +# Compares to analytic solution. +let + # Model with easily computable steady states. + steady_state_network_1 = @reaction_network begin + (k1, k2), ∅ ↔ X1 + (k3, k4), ∅ ↔ 3X2 + (k5, k6*X1), ∅ ↔ X3 + end + + # Creates NonlinearProblem. + u0 = rnd_u0(steady_state_network_1, rng) + ps = rnd_ps(steady_state_network_1, rng) + nlprob = NonlinearProblem(steady_state_network_1, u0, ps) + + # Solves it using standard algorithm and simulation based algorithm. + sol1 = solve(nlprob; abstol=1e-12, reltol=1e-12) + sol2 = solve(nlprob, DynamicSS(Rosenbrock23()); abstol=1e-12, reltol=1e-12) + + # Tests solutions are correct. + @test sol1[:X1] ≈ nlprob.ps[:k1] / nlprob.ps[:k2] atol=1e-10 + @test sol1[:X2]^3 / factorial(3) ≈ nlprob.ps[:k3] / nlprob.ps[:k4] atol=1e-10 + @test sol1[:X3] ≈ (nlprob.ps[:k5] * nlprob.ps[:k2]) / (nlprob.ps[:k6] * nlprob.ps[:k1]) atol=1e-10 + @test sol2[:X1] ≈ nlprob.ps[:k1] / nlprob.ps[:k2] atol=1e-10 + @test sol2[:X2]^3 / factorial(3) ≈ nlprob.ps[:k3] / nlprob.ps[:k4] atol=1e-10 + @test sol2[:X3] ≈ (nlprob.ps[:k5] * nlprob.ps[:k2]) / (nlprob.ps[:k6] * nlprob.ps[:k1]) atol=1e-10 +end + +# Creates a system with multiple steady states. +# Checks that corresponding ODEFunction return 0.0 in both cases. +# Checks for manually computed function as well. +let + # Creates steady state network. + steady_state_network_2 = @reaction_network begin + v / 10 + hill(X, v, K, n), ∅ → X + d, X → ∅ + end + + # Creates NonlinearProblem. + u0 = rnd_u0(steady_state_network_2, rng; min = 1.0) + ps = rnd_ps(steady_state_network_2, rng) + nlprob = NonlinearProblem(steady_state_network_2, u0, ps) + + # Solves it using standard algorithm and simulation based algorithm. + sol1 = solve(nlprob; abstol=1e-12, reltol=1e-12) + sol2 = solve(nlprob, DynamicSS(Rosenbrock23()); abstol=1e-12, reltol=1e-12) + + # Computes NonlinearFunction (manually and automatically). + nlprob_eval_1 = remake(nlprob; u0 = [:X => sol1[1]]) + nlprob_eval_2 = remake(nlprob; u0 = [:X => sol2[1]]) + function nf_manual(u,p) + X = u[:X] + v, K, n, d = p + return v/10 + v * X^n/(X^n + K^n) - d*X + end + + # Tests solutions are correct. + @test nlprob_eval_1.f(nlprob_eval_1.u0, nlprob_eval_1.p)[1] ≈ 0.0 atol=1e-10 + @test nlprob_eval_2.f(nlprob_eval_2.u0, nlprob_eval_2.p)[1] ≈ 0.0 atol=1e-10 + @test nf_manual(sol1, last.(ps)) ≈ 0.0 atol=1e-10 + @test nf_manual(sol2, last.(ps)) ≈ 0.0 atol=1e-10 +end + +# Checks for system with conservation laws. +# Checks using interfacing with output solution. +let + # Creates steady state network, unpack the parameter values. + steady_state_network_3 = @reaction_network begin + (p,d), 0 <--> X + (k1, k2), 2Y <--> Y2 + (k3, k4), X + Y2 <--> XY2 + end + @unpack X, Y, Y2, XY2 = steady_state_network_3 + + # Creates NonlinearProblem. + u0 = [steady_state_network_3.X => rand(), steady_state_network_3.Y => rand() + 1.0, steady_state_network_3.Y2 => rand() + 3.0, steady_state_network_3.XY2 => 0.0] + p = [:p => rand()+1.0, :d => 0.5, :k1 => 1.0, :k2 => 2.0, :k3 => 3.0, :k4 => 4.0] + nl_prob_1 = NonlinearProblem(steady_state_network_3, u0, p; remove_conserved = true) + nl_prob_2 = NonlinearProblem(steady_state_network_3, u0, p) + + # Solves it using standard algorithm and simulation based algorithm. + sol1 = solve(nl_prob_1; abstol=1e-12, reltol=1e-12) + sol2 = solve(nl_prob_2, DynamicSS(Rosenbrock23()); abstol=1e-12, reltol=1e-12) + + # Checks output using the ODE's drift function + @test f_eval(steady_state_network_3, [:X => sol1[X], :Y => sol1[Y], :Y2 => sol1[Y2], :XY2 => sol1[XY2]], p, 0.0) ≈ [0.0, 0.0, 0.0, 0.0] atol=1e-10 + @test f_eval(steady_state_network_3, [:X => sol2[X], :Y => sol2[Y], :Y2 => sol2[Y2], :XY2 => sol2[XY2]], p, 0.0) ≈ [0.0, 0.0, 0.0, 0.0] atol=1e-10 +end \ No newline at end of file From 5da15fedcddc91a6e6fd8e2fb92431f1cc4d1a23 Mon Sep 17 00:00:00 2001 From: Torkel Date: Mon, 20 May 2024 18:27:01 -0400 Subject: [PATCH 14/28] remove nonliear test on non-autonomous system --- test/reactionsystem_core/reactionsystem.jl | 2 -- 1 file changed, 2 deletions(-) diff --git a/test/reactionsystem_core/reactionsystem.jl b/test/reactionsystem_core/reactionsystem.jl index 2bdb3df42e..53a93eb66f 100644 --- a/test/reactionsystem_core/reactionsystem.jl +++ b/test/reactionsystem_core/reactionsystem.jl @@ -118,13 +118,11 @@ let odesys = complete(convert(ODESystem, rs)) sdesys = complete(convert(SDESystem, rs)) js = complete(convert(JumpSystem, rs)) - nlsys = complete(convert(NonlinearSystem, rs)) @test ModelingToolkit.get_defaults(rs) == ModelingToolkit.get_defaults(odesys) == ModelingToolkit.get_defaults(sdesys) == ModelingToolkit.get_defaults(js) == - ModelingToolkit.get_defaults(nlsys) == defs u0map = [A => 5.0] From 20cbc75a569d6076acd0eb7af7a36ae5b4509562 Mon Sep 17 00:00:00 2001 From: Torkel Date: Mon, 20 May 2024 19:59:48 -0400 Subject: [PATCH 15/28] up --- test/reactionsystem_core/reactionsystem.jl | 20 +++++++++++--------- 1 file changed, 11 insertions(+), 9 deletions(-) diff --git a/test/reactionsystem_core/reactionsystem.jl b/test/reactionsystem_core/reactionsystem.jl index 53a93eb66f..24a0d0ad70 100644 --- a/test/reactionsystem_core/reactionsystem.jl +++ b/test/reactionsystem_core/reactionsystem.jl @@ -156,15 +156,17 @@ let G2 = sf.g(sprob.u0, sprob.p, 0.0) @test norm(G - G2) < 100 * eps() - # Test conversion to NonlinearSystem. - ns = convert(NonlinearSystem, rs) - nlprob = NonlinearProblem(rs, u, p) - fnl = eval(generate_function(ns)[2]) - dunl = similar(du) - @test_broken let # The next line throws an error. - fnl(dunl, nlprob.u0, nlprob.p) - @test norm(du - dunl) < 100 * eps() - end + # Test conversion to NonlinearSystem. Permanently broken as we no longer permit non-autonomous + # `ReactionSystem`s to be converted to `NonlinearSystem`s. + @test_broken false + # ns = convert(NonlinearSystem, rs) + # nlprob = NonlinearProblem(rs, u, p) + # fnl = eval(generate_function(ns)[2]) + # dunl = similar(du) + # @test_broken let # The next line throws an error. + # fnl(dunl, nlprob.u0, nlprob.p) + # @test norm(du - dunl) < 100 * eps() + # end end # Test with JumpSystem. From 9c45cd25e6581f72e843a0cf82c4cab3157a670c Mon Sep 17 00:00:00 2001 From: Torkel Date: Mon, 20 May 2024 21:28:06 -0400 Subject: [PATCH 16/28] fix --- src/reactionsystem.jl | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/src/reactionsystem.jl b/src/reactionsystem.jl index 7292042b7f..976ed1fd77 100644 --- a/src/reactionsystem.jl +++ b/src/reactionsystem.jl @@ -1224,7 +1224,11 @@ 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]) + if isempty(dep_var_param_rxs) && isempty(dep_var_param_eqs) + dep_var_param = [] + else + dep_var_param = reduce(vcat,[dep_var_param_rxs; dep_var_param_eqs]) + end # Checks for iv and spatial ivs any(isequal(get_iv(rs), var) for var in dep_var_param) && (return false) From 522b6fa64cdde949a52bb8104236cdeafbbe4d92 Mon Sep 17 00:00:00 2001 From: Torkel Date: Tue, 21 May 2024 11:30:22 -0400 Subject: [PATCH 17/28] add tolerances --- .../steady_state_stability_computation.md | 2 +- src/steady_state_stability.jl | 55 +++++++++++++------ 2 files changed, 40 insertions(+), 17 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 18f2ad0c47..7f84c78899 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 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. + 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: diff --git a/src/steady_state_stability.jl b/src/steady_state_stability.jl index 95add620f9..da65b9bc38 100644 --- a/src/steady_state_stability.jl +++ b/src/steady_state_stability.jl @@ -1,19 +1,25 @@ ### Stability Analysis ### """ - stability(u::Vector{T}, rs::ReactionSystem, p; sparse = false, - ss_jac = steady_state_jac(u, rs, p; sparse=sparse)) + stability(u::Vector{T}, rs::ReactionSystem, p; tol = 10*sqrt(eps()) + 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 `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. +- `u`: The steady state 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. +- `tol = 10*sqrt(eps())`: The tolerance used for determining whether computed eigenvalues real +parts can reliably be considered negative/positive. In practise, a steady state is considered +stable if the corresponding Jacobian's maximum eigenvalue real part is < 0. However, if this maximum +eigenvalue is in the range `-tol< eig < tol`, and error is throw, as we do not deem that stability +can be ensured with enough certainty. The choice `tol = 10*sqrt(eps())` has *not* been subject +to much analysis. +- `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 @@ -28,14 +34,20 @@ steady_state = [2.0] steady_state_stability(steady_state, rn, p) Notes: - - 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. +- 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. +- 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). ``` """ -function steady_state_stability(u::Vector, rs::ReactionSystem, ps; sparse = false, - ss_jac = steady_state_jac(rs; u0 = u, sparse = sparse)) +function steady_state_stability(u::Vector, rs::ReactionSystem, ps; tol = 10*sqrt(eps(ss_val_type(u))), + sparse = false, ss_jac = steady_state_jac(rs; u0 = u, sparse = sparse)) # Warning checks. 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.") @@ -53,9 +65,20 @@ function steady_state_stability(u::Vector, rs::ReactionSystem, ps; sparse = fals 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 + max_eig = maximum(real(ev) for ev in (eigvals(J))) + if abs(max_eig) < tol + error("The system Jacobian's maximum eigenvalue at the steady state is within the tolerance range (abs($max_eig) < $tol). Hence, stability could not be reliably determined. If you still wish to compute the stability, reduce the `tol` argument.") + end + return max_eig < 0 end +# Used to determine the type of the steady states values, which is then used to set the tolerance's +# type. +ss_val_type(u::Vector{T}) where {T} = T +ss_val_type(u::Vector{Pair{S,T}}) where {S,T} = T +ss_val_type(u::Dict{S,T}) where {S,T} = T + + """ steady_state_jac(rs::ReactionSystem; u0 = [], sparse = false) From 834af284c4a0b787e981e0daa1dc3cc5fdfe44e6 Mon Sep 17 00:00:00 2001 From: Torkel Date: Tue, 21 May 2024 11:48:07 -0400 Subject: [PATCH 18/28] add tolerance test --- src/steady_state_stability.jl | 1 - .../stability_computation.jl | 17 +++++++++++++++-- 2 files changed, 15 insertions(+), 3 deletions(-) diff --git a/src/steady_state_stability.jl b/src/steady_state_stability.jl index da65b9bc38..e4b5c4a5cb 100644 --- a/src/steady_state_stability.jl +++ b/src/steady_state_stability.jl @@ -78,7 +78,6 @@ ss_val_type(u::Vector{T}) where {T} = T ss_val_type(u::Vector{Pair{S,T}}) where {S,T} = T ss_val_type(u::Dict{S,T}) where {S,T} = T - """ steady_state_jac(rs::ReactionSystem; u0 = [], sparse = false) diff --git a/test/miscellaneous_tests/stability_computation.jl b/test/miscellaneous_tests/stability_computation.jl index 1cdc3304b7..75069bc9da 100644 --- a/test/miscellaneous_tests/stability_computation.jl +++ b/test/miscellaneous_tests/stability_computation.jl @@ -8,7 +8,7 @@ import HomotopyContinuation using StableRNGs rng = StableRNG(12345) -### Run Tests ### +### Basic Tests ### # Tests that stability is correctly assessed (using simulation) in multi stable system. # Tests that `steady_state_jac` function works. @@ -86,4 +86,17 @@ let # Confirms error when computing Jacobian with wrong length of u0. @test_throws Exception steady_state_jac(rn; u0 = [1.0, 1.0]) -end \ No newline at end of file +end + +### Other Tests ### + +# Tests `tol` option. In this case, the maximum eigenvalue real part is -1.0, which generates +# and error with `tol = 100`. +let + rn = @reaction_network begin + (p,d), 0 <--> X + end + p = [:p => 1.0, :d => 1.0] + u = [1.0] + @test_throws Exception steady_state_stability(u, rn, p; tol = 1e2) +end From 9cc5f9cd63fb8994f2783b1ac05f9163928cb3f5 Mon Sep 17 00:00:00 2001 From: Torkel Date: Tue, 21 May 2024 13:29:38 -0400 Subject: [PATCH 19/28] 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. From cf46e42262084f6108c344677b6bc45b2de331f2 Mon Sep 17 00:00:00 2001 From: Torkel Date: Tue, 21 May 2024 19:11:54 -0400 Subject: [PATCH 20/28] remove option for sparsity in Jacobian (requires a special package) --- .../steady_state_stability_computation.md | 6 ------ src/steady_state_stability.jl | 19 +++++++++++-------- .../stability_computation.jl | 16 ++++++---------- 3 files changed, 17 insertions(+), 24 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 fdd12530e8..affbb5dec4 100644 --- a/docs/src/steady_state_functionality/steady_state_stability_computation.md +++ b/docs/src/steady_state_functionality/steady_state_stability_computation.md @@ -53,11 +53,5 @@ stability_2 = steady_state_stability(steady_states_2, sa_loop, ps_2; ss_jac=ss_j 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 -``` - !!! 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 dc901f577b..26b011cff8 100644 --- a/src/steady_state_stability.jl +++ b/src/steady_state_stability.jl @@ -2,7 +2,7 @@ """ stability(u::Vector{T}, rs::ReactionSystem, p; tol = 10*sqrt(eps()) - sparse = false, ss_jac = steady_state_jac(u, rs, p; sparse=sparse)) + ss_jac = steady_state_jac(u, rs, p)) Compute the stability of a steady state (Returned as a `Bool`, with `true` indicating stability). @@ -16,8 +16,7 @@ stable if the corresponding Jacobian's maximum eigenvalue real part is < 0. Howe eigenvalue is in the range `-tol< eig < tol`, and error is throw, as we do not deem that stability can be ensured with enough certainty. The choice `tol = 10*sqrt(eps())` has *not* been subject to much analysis. -- `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 +- `ss_jac = steady_state_jac(u, rs)`: 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. @@ -47,7 +46,7 @@ however, have not been subject to further analysis (and can be changed through t ``` """ function steady_state_stability(u::Vector, rs::ReactionSystem, ps; tol = 10*sqrt(eps(ss_val_type(u))), - sparse = false, ss_jac = steady_state_jac(rs; u0 = u, sparse = sparse)) + ss_jac = steady_state_jac(rs; u0 = u)) # Warning checks. 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.") @@ -62,6 +61,11 @@ function steady_state_stability(u::Vector, rs::ReactionSystem, ps; tol = 10*sqrt # Computes stability (by checking that the real part of all eigenvalues are negative). # Here, `ss_jac` is a `ODEProblem` with dummy values for `u0` and `p`. + + if isdefined(Main, :Infiltrator) + Main.infiltrate(@__MODULE__, Base.@locals, @__FILE__, @__LINE__) + end + 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) @@ -79,7 +83,7 @@ ss_val_type(u::Vector{Pair{S,T}}) where {S,T} = T ss_val_type(u::Dict{S,T}) where {S,T} = T """ - steady_state_jac(rs::ReactionSystem; u0 = [], sparse = false) + steady_state_jac(rs::ReactionSystem; u0 = []) 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. @@ -87,7 +91,6 @@ a large number of stability computation has to be carried out in a performant ma 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 @@ -110,7 +113,7 @@ Notes: ``` """ function steady_state_jac(rs::ReactionSystem; u0 = [sp => 0.0 for sp in unknowns(rs)], - sparse = false, combinatoric_ratelaws = get_combinatoric_ratelaws(rs)) + combinatoric_ratelaws = get_combinatoric_ratelaws(rs)) # If u0 is a vector of values, must be converted to something MTK understands. # Converts u0 to values MTK understands, and checks that potential conservation laws are accounted for. @@ -119,7 +122,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, remove_conserved = true, sparse = sparse, + return ODEProblem(rs, u0, 0, ps; jac = true, remove_conserved = true, combinatoric_ratelaws = combinatoric_ratelaws) end diff --git a/test/miscellaneous_tests/stability_computation.jl b/test/miscellaneous_tests/stability_computation.jl index 75069bc9da..36f6125f02 100644 --- a/test/miscellaneous_tests/stability_computation.jl +++ b/test/miscellaneous_tests/stability_computation.jl @@ -12,7 +12,6 @@ rng = StableRNG(12345) # Tests that stability is correctly assessed (using simulation) in multi stable system. # Tests that `steady_state_jac` function works. -# Tests with and without sparsity. # Tests using symbolic input. let # System which may have between 1 and 7 fixed points. @@ -22,10 +21,9 @@ let d, (X,Y) --> 0 end ss_jac = steady_state_jac(rn) - ss_jac_sparse = steady_state_jac(rn; sparse = true) # Repeats several times, most steady state stability cases should be encountered several times. - for repeat = 1:50 + for repeat = 1:20 # Generates random parameter values (which can generate all steady states cases). ps = (:v => 1.0 + 3*rand(rng), :K => 0.5 + 2*rand(rng), :n => rand(rng,[1,2,3,4]), :d => 0.5 + rand(rng)) @@ -33,20 +31,18 @@ let # Computes stability using various jacobian options. sss = hc_steady_states(rn, ps) stabs_1 = [steady_state_stability(ss, rn, ps) for ss in sss] - stabs_2 = [steady_state_stability(ss, rn, ps; sparse = true) for ss in sss] - stabs_3 = [steady_state_stability(ss, rn, ps; ss_jac = ss_jac) for ss in sss] - stabs_4 = [steady_state_stability(ss, rn, ps; ss_jac = ss_jac_sparse) for ss in sss] + stabs_2 = [steady_state_stability(ss, rn, ps; ss_jac = ss_jac) for ss in sss] # Confirms stability using simulations. for (idx, ss) in enumerate(sss) ssprob = SteadyStateProblem(rn, [1.001, 0.999] .* ss, ps) sol = solve(ssprob, DynamicSS(Vern7()); abstol = 1e-8, reltol = 1e-8) - stabs_5 = isapprox(ss, sol.u; atol = 1e-6) - @test stabs_1[idx] == stabs_2[idx] == stabs_3[idx] == stabs_4[idx] == stabs_5 + stabs_3 = isapprox(ss, sol.u; atol = 1e-6) + @test stabs_1[idx] == stabs_2[idx] == stabs_3 # Checks stability when steady state is given on a pair form ([X => x_val, Y => y_val]). - stabs_6 = steady_state_stability(Pair.(unknowns(rn), ss), rn, ps) - @test stabs_5 == stabs_6 + stabs_4 = steady_state_stability(Pair.(unknowns(rn), ss), rn, ps) + @test stabs_3 == stabs_4 end end end From 0226cc607ebba1c2a38562003e458149988200b8 Mon Sep 17 00:00:00 2001 From: Torkel Date: Tue, 21 May 2024 21:42:38 -0400 Subject: [PATCH 21/28] up --- src/steady_state_stability.jl | 5 ----- 1 file changed, 5 deletions(-) diff --git a/src/steady_state_stability.jl b/src/steady_state_stability.jl index 26b011cff8..cfcf606afe 100644 --- a/src/steady_state_stability.jl +++ b/src/steady_state_stability.jl @@ -61,11 +61,6 @@ function steady_state_stability(u::Vector, rs::ReactionSystem, ps; tol = 10*sqrt # Computes stability (by checking that the real part of all eigenvalues are negative). # Here, `ss_jac` is a `ODEProblem` with dummy values for `u0` and `p`. - - if isdefined(Main, :Infiltrator) - Main.infiltrate(@__MODULE__, Base.@locals, @__FILE__, @__LINE__) - end - 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) From b9e6aaa566d46fcc27ad36bf9b39554498ea58ed Mon Sep 17 00:00:00 2001 From: Torkel Date: Tue, 21 May 2024 22:13:18 -0400 Subject: [PATCH 22/28] up --- src/steady_state_stability.jl | 10 +++++++++- 1 file changed, 9 insertions(+), 1 deletion(-) diff --git a/src/steady_state_stability.jl b/src/steady_state_stability.jl index cfcf606afe..c7af2e00e6 100644 --- a/src/steady_state_stability.jl +++ b/src/steady_state_stability.jl @@ -64,7 +64,15 @@ function steady_state_stability(u::Vector, rs::ReactionSystem, ps; tol = 10*sqrt 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) - max_eig = maximum(real(ev) for ev in (eigvals(J))) + println() + println() + println("Got here") + println(J) + println(typeof(J)) + println(eigvals(J)) + println(typeof(eigvals(J))) + max_eig = maximum(real(ev) for ev in eigvals(J)) + println("Didn't get here") if abs(max_eig) < tol error("The system Jacobian's maximum eigenvalue at the steady state is within the tolerance range (abs($max_eig) < $tol). Hence, stability could not be reliably determined. If you still wish to compute the stability, reduce the `tol` argument.") end From f6bae74764eead92b89aa15acf5de89de2f4a6e7 Mon Sep 17 00:00:00 2001 From: Torkel Loman Date: Wed, 22 May 2024 05:01:36 +0200 Subject: [PATCH 23/28] Update steady_state_stability.jl up --- src/steady_state_stability.jl | 15 ++++----------- 1 file changed, 4 insertions(+), 11 deletions(-) diff --git a/src/steady_state_stability.jl b/src/steady_state_stability.jl index c7af2e00e6..e65e4ad921 100644 --- a/src/steady_state_stability.jl +++ b/src/steady_state_stability.jl @@ -59,20 +59,13 @@ function steady_state_stability(u::Vector, rs::ReactionSystem, ps; tol = 10*sqrt 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). - # Here, `ss_jac` is a `ODEProblem` with dummy values for `u0` and `p`. + # Generates the Jacobian at the steady state (technically, `ss_jac` is an `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) - println() - println() - println("Got here") - println(J) - println(typeof(J)) - println(eigvals(J)) - println(typeof(eigvals(J))) + + # Computes stability (by checking that the real part of all eigenvalues is negative). max_eig = maximum(real(ev) for ev in eigvals(J)) - println("Didn't get here") if abs(max_eig) < tol error("The system Jacobian's maximum eigenvalue at the steady state is within the tolerance range (abs($max_eig) < $tol). Hence, stability could not be reliably determined. If you still wish to compute the stability, reduce the `tol` argument.") end @@ -139,4 +132,4 @@ function steady_state_u_conversion(u, rs::ReactionSystem) end end return symmap_to_varmap(rs, u) -end \ No newline at end of file +end From 2bae4a4dbc011717919f336dae144e6119067d3b Mon Sep 17 00:00:00 2001 From: Torkel Loman Date: Thu, 23 May 2024 15:16:04 -0400 Subject: [PATCH 24/28] Update docs/src/steady_state_functionality/steady_state_stability_computation.md Co-authored-by: Sam Isaacson --- .../steady_state_stability_computation.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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 affbb5dec4..df791ac87f 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 tolerance `tol = 10*sqrt(eps())` to determine whether a computed eigenvalue is far away enough from 0 to be reliably used. This threshold can be changed through the `tol` keyword argument. ## Basic examples Let us consider the following basic example: From 2e56390c76ec24de94c097986d82d69371d977b1 Mon Sep 17 00:00:00 2001 From: Torkel Loman Date: Thu, 23 May 2024 15:16:17 -0400 Subject: [PATCH 25/28] Update src/steady_state_stability.jl Co-authored-by: Sam Isaacson --- src/steady_state_stability.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/steady_state_stability.jl b/src/steady_state_stability.jl index e65e4ad921..8ac4e27a83 100644 --- a/src/steady_state_stability.jl +++ b/src/steady_state_stability.jl @@ -67,7 +67,7 @@ function steady_state_stability(u::Vector, rs::ReactionSystem, ps; tol = 10*sqrt # Computes stability (by checking that the real part of all eigenvalues is negative). max_eig = maximum(real(ev) for ev in eigvals(J)) if abs(max_eig) < tol - error("The system Jacobian's maximum eigenvalue at the steady state is within the tolerance range (abs($max_eig) < $tol). Hence, stability could not be reliably determined. If you still wish to compute the stability, reduce the `tol` argument.") + error("The system Jacobian's maximum eigenvalue at the steady state is within the tolerance range (abs($max_eig) < $tol). Hence, stability could not be reliably determined. If you still wish to compute the stability, reduce the `tol` argument, but note that floating point error in the eigenvalue calculation could lead to incorrect results.") end return max_eig < 0 end From 06a2937209fe84f84bb6274fda1bb962522123f4 Mon Sep 17 00:00:00 2001 From: Torkel Loman Date: Thu, 23 May 2024 15:16:36 -0400 Subject: [PATCH 26/28] Update src/steady_state_stability.jl Co-authored-by: Sam Isaacson --- src/steady_state_stability.jl | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/steady_state_stability.jl b/src/steady_state_stability.jl index 8ac4e27a83..8f97986c51 100644 --- a/src/steady_state_stability.jl +++ b/src/steady_state_stability.jl @@ -41,8 +41,7 @@ form of maps. 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). +computed eigenvalue is far away enough from 0 to be reliably used. This selected threshold can be changed through the `tol` argument. ``` """ function steady_state_stability(u::Vector, rs::ReactionSystem, ps; tol = 10*sqrt(eps(ss_val_type(u))), From f839ff96bab948c9dcbcc97f5477f34fa5bceaeb Mon Sep 17 00:00:00 2001 From: Torkel Loman Date: Thu, 23 May 2024 15:16:45 -0400 Subject: [PATCH 27/28] Update src/steady_state_stability.jl Co-authored-by: Sam Isaacson --- src/steady_state_stability.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/steady_state_stability.jl b/src/steady_state_stability.jl index 8f97986c51..d544a8c638 100644 --- a/src/steady_state_stability.jl +++ b/src/steady_state_stability.jl @@ -40,7 +40,7 @@ form of maps. - 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 +these. Furthermore, Catalyst uses a tolerance `tol = 10*sqrt(eps())` to determine whether a computed eigenvalue is far away enough from 0 to be reliably used. This selected threshold can be changed through the `tol` argument. ``` """ From b466c1075b9ff6d92a27b3d8a17b8db511523cc8 Mon Sep 17 00:00:00 2001 From: Torkel Date: Thu, 23 May 2024 15:40:48 -0400 Subject: [PATCH 28/28] up --- HISTORY.md | 3 ++- docs/src/api.md | 2 +- .../steady_state_stability_computation.md | 6 +++++ .../bifurcation_kit_extension.jl | 2 +- .../homotopy_continuation_extension.jl | 2 +- src/Catalyst.jl | 5 ++-- src/reactionsystem.jl | 23 +++++++--------- src/reactionsystem_conversions.jl | 2 +- src/steady_state_stability.jl | 2 +- test/miscellaneous_tests/api.jl | 26 ++++++++++++------- test/reactionsystem_core/reactionsystem.jl | 12 --------- test/runtests.jl | 2 +- 12 files changed, 43 insertions(+), 44 deletions(-) diff --git a/HISTORY.md b/HISTORY.md index 41b23916d5..c568dc8ddd 100644 --- a/HISTORY.md +++ b/HISTORY.md @@ -85,7 +85,7 @@ plot(bif_dia; xguide="k1", yguide="X") ``` - Automatically handles elimination of conservation laws for computing bifurcation diagrams. - Updated Bifurcation documentation with respect to this new feature. -- Added function `is_autonomous` to check if a `ReactionSystem` is autonomous. +- Added function `isautonomous` to check if a `ReactionSystem` is autonomous. - Added function `steady_state_stability` to compute stability for steady states. Example: ```julia # Creates model. @@ -98,6 +98,7 @@ p = [:p => 1.0, :d => 0.5] steady_state = [2.0] steady_state_stability(steady_state, rn, p) ``` +Here, `steady_state_stability` take an optional argument `tol = 10*sqrt(eps())`, which is used to determine whether a eigenvalue real part is reliably less that 0. ## Catalyst 13.5 - Added a CatalystHomotopyContinuationExtension extension, which exports the `hc_steady_state` function if HomotopyContinuation is exported. `hc_steady_state` finds the steady states of a reaction system using the homotopy continuation method. This feature is only available for julia versions 1.9+. Example: diff --git a/docs/src/api.md b/docs/src/api.md index 5e267791a3..741340d0c7 100644 --- a/docs/src/api.md +++ b/docs/src/api.md @@ -166,7 +166,7 @@ speciesmap paramsmap reactionparamsmap isspecies -is_autonomous +isautonomous Catalyst.isconstant Catalyst.isbc ``` 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 df791ac87f..4cbe598758 100644 --- a/docs/src/steady_state_functionality/steady_state_stability_computation.md +++ b/docs/src/steady_state_functionality/steady_state_stability_computation.md @@ -36,6 +36,12 @@ Next, we can apply `steady_state_stability` to each steady state yielding a vect [steady_state_stability(sstate, sa_loop, ps) for sstate in steady_states] ``` +Finally, as described above, Catalyst uses an optional argument, `tol`, to determine how strict to make the stability check. I.e. below we set the tolerance to `1e-6` (a larger value, that is stricter, than the default of `10*sqrt(eps())`) +```@example stability_1 +[steady_state_stability(sstate, sa_loop, ps; tol = 1e-6) for sstate in steady_states] +nothing# hide +``` + ## 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. diff --git a/ext/CatalystBifurcationKitExtension/bifurcation_kit_extension.jl b/ext/CatalystBifurcationKitExtension/bifurcation_kit_extension.jl index 3e2ea2effc..ef62ee3bc1 100644 --- a/ext/CatalystBifurcationKitExtension/bifurcation_kit_extension.jl +++ b/ext/CatalystBifurcationKitExtension/bifurcation_kit_extension.jl @@ -3,7 +3,7 @@ # 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...) - if !is_autonomous(rs) + if !isautonomous(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 diff --git a/ext/CatalystHomotopyContinuationExtension/homotopy_continuation_extension.jl b/ext/CatalystHomotopyContinuationExtension/homotopy_continuation_extension.jl index 926e7ff0eb..1dcc64d9ba 100644 --- a/ext/CatalystHomotopyContinuationExtension/homotopy_continuation_extension.jl +++ b/ext/CatalystHomotopyContinuationExtension/homotopy_continuation_extension.jl @@ -36,7 +36,7 @@ Notes: ``` """ function Catalyst.hc_steady_states(rs::ReactionSystem, ps; filter_negative=true, neg_thres=-1e-20, u0=[], kwargs...) - if !is_autonomous(rs) + if !isautonomous(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) diff --git a/src/Catalyst.jl b/src/Catalyst.jl index 33ee2b87fb..469003af26 100644 --- a/src/Catalyst.jl +++ b/src/Catalyst.jl @@ -18,7 +18,7 @@ using DynamicQuantities#, Unitful # Having Unitful here as well currently gives @reexport using ModelingToolkit using Symbolics - +using LinearAlgebra using RuntimeGeneratedFunctions RuntimeGeneratedFunctions.init(@__MODULE__) @@ -42,7 +42,6 @@ import Base: (==), hash, size, getindex, setindex, isless, Sort.defalg, length, import MacroTools, Graphs import Graphs: DiGraph, SimpleGraph, SimpleDiGraph, vertices, edges, add_vertices!, nv, ne import DataStructures: OrderedDict, OrderedSet -import LinearAlgebra.eigvals import Parameters: @with_kw_noshow import Symbolics: occursin, wrap @@ -103,7 +102,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 isautonomous export reactionrates export isequivalent export set_default_noise_scaling diff --git a/src/reactionsystem.jl b/src/reactionsystem.jl index 976ed1fd77..98a808fb5a 100644 --- a/src/reactionsystem.jl +++ b/src/reactionsystem.jl @@ -1204,7 +1204,7 @@ function dependants(rx, network) end """ - is_autonomous(rs::ReactionSystem) +isautonomous(rs::ReactionSystem) Checks if a system is autonomous (i.e. no rate or equation depend on the independent variable(s)). Example: @@ -1212,27 +1212,24 @@ Example: rs1 = @reaction_system (p,d), 0 <--> X end -is_autonomous(rs1) # Returns `true`. +isautonomous(rs1) # Returns `true`. rs2 = @reaction_system (p/t,d), 0 <--> X end -is_autonomous(rs2) # Returns `false`. +isautonomous(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))] - if isempty(dep_var_param_rxs) && isempty(dep_var_param_eqs) - dep_var_param = [] - else - dep_var_param = reduce(vcat,[dep_var_param_rxs; dep_var_param_eqs]) +function isautonomous(rs::ReactionSystem) + # Get all variables occurring in reactions and equations. + vars = Set() + for eq in equations(rs) + (eq isa Reaction) ? get_variables!(vars, eq.rate) : get_variables!(vars, eq) end # 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) + (get_iv(rs) in vars) && return false + any(in(vars), get_sivs(rs)) && return false return true end diff --git a/src/reactionsystem_conversions.jl b/src/reactionsystem_conversions.jl index b7c7e656f3..38a634ff77 100644 --- a/src/reactionsystem_conversions.jl +++ b/src/reactionsystem_conversions.jl @@ -521,7 +521,7 @@ function Base.convert(::Type{<:NonlinearSystem}, rs::ReactionSystem; name = name # Error checks. iscomplete(rs) || error(COMPLETENESS_ERROR) spatial_convert_err(rs::ReactionSystem, NonlinearSystem) - if !is_autonomous(rs) + if !isautonomous(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 diff --git a/src/steady_state_stability.jl b/src/steady_state_stability.jl index d544a8c638..f8512116e4 100644 --- a/src/steady_state_stability.jl +++ b/src/steady_state_stability.jl @@ -47,7 +47,7 @@ computed eigenvalue is far away enough from 0 to be reliably used. This selected function steady_state_stability(u::Vector, rs::ReactionSystem, ps; tol = 10*sqrt(eps(ss_val_type(u))), ss_jac = steady_state_jac(rs; u0 = u)) # Warning checks. - if !is_autonomous(rs) + if !isautonomous(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 diff --git a/test/miscellaneous_tests/api.jl b/test/miscellaneous_tests/api.jl index ccf0ce2d26..78627a8b2b 100644 --- a/test/miscellaneous_tests/api.jl +++ b/test/miscellaneous_tests/api.jl @@ -447,7 +447,7 @@ let @test_throws MethodError Catalyst.to_multivariate_poly(neweqs) end -# Tests `is_autonomous` function. +# Tests `isautonomous` function. let # Using default iv. rn1 = @reaction_network begin @@ -462,9 +462,9 @@ 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) + @test !isautonomous(rn1) + @test !isautonomous(rn2) + @test isautonomous(rn3) # Using non-default iv. rn4 = @reaction_network begin @@ -487,15 +487,23 @@ 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) + @test !isautonomous(rn4) + @test !isautonomous(rn5) + @test !isautonomous(rn6) + @test isautonomous(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) + @test !isautonomous(rn7) + + # Using a registered function. + f(d,t) = d/(1 + t) + Symbolics.@register_symbolic f(d,t) + rn8 = @reaction_network begin + f(d,t), X --> 0 + end + @test !isautonomous(rn8) end \ No newline at end of file diff --git a/test/reactionsystem_core/reactionsystem.jl b/test/reactionsystem_core/reactionsystem.jl index 24a0d0ad70..7e854efcc7 100644 --- a/test/reactionsystem_core/reactionsystem.jl +++ b/test/reactionsystem_core/reactionsystem.jl @@ -155,18 +155,6 @@ let @test norm(du - du2) < 100 * eps() G2 = sf.g(sprob.u0, sprob.p, 0.0) @test norm(G - G2) < 100 * eps() - - # Test conversion to NonlinearSystem. Permanently broken as we no longer permit non-autonomous - # `ReactionSystem`s to be converted to `NonlinearSystem`s. - @test_broken false - # ns = convert(NonlinearSystem, rs) - # nlprob = NonlinearProblem(rs, u, p) - # fnl = eval(generate_function(ns)[2]) - # dunl = similar(du) - # @test_broken let # The next line throws an error. - # fnl(dunl, nlprob.u0, nlprob.p) - # @test norm(du - dunl) < 100 * eps() - # end end # Test with JumpSystem. diff --git a/test/runtests.jl b/test/runtests.jl index d7953daf40..8f255921d7 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -33,10 +33,10 @@ using SafeTestsets, Test #if GROUP == "All" || GROUP == "Miscellaneous-NetworkAnalysis" # Tests various miscellaneous features. @time @safetestset "API" begin include("miscellaneous_tests/api.jl") end + @time @safetestset "Units" begin include("miscellaneous_tests/units.jl") end @time @safetestset "Steady State Stability Computations" begin include("miscellaneous_tests/stability_computation.jl") end @time @safetestset "Compound Species" begin include("miscellaneous_tests/compound_macro.jl") end @time @safetestset "Reaction Balancing" begin include("miscellaneous_tests/reaction_balancing.jl") end - @time @safetestset "Units" begin include("miscellaneous_tests/units.jl") end # Tests reaction network analysis features. @time @safetestset "Conservation Laws" begin include("network_analysis/conservation_laws.jl") end