diff --git a/HISTORY.md b/HISTORY.md index 3b51630da3..c568dc8ddd 100644 --- a/HISTORY.md +++ b/HISTORY.md @@ -85,6 +85,20 @@ 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 `isautonomous` 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) +``` +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/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/api.md b/docs/src/api.md index 51a403c1a5..741340d0c7 100644 --- a/docs/src/api.md +++ b/docs/src/api.md @@ -166,6 +166,7 @@ speciesmap paramsmap reactionparamsmap isspecies +isautonomous Catalyst.isconstant Catalyst.isbc ``` 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 diff --git a/docs/src/steady_state_functionality/steady_state_stability_computation.md b/docs/src/steady_state_functionality/steady_state_stability_computation.md new file mode 100644 index 0000000000..4cbe598758 --- /dev/null +++ b/docs/src/steady_state_functionality/steady_state_stability_computation.md @@ -0,0 +1,63 @@ +# 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 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: +```@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](@ref ref): +```@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` to each steady state yielding a vector of their stabilities: +```@example stability_1 +[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. + +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) + +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 +``` + +!!! 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/ext/CatalystBifurcationKitExtension/bifurcation_kit_extension.jl b/ext/CatalystBifurcationKitExtension/bifurcation_kit_extension.jl index f99389c6df..ef62ee3bc1 100644 --- a/ext/CatalystBifurcationKitExtension/bifurcation_kit_extension.jl +++ b/ext/CatalystBifurcationKitExtension/bifurcation_kit_extension.jl @@ -3,6 +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...) + 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 # 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..1dcc64d9ba 100644 --- a/ext/CatalystHomotopyContinuationExtension/homotopy_continuation_extension.jl +++ b/ext/CatalystHomotopyContinuationExtension/homotopy_continuation_extension.jl @@ -36,6 +36,9 @@ Notes: ``` """ function Catalyst.hc_steady_states(rs::ReactionSystem, ps; filter_negative=true, neg_thres=-1e-20, u0=[], kwargs...) + 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) sols = HC.real_solutions(HC.solve(ss_poly; kwargs...)) reorder_sols!(sols, ss_poly, rs) @@ -57,14 +60,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) diff --git a/src/Catalyst.jl b/src/Catalyst.jl index b347f0ab31..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__) @@ -102,6 +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 isautonomous export reactionrates export isequivalent export set_default_noise_scaling @@ -149,6 +150,10 @@ export @compound, @compounds export iscompound, components, coefficients, component_coefficients export balance_reaction +# Functionality for computing the stability of system steady states. +include("steady_state_stability.jl") +export steady_state_stability, steady_state_jac + ### Extensions ### # HomotopyContinuation diff --git a/src/reactionsystem.jl b/src/reactionsystem.jl index 47a0af56b4..98a808fb5a 100644 --- a/src/reactionsystem.jl +++ b/src/reactionsystem.jl @@ -1203,6 +1203,36 @@ function dependants(rx, network) dependents(rx, network) end +""" +isautonomous(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 +isautonomous(rs1) # Returns `true`. + +rs2 = @reaction_system + (p/t,d), 0 <--> X +end +isautonomous(rs2) # Returns `false`. +``` +""" +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 + (get_iv(rs) in vars) && return false + any(in(vars), get_sivs(rs)) && return false + return true +end + ### `ReactionSystem` Remaking ### diff --git a/src/reactionsystem_conversions.jl b/src/reactionsystem_conversions.jl index f122d3f62d..38a634ff77 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 !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 + + # 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 new file mode 100644 index 0000000000..f8512116e4 --- /dev/null +++ b/src/steady_state_stability.jl @@ -0,0 +1,134 @@ +### Stability Analysis ### + +""" + stability(u::Vector{T}, rs::ReactionSystem, p; tol = 10*sqrt(eps()) + ss_jac = steady_state_jac(u, rs, p)) + +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. +- `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. +- `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. + +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 its stability. +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. +- 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 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))), + ss_jac = steady_state_jac(rs; u0 = u)) + # Warning checks. + 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 + + # 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 + + # 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) + + # 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, but note that floating point error in the eigenvalue calculation could lead to incorrect results.") + 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 = []) + +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. + +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) + +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 = [sp => 0.0 for sp in unknowns(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. + u0 = steady_state_u_conversion(u0, rs) + conservationlaw_errorcheck(rs, u0) + + # Creates an `ODEProblem` with a Jacobian. Dummy values for `u0` and `ps` must be provided. + ps = [p => 0.0 for p in parameters(rs)] + return ODEProblem(rs, u0, 0, ps; jac = true, remove_conserved = true, + combinatoric_ratelaws = combinatoric_ratelaws) +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 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/extensions/homotopy_continuation.jl b/test/extensions/homotopy_continuation.jl index bc0f30e3a4..7a260e3b59 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 diff --git a/test/miscellaneous_tests/api.jl b/test/miscellaneous_tests/api.jl index 2c42475628..78627a8b2b 100644 --- a/test/miscellaneous_tests/api.jl +++ b/test/miscellaneous_tests/api.jl @@ -446,3 +446,64 @@ let neweqs = getfield.(equations(ns), :rhs) @test_throws MethodError Catalyst.to_multivariate_poly(neweqs) end + +# Tests `isautonomous` 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 !isautonomous(rn1) + @test !isautonomous(rn2) + @test isautonomous(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 !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 !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/miscellaneous_tests/stability_computation.jl b/test/miscellaneous_tests/stability_computation.jl new file mode 100644 index 0000000000..36f6125f02 --- /dev/null +++ b/test/miscellaneous_tests/stability_computation.jl @@ -0,0 +1,98 @@ +### Fetch Packages ### + +# Fetch packages. +using Catalyst, OrdinaryDiffEq, SteadyStateDiffEq, Test +import HomotopyContinuation + +# Sets rnd number. +using StableRNGs +rng = StableRNG(12345) + +### Basic Tests ### + +# Tests that stability is correctly assessed (using simulation) in multi stable system. +# Tests that `steady_state_jac` function works. +# 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) + + # Repeats several times, most steady state stability cases should be encountered several times. + 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)) + + # 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; 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_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_4 = steady_state_stability(Pair.(unknowns(rn), ss), rn, ps) + @test stabs_3 == stabs_4 + 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 + # Creates model. + 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 + + # 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] + 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] + + # 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] + 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]) +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 diff --git a/test/reactionsystem_core/reactionsystem.jl b/test/reactionsystem_core/reactionsystem.jl index 2bdb3df42e..7e854efcc7 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] @@ -157,16 +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. - 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 744ade8cdb..8f255921d7 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -33,9 +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