up

Author: TorkelE

ext/CatalystHomotopyContinuationExtension/homotopy_continuation_extension.jl

@@ -32,6 +32,7 @@ gives
 ```
 
 Notes:
+- Homotopy-based steady state finding only works when all rates are rational polynomials (e.g. constant, linear, mm, or hill functions).
 ```
   """
 function Catalyst.hc_steady_states(rs::ReactionSystem, ps; filter_negative=true, neg_thres=-1e-20, u0=[], kwargs...)
@@ -42,8 +43,9 @@ function Catalyst.hc_steady_states(rs::ReactionSystem, ps; filter_negative=true,
 end
 
 # 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, expand_functions = true)
+function steady_state_polynomial(rs_in::ReactionSystem, ps, u0)
+    rs = ModelingToolkit.@set rs_in.rxs = [expand_registered_functions(rx) for rx in rn.rxs]
+    ns = convert(NonlinearSystem, rs; remove_conserved = true)
     pre_varmap = [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))
@@ -102,27 +104,4 @@ function filter_negative_f(sols; neg_thres=-1e-20)
         (neg_thres < sol[idx] < 0) && (sol[idx] = 0)
     end
     return filter(sol -> all(>=(0), sol), sols)
-end
-
-### Archived ###
-
-# # Unfolds a function (like mm or hill). 
-# function deregister(fs::Vector{T}, expr) where T
-#     for f in fs
-#         expr = deregister(f, expr) 
-#     end
-#     return expr
-# end
-# # Provided by Shashi Gowda.
-# deregister(f, expr) = wrap(Rewriters.Postwalk(Rewriters.PassThrough(___deregister(f)))(unwrap(expr)))
-# function ___deregister(f)
-#     (expr) ->
-#     if istree(expr) && operation(expr) == f
-#         args = arguments(expr)
-#         invoke_with = map(args) do a
-#             t = typeof(a)
-#             issym(a) || istree(a) ? wrap(a) => symtype(a) : a => typeof(a)
-#         end 
-#         invoke(f, Tuple{last.(invoke_with)...}, first.(invoke_with)...)
-#     end
-# end
+end

src/reactionsystem.jl

@@ -538,11 +538,6 @@ struct ReactionSystem{V <: NetworkProperties} <:
         checks && validate(rs)
         rs
     end
-
-    # Copies a reaction system, but with the option of having some fields replaced
-    function ReactionSystem(rs::ReactionSystem; eqs = rs.eqs, rxs = rs.rxs, iv = rs.iv, sivs = rs.sivs, states = rs.states, species = rs.species, ps = rs.ps, var_to_name = rs.var_to_name, observed = rs.observed, name = rs.name, systems = rs.systems, defaults = rs.defaults, connection_type = rs.connection_type, networkproperties = rs.networkproperties, combinatoric_ratelaws = rs.combinatoric_ratelaws, continuous_events = rs.continuous_events, discrete_events = rs.discrete_events, complete = rs.complete)
-        new{typeof(networkproperties)}(eqs, rxs, ModelingToolkit.unwrap(iv), ModelingToolkit.unwrap.(sivs), ModelingToolkit.unwrap.(states), ModelingToolkit.unwrap.(species), ModelingToolkit.unwrap.(ps), var_to_name, observed, name, systems, defaults, connection_type, networkproperties, combinatoric_ratelaws, continuous_events, discrete_events, complete)
-    end
 end
 
 function get_speciestype(iv, states, systems)
@@ -1292,17 +1287,14 @@ Keyword args and default values:
 function Base.convert(::Type{<:ODESystem}, rs::ReactionSystem; name = nameof(rs),
                       combinatoric_ratelaws = get_combinatoric_ratelaws(rs),
                       include_zero_odes = true, remove_conserved = false, checks = false,
-                      expand_functions = false, kwargs...)
+                      kwargs...)
     spatial_convert_err(rs::ReactionSystem, ODESystem)
     fullrs = Catalyst.flatten(rs)
     remove_conserved && conservationlaws(fullrs)
     ists, ispcs = get_indep_sts(fullrs, remove_conserved)
     eqs = assemble_drift(fullrs, ispcs; combinatoric_ratelaws, remove_conserved,
                          include_zero_odes)
     eqs, sts, ps, obs, defs = addconstraints!(eqs, fullrs, ists, ispcs; remove_conserved)
-    
-    # Converts expressions like mm(X,v,K) to v*X/(X+K).
-    expand_functions  && (eqs = [eq.lhs ~ expand_registered_functions!(eq.rhs) for eq in eqs])
 
     ODESystem(eqs, get_iv(fullrs), sts, ps;
               observed = obs,
@@ -1334,7 +1326,7 @@ Keyword args and default values:
 function Base.convert(::Type{<:NonlinearSystem}, rs::ReactionSystem; name = nameof(rs),
                       combinatoric_ratelaws = get_combinatoric_ratelaws(rs),
                       include_zero_odes = true, remove_conserved = false, checks = false,
-                      expand_functions = false, kwargs...)
+                      kwargs...)
     spatial_convert_err(rs::ReactionSystem, NonlinearSystem)
     fullrs = Catalyst.flatten(rs)
     remove_conserved && conservationlaws(fullrs)
@@ -1344,9 +1336,6 @@ function Base.convert(::Type{<:NonlinearSystem}, rs::ReactionSystem; name = name
     error_if_constraint_odes(NonlinearSystem, fullrs)
     eqs, sts, ps, obs, defs = addconstraints!(eqs, fullrs, ists, ispcs; remove_conserved)
 
-    # Converts expressions like mm(X,v,K) to v*X/(X+K).
-    expand_functions  && (eqs = [eq.lhs ~ expand_registered_functions!(eq.rhs) for eq in eqs])
-
     NonlinearSystem(eqs, sts, ps;
                     name,
                     observed = obs,
@@ -1386,7 +1375,7 @@ function Base.convert(::Type{<:SDESystem}, rs::ReactionSystem;
                       noise_scaling = nothing, name = nameof(rs),
                       combinatoric_ratelaws = get_combinatoric_ratelaws(rs),
                       include_zero_odes = true, checks = false, remove_conserved = false,
-                      expand_functions = false, kwargs...)
+                      kwargs...)
     spatial_convert_err(rs::ReactionSystem, SDESystem)
 
     flatrs = Catalyst.flatten(rs)
@@ -1412,12 +1401,6 @@ function Base.convert(::Type{<:SDESystem}, rs::ReactionSystem;
     eqs, sts, ps, obs, defs = addconstraints!(eqs, flatrs, ists, ispcs; remove_conserved)
     ps = (noise_scaling === nothing) ? ps : vcat(ps, toparam(noise_scaling))
 
-    # Converts expressions like mm(X,v,K) to v*X/(X+K).
-    if expand_functions
-        eqs = [eq.lhs ~ expand_registered_functions!(eq.rhs) for eq in eqs]
-        noiseeqs = [expand_registered_functions!(neq) for neq in noiseeqs]
-    end
-
     if any(isbc, get_states(flatrs))
         @info "Boundary condition species detected. As constraint equations are not currently supported when converting to SDESystems, the resulting system will be undetermined. Consider using constant species instead."
     end