up

Author: TorkelE

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:

docs/src/api.md

@@ -166,7 +166,7 @@ speciesmap
 paramsmap
 reactionparamsmap
 isspecies
-is_autonomous
+isautonomous
 Catalyst.isconstant
 Catalyst.isbc
 ```

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. 
 

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
 

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)

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

src/reactionsystem.jl

@@ -1204,35 +1204,32 @@ 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:
 ```julia
 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
 

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
 

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
 

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