From 4afd0a767a7c75cbe39abff101375bd4cc326d3e Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Wed, 11 Dec 2024 12:01:33 +0530 Subject: [PATCH 1/2] fix: handle bug in `remove_denominators` --- ext/CatalystHomotopyContinuationExtension.jl | 2 +- .../homotopy_continuation_extension.jl | 41 ++++++++++++++++++- 2 files changed, 41 insertions(+), 2 deletions(-) diff --git a/ext/CatalystHomotopyContinuationExtension.jl b/ext/CatalystHomotopyContinuationExtension.jl index a4785992fe..cf376db84a 100644 --- a/ext/CatalystHomotopyContinuationExtension.jl +++ b/ext/CatalystHomotopyContinuationExtension.jl @@ -6,7 +6,7 @@ import DynamicPolynomials import ModelingToolkit as MT import HomotopyContinuation as HC import Setfield: @set -import Symbolics: unwrap, wrap, Rewriters, symtype, issym +import Symbolics: unwrap, wrap, Rewriters, symtype, issym, maketerm, BasicSymbolic, metadata using Symbolics: iscall # Creates and exports hc_steady_states function. diff --git a/ext/CatalystHomotopyContinuationExtension/homotopy_continuation_extension.jl b/ext/CatalystHomotopyContinuationExtension/homotopy_continuation_extension.jl index b0093e5a03..2de2603c69 100644 --- a/ext/CatalystHomotopyContinuationExtension/homotopy_continuation_extension.jl +++ b/ext/CatalystHomotopyContinuationExtension/homotopy_continuation_extension.jl @@ -59,7 +59,7 @@ function steady_state_polynomial(rs::ReactionSystem, ps, u0) eqs_pars_funcs = vcat(equations(ns), conservedequations(rs)) eqs = map(eq -> substitute(eq.rhs - eq.lhs, p_dict), eqs_pars_funcs) eqs_intexp = make_int_exps.(eqs) - ss_poly = Catalyst.to_multivariate_poly(remove_denominators.(eqs_intexp)) + ss_poly = Catalyst.to_multivariate_poly(remove_denominators.(common_denominator.(eqs_intexp))) return poly_type_convert(ss_poly) end @@ -78,6 +78,45 @@ function ___make_int_exps(expr) end end +function common_denominator(expr) + iscall(expr) || return expr + if operation(expr) == / + num, den = arguments(expr) + num = common_denominator(num) + den = common_denominator(den) + return num / den + end + if operation(expr) == + + num = 0 + den = 1 + for arg in arguments(expr) + arg = common_denominator(arg) + if iscall(arg) && operation(arg) == / + n, d = arguments(arg) + else + n = arg + d = 1 + end + num = num * d + den * n + den *= d + end + return num / den + end + if operation(expr) == ^ + base, pow = arguments(expr) + base = common_denominator(base) + if iscall(base) && operation(base) == / + num, den = arguments(base) + else + num, den = base, 1 + end + num ^= pow + den ^= pow + return num / den + end + return maketerm(BasicSymbolic, operation(expr), map(common_denominator, arguments(expr)), metadata(expr)) +end + # If the input is a fraction, removes the denominator. function remove_denominators(expr) s_expr = simplify_fractions(expr) From 3d08d21cb54a65d73f4fd45b982f874128d1076b Mon Sep 17 00:00:00 2001 From: Torkel Loman Date: Mon, 16 Dec 2024 16:20:19 +0000 Subject: [PATCH 2/2] add test --- .../homotopy_continuation_extension.jl | 21 +++++--- test/extensions/homotopy_continuation.jl | 48 +++++++++++++------ 2 files changed, 49 insertions(+), 20 deletions(-) diff --git a/ext/CatalystHomotopyContinuationExtension/homotopy_continuation_extension.jl b/ext/CatalystHomotopyContinuationExtension/homotopy_continuation_extension.jl index 2de2603c69..ec980615a2 100644 --- a/ext/CatalystHomotopyContinuationExtension/homotopy_continuation_extension.jl +++ b/ext/CatalystHomotopyContinuationExtension/homotopy_continuation_extension.jl @@ -10,8 +10,8 @@ Arguments: - `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. 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. +- `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 ```@repl @@ -48,6 +48,8 @@ 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) + # Creates the appropriate nonlinear system, and converts parameters to a form that can + # be substituted in later. rs = Catalyst.expand_registered_functions(rs) ns = complete(convert(NonlinearSystem, rs; remove_conserved = true, remove_conserved_warn = false)) @@ -56,10 +58,15 @@ function steady_state_polynomial(rs::ReactionSystem, ps, u0) p_vals = ModelingToolkit.varmap_to_vars(pre_varmap, parameters(ns); defaults = ModelingToolkit.defaults(ns)) p_dict = Dict(parameters(ns) .=> p_vals) - eqs_pars_funcs = vcat(equations(ns), conservedequations(rs)) - eqs = map(eq -> substitute(eq.rhs - eq.lhs, p_dict), eqs_pars_funcs) - eqs_intexp = make_int_exps.(eqs) - ss_poly = Catalyst.to_multivariate_poly(remove_denominators.(common_denominator.(eqs_intexp))) + + # Step by step convert the equation to something HC can work on (adds conserved equations, + # inserts parameter values and put everything on a side, converts e.g. 2.0 to 2, remove + # denominators to avoid rational polynomials). + eqs = vcat(equations(ns), conservedequations(rs)) + eqs = [substitute(eq.rhs - eq.lhs, p_dict) for eq in eqs] + eqs = make_int_exps.(eqs) + eqs = [remove_denominators(common_denominator(eq)) for eq in eqs] + ss_poly = Catalyst.to_multivariate_poly(eqs) return poly_type_convert(ss_poly) end @@ -78,6 +85,8 @@ function ___make_int_exps(expr) end end +# Converts an expression of the form p(x)/q(x) + r(x)/s(x) to t(x)/u(x) (i.e. puts everything +# above a single denominator, which is what `remove_denominators` is able to simplify away). function common_denominator(expr) iscall(expr) || return expr if operation(expr) == / diff --git a/test/extensions/homotopy_continuation.jl b/test/extensions/homotopy_continuation.jl index 1cb31a1bfc..409441b64e 100644 --- a/test/extensions/homotopy_continuation.jl +++ b/test/extensions/homotopy_continuation.jl @@ -14,7 +14,7 @@ include("../test_functions.jl") # 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 +let # Creates the model. rs = @reaction_network begin (k1,k2), X1 <--> X2 @@ -37,7 +37,7 @@ end # Tests for Symbol parameter input. # Tests that passing kwargs to HC.solve does not error and have an effect (i.e. modifying the seed # slightly modified the output in some way). -let +let wilhelm_2009_model = @reaction_network begin k1, Y --> 2X k2, 2X --> X + Y @@ -61,7 +61,7 @@ end # 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 + @parameters kX1=1.0 kX2=2.0 kY1=12345.0 @species X1(t)=0.1 X2(t)=0.2 Y1(t)=12345.0 (kX1,kX2), X1 <--> X2 (kY1,kY2), Y1 <--> Y2 @@ -69,36 +69,36 @@ let end 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, seed = 0x000004d1), by = sol->sol[1]) @test ss_1 ≈ [[0.2, 0.1, 3.0, 1.0, 40.0, 10.0]] - + rs_2 = @reaction_network begin - @parameters kX1=1.0 kX2=2.0 kY1=12345.0 + @parameters kX1=1.0 kX2=2.0 kY1=12345.0 @species C2(t)=0.1 C1(t)=0.2 B2(t)=12345.0 (kX1,kX2), C2 <--> C1 (kY1,kY2), B2 <--> B1 (kZ1,kZ2), A2 <--> A1 end u0_2 = [:B2 => 1.0, :B1 => 3.0, :A2 => 10.0, :A1 =>40.0] - + ss_2 = sort(hc_steady_states(rs_2, ps; u0 = u0_2, show_progress = false, seed = 0x000004d1), by = sol->sol[1]) @test ss_1 ≈ ss_2 end # Tests that non-scalar reaction rates work. -# Tests that rational polynomial steady state systems work. +# Tests that rational polynomial steady state systems work. # Tests that Hill function is correctly expanded even if nested. # Test filter_negative=false works. # Tests than non-integer exponents throws an error. -let +let rs = @reaction_network begin v*(0.1/v + hill(X,1,K,n)), 0 --> X d, X --> 0 end ps = Dict([:v => 5.0, :K => 2.5, :n => 3, :d => 1.0]) sss = hc_steady_states(rs, ps; filter_negative = false, show_progress = false, seed = 0x000004d1) - + @test length(sss) == 4 for ss in sss @test ps[:v]*(0.1/ps[:v] + ss[1]^ps[:n]/(ss[1]^ps[:n] + ps[:K]^ps[:n])) - ps[:d]*ss[1]≈ 0.0 atol = 1e-12 @@ -108,6 +108,26 @@ let @test_throws Exception hc_steady_states(rs, ps; show_progress = false, seed = 0x000004d1) end +# Checks that the correct steady states are found for system with multiple, known, steady states. +# Checks where (activating/repressing) Hill function is written out/using `hillar`. +# The system is known to have (exactly five steady states for the given parameter set. +let + # Finds the model steady states. + rs = @reaction_network begin + 0.01 + hillar(X,Y,1.0,Kx,3), ∅ --> X + 0.01 + (Y^3) / (X^3 + Y^3 + Ky^3), ∅ --> Y + 1.0, (X, Y) --> ∅ + end + ps = [:Kx => 0.3, :Ky => 0.5] + sss = hc_steady_states(rs, ps) + + # Checks that the steady states are correct, all unique, and that 5 (known number) is found. + @test length(sss) == 5 + @test allunique(sss) + for ss in sss + @test f_eval(rs, ss, ps, 0.0) ≈ [0.0, 0.0] atol = 1e-12 rtol = 1e-12 + end +end ### Other Tests ### @@ -131,22 +151,22 @@ let end # Checks that `hc_steady_states` cannot be applied to non-complete `ReactionSystems`s. -let +let # Create model. incomplete_network = @network_component begin (p, d), 0 <--> X end p_start = [:p => 1.0, :d => 0.2] - + # Computes bifurcation diagram. @test_throws Exception hc_steady_states(incomplete_network, p_start; show_progress = false, seed = 0x000004d1) end # Tests that non-autonomous system throws an error -let +let rs = @reaction_network begin (k,t), 0 <--> X end ps = [:k => 1.0] @test_throws Exception hc_steady_states(rs, ps; show_progress = false, seed = 0x000004d1) -end \ No newline at end of file +end