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