init

Author: TorkelE

ext/CatalystHomotopyContinuationExtension/homotopy_continuation_extension.jl

@@ -43,7 +43,7 @@ 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)
+    ns = convert(NonlinearSystem, rs; remove_conserved = true, expand_functions = 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))
@@ -63,27 +63,6 @@ function conservationlaw_errorcheck(rs, pre_varmap)
         error("The system has conservation laws but initial conditions were not provided for some species.")
 end
 
-# 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
-
 # Parses and expression and return a version where any exponents that are Float64 (but an int, like 2.0) are turned into Int64s.
 make_int_exps(expr) = wrap(Rewriters.Postwalk(Rewriters.PassThrough(___make_int_exps))(unwrap(expr))).val
 function ___make_int_exps(expr)
@@ -124,4 +103,27 @@ 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
+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

src/reactionsystem.jl

@@ -1287,14 +1287,19 @@ 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,
-                      kwargs...)
+                      expand_functions = false, 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 && for eq in eqs
+        eq.rhs = expand_registered_functions!(eq.rhs)
+    end
 
     ODESystem(eqs, get_iv(fullrs), sts, ps;
               observed = obs,
@@ -1326,7 +1331,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,
-                      kwargs...)
+                      expand_functions = false, kwargs...)
     spatial_convert_err(rs::ReactionSystem, NonlinearSystem)
     fullrs = Catalyst.flatten(rs)
     remove_conserved && conservationlaws(fullrs)
@@ -1335,6 +1340,12 @@ function Base.convert(::Type{<:NonlinearSystem}, rs::ReactionSystem; name = name
                          as_odes = false, include_zero_odes)
     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 && for eq in eqs
+        eq.rhs = expand_registered_functions!(eq.rhs)
+    end
+
     NonlinearSystem(eqs, sts, ps;
                     name,
                     observed = obs,
@@ -1374,7 +1385,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,
-                      kwargs...)
+                      expand_functions = false, kwargs...)
     spatial_convert_err(rs::ReactionSystem, SDESystem)
 
     flatrs = Catalyst.flatten(rs)
@@ -1400,6 +1411,11 @@ 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).
+    expand_functions && for eq in eqs
+        eq.rhs = expand_registered_functions!(eq.rhs)
+    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

src/registered_functions.jl

@@ -107,3 +107,28 @@ function Symbolics.derivative(::typeof(hillar), args::NTuple{5, Any}, ::Val{5})
      (args[4]^args[5]) * log(args[4])) /
     (args[1]^args[5] + args[2]^args[5] + args[4]^args[5])^2
 end
+
+"""
+expand_registered_functions!(expr)
+
+Takes an expression, and expands registered function expressions. E.g. `mm(X,v,K)` is replaced with v*X/(X+K). Currently supported functions: `mm`, `mmr`, `hill`, `hillr`, and `hill`.
+"""
+function expand_registered_functions!(expr)
+    istree(expr) || return expr
+    args = arguments(expr)
+    if operation(expr) == Catalyst.mm
+        return args[2]*args[1]/(args[1] + args[3])
+    elseif operation(expr) == Catalyst.mmr
+        return args[2]*args[3]/(args[1] + args[3])
+    elseif operation(expr) == Catalyst.hill
+        return args[2]*(args[1]^args[4])/((args[1])^args[4] + (args[3])^args[4])
+    elseif operation(expr) == Catalyst.hillr
+        return args[2]*(args[3]^args[4])/((args[1])^args[4] + (args[3])^args[4])
+    elseif operation(expr) == Catalyst.hillar
+        return args[3]*(args[1]^args[5])/((args[1])^args[5] + (args[2])^args[5] + (args[4])^args[5])
+    end
+    for i = 1:length(args)
+        args[i] = expand_registered_functions!(args[i])
+    end
+    return expr
+end