update

Author: TorkelE

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},

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 is removed from the output.
-- `neg_thres=-1e-20`: Determine the minimum values for which a species concentration is to be considred non-negative. Species conentrations ``> 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]

src/miscellaneous_functions.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]]