From 168044059f44196ac28557a12edd37c1ce7a292a Mon Sep 17 00:00:00 2001 From: Torkel Date: Wed, 4 Oct 2023 09:37:39 -0400 Subject: [PATCH] update --- .../homotopy_continuation.md | 4 +- .../homotopy_continuation_extension.jl | 16 ++++---- src/miscellaneous_functions.jl | 38 ------------------- test/extensions/homotopy_continuation.jl | 11 +++--- 4 files changed, 16 insertions(+), 53 deletions(-) delete mode 100644 src/miscellaneous_functions.jl diff --git a/docs/src/catalyst_applications/homotopy_continuation.md b/docs/src/catalyst_applications/homotopy_continuation.md index 0ef5c44b36..bda625a9d8 100644 --- a/docs/src/catalyst_applications/homotopy_continuation.md +++ b/docs/src/catalyst_applications/homotopy_continuation.md @@ -60,12 +60,12 @@ hc_steady_states(wilhelm_2009_model, ps; u0) ## Final notes - `hc_steady_states` supports any systems where all rates are systems of rational polynomials (such as Hill functions with integer Hill coefficients). - +- When providing initial conditions to compute conservation laws, values are only required for those species that are part of conserved quantities. If this set of species is unknown, it is recommended to provide initial conditions for all species. - Additional arguments provided to `hc_steady_states` are automatically passed to HomotopyContinuation's `solve` command. Use e.g. `show_progress=false` to disable the progress bar. --- ## [Citation](@id homotopy_continuation_citation) -If you use this functionality in your reasearch, please cite the following paper to support the authors of the HomotopyContinuation package: +If you use this functionality in your research, please cite the following paper to support the authors of the HomotopyContinuation package: ``` @inproceedings{HomotopyContinuation.jl, title={{H}omotopy{C}ontinuation.jl: {A} {P}ackage for {H}omotopy {C}ontinuation in {J}ulia}, diff --git a/ext/CatalystHomotopyContinuationExtension/homotopy_continuation_extension.jl b/ext/CatalystHomotopyContinuationExtension/homotopy_continuation_extension.jl index 40f74debb5..be29fa28d1 100644 --- a/ext/CatalystHomotopyContinuationExtension/homotopy_continuation_extension.jl +++ b/ext/CatalystHomotopyContinuationExtension/homotopy_continuation_extension.jl @@ -9,8 +9,8 @@ Arguments: - `rs::ReactionSystem`: The reaction system for which we want to find the steady states. - `ps`: The parameter values for which we want to find the steady states. - `filter_negative=true`: If set to true, solutions with any species concentration neg_thres`` but `< 0.0` are set to `0.0`. -- `u0=nothing`: Initial conditions for which we want to find the steady states. For systems with conservation laws this are required to compute conserved quantities. +- `neg_thres=-1e-20`: Determine the minimum values for which a species concentration is to be considered non-negative. Species concentrations ``> neg_thres`` but `< 0.0` are set to `0.0`. +- `u0=nothing`: Initial conditions for which we want to find the steady states. For systems with conservation laws this are required to compute conserved quantities. Initial conditions are not required for all species, only those involved in conserved quantities (if this set is unknown, it is recommended to provide initial conditions for all species). - `kwargs...`: any additional arguments (like `show_progress= true`) are passed into HomotopyContinuation.jl's `solve` call. Examples @@ -41,10 +41,10 @@ function Catalyst.hc_steady_states(rs::ReactionSystem, ps; filter_negative=true, return (filter_negative ? filter_negative_f(sols; neg_thres) : sols) end -# For a given reaction system, paraemter values, and initial conditions, find the polynomial that HC solves to find steady states. +# For a given reaction system, parameter values, and initial conditions, find the polynomial that HC solves to find steady states. function steady_state_polynomial(rs::ReactionSystem, ps, u0) ns = convert(NonlinearSystem, rs; remove_conserved = true) - pre_varmap = map(pair -> pair, [symmap_to_varmap(rs,u0); symmap_to_varmap(rs,ps)]) + pre_varmap = map(pair -> pair, [symmap_to_varmap(rs,u0)..., symmap_to_varmap(rs,ps)...]) conservationlaw_errorcheck(rs, pre_varmap) p_vals = ModelingToolkit.varmap_to_vars(pre_varmap, parameters(ns); defaults = ModelingToolkit.defaults(ns)) p_dict = Dict(parameters(ns) .=> p_vals) @@ -58,9 +58,9 @@ 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) - isempty(intersect(species(rs), vars_with_vals)) || return + 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 all species.") + error("The system has conservation laws but initial conditions were not provided for some species.") end # Unfolds a function (like mm or hill). @@ -105,12 +105,12 @@ function remove_denominators(expr) return remove_denominators(arguments(s_expr)[1]) end if operation(s_expr) == + - return +(remove_denominators.(arguments(s_expr))...) + return sum(remove_denominators(arg) for arg in arguments(s_expr)) end return s_expr end -# HC orders the solution vector according to the lexiographic values of the variable names. This reorders the output acording to the species index in the reaction system species vector. +# HC orders the solution vector according to the lexicographic values of the variable names. This reorders the output according to the species index in the reaction system species vector. function reorder_sols!(sols, ss_poly, rs::ReactionSystem) var_names_extended = String.(Symbol.(HC.variables(ss_poly))) var_names = [Symbol(s[1:prevind(s, findlast('_', s))]) for s in var_names_extended] diff --git a/src/miscellaneous_functions.jl b/src/miscellaneous_functions.jl deleted file mode 100644 index 3675fe5fad..0000000000 --- a/src/miscellaneous_functions.jl +++ /dev/null @@ -1,38 +0,0 @@ -# Contains various functions that does not clearly belong to any other files. Hopefully these can eventually be moved to a better place. - -### Steady State Stability ### - -""" - steady_state_stability(s_state, ps, sys) - -Computes the stability of steady states. For `Vector{Float64}` input, returns a `Bool` corresponding the state's stability. For `Vector{Vector{Float64}}` input, returns a `Vector{Bool}` corresponding the states' stability. - -Arguments: -- `s_state`: The steate(s) for which we want to compute the stability. These are assumed to be system steady states. -- `ps` the parameter values for which we wish to compute stability. -- `sys`: The system for which we wish to comptue stability. Can be a `ReactionSystem` or an `ODEFunction`. - -Example -```@repl -rs = @reaction_network begin - (p,d), 0 <--> X -end -ps = [:p => 1.0, :d => 0.2] -steady_state_stability([5.0], ps, rs) -``` -gives -``` -true -``` - -Notes: -- If a `ReactionSystem` is given, it is coverted to a `ODEFunction` internally. If you plan to compute the stability for a large number of steady states. First covert your ReactionSystem to a ODEFunction using `ODEFunction(convert(ODESystem, rs); jac=true)` and then use this in your input. - - -""" -steady_state_stability(s_states::Vector{Float64}, ps, rs::ReactionSystem) = steady_state_stability(s_states, ODEFunction(convert(ODESystem, rs); jac=true)) -# If a ODESystem is given directyl, no conversion is needed. -function steady_state_stability(s_states::Vector{Float64}, ps, ofun::ODEFunction) - maximum(eigen(ofun.jac(ss, last.(ps), Inf)).values) < 0.0 -end -steady_state_stability(s_states::Vector{Vector{Float64}}, args...) = steady_state_stability.(s_states, args...) diff --git a/test/extensions/homotopy_continuation.jl b/test/extensions/homotopy_continuation.jl index 8174f6ec2e..144774db80 100644 --- a/test/extensions/homotopy_continuation.jl +++ b/test/extensions/homotopy_continuation.jl @@ -6,7 +6,7 @@ import HomotopyContinuation # Tests for network without conservation laws. # Tests for Symbol parameter input. -# Tests for Symbolics initial condiiton input. +# Tests for Symbolics initial condition input. # Tests for different types (Symbol/Symbolics) for parameters and initial conditions. # Tests that attempts to find steady states of system with conservation laws, while u0 is not provided, gives an error. let @@ -47,8 +47,9 @@ let end # Tests that reordering is correct. -# Tests corectness in presence of default values. -# Tests where some defaul values are overwritten with other values +# Tests correctness in presence of default values. +# Tests where some default values are overwritten with other values. +# Tests where input ps/u0 are tuples with mixed types. let rs_1 = @reaction_network begin @parameters kX1=1.0 kX2=2.0 kY1=12345.0 @@ -57,8 +58,8 @@ let (kY1,kY2), Y1 <--> Y2 (kZ1,kZ2), Z1 <--> Z2 end - ps = [:kY1 => 1.0, :kY2 => 3.0, :kZ1 => 1.0, :kZ2 => 4.0] - u0_1 = [:Y1 => 1.0, :Y2 => 3.0, :Z1 => 10.0, :Z2 =>40.0] + ps = (:kY1 => 1.0, :kY2 => 3, :kZ1 => 1.0, :kZ2 => 4.0) + u0_1 = (:Y1 => 1.0, :Y2 => 3, :Z1 => 10, :Z2 =>40.0) ss_1 = sort(hc_steady_states(rs_1, ps; u0=u0_1, show_progress=false), by=sol->sol[1]) @test ss_1 ≈ [[0.2, 0.1, 3.0, 1.0, 40.0, 10.0]]