From 4f8be2975b664e5d9ca56641c0b4e0a32aed644f Mon Sep 17 00:00:00 2001 From: Sam Isaacson Date: Mon, 13 May 2024 11:18:37 -0400 Subject: [PATCH 1/2] add combinatoric_ratelaws to DSL --- src/reaction_network.jl | 16 ++++++++-- test/dsl/dsl_options.jl | 66 +++++++++++++++++++++++++++++++---------- 2 files changed, 63 insertions(+), 19 deletions(-) diff --git a/src/reaction_network.jl b/src/reaction_network.jl index df66043ca2..a8ee74b718 100644 --- a/src/reaction_network.jl +++ b/src/reaction_network.jl @@ -82,8 +82,9 @@ const forbidden_variables_error = let end # Declares the keys used for various options. -const option_keys = (:species, :parameters, :variables, :ivs, :compounds, :observables, :default_noise_scaling, - :differentials, :equations, :continuous_events, :discrete_events) +const option_keys = (:species, :parameters, :variables, :ivs, :compounds, :observables, + :default_noise_scaling, :differentials, :equations, + :continuous_events, :discrete_events, :combinatoric_ratelaws) ### The @species macro, basically a copy of the @variables macro. ### macro species(ex...) @@ -364,6 +365,13 @@ function make_reaction_system(ex::Expr; name = :(gensym(:ReactionSystem))) sivs = (length(ivs) > 1) ? Expr(:vect, ivs[2:end]...) : nothing all_ivs = (isnothing(sivs) ? [tiv] : [tiv; sivs.args]) + if haskey(options, :combinatoric_ratelaws) + combinatoric_ratelaws = options[:combinatoric_ratelaws].args[end] + else + combinatoric_ratelaws = true + end + @show combinatoric_ratelaws + # Reads more options. observed_vars, observed_eqs, obs_syms = read_observed_options(options, [species_declared; variables], all_ivs) @@ -415,7 +423,9 @@ function make_reaction_system(ex::Expr; name = :(gensym(:ReactionSystem))) Catalyst.make_ReactionSystem_internal( $rxexprs, $tiv, setdiff(union($spssym, $varssym, $compssym), $obs_syms), $pssym; name = $name, spatial_ivs = $sivs, observed = $observed_eqs, - continuous_events = $continuous_events_expr, discrete_events = $discrete_events_expr); + continuous_events = $continuous_events_expr, + discrete_events = $discrete_events_expr, + combinatoric_ratelaws = $combinatoric_ratelaws); default_reaction_metadata = $default_reaction_metadata ) end diff --git a/test/dsl/dsl_options.jl b/test/dsl/dsl_options.jl index de798e46a2..a9a6026317 100644 --- a/test/dsl/dsl_options.jl +++ b/test/dsl/dsl_options.jl @@ -398,7 +398,7 @@ end # Test basic functionality. # Tests various types of indexing. -let +let rn = @reaction_network begin @observables begin X ~ Xi + Xa @@ -489,23 +489,23 @@ end # Tests using a single observable (without begin/end statement). # Tests using observable component not part of reaction. # Tests using parameters in observables formula. -let +let rn = @reaction_network begin @parameters op_1 op_2 @species X4(t) - @observables X ~ X1^2 + op_1*(X2 + 2X3) + X1*X4/op_2 + p + @observables X ~ X1^2 + op_1*(X2 + 2X3) + X1*X4/op_2 + p (p,d), 0 <--> X1 (k1,k2), X1 <--> X2 (k3,k4), X2 <--> X3 end - + u0 = Dict([:X1 => 1.0, :X2 => 2.0, :X3 => 3.0, :X4 => 4.0]) ps = Dict([:p => 1.0, :d => 0.2, :k1 => 1.5, :k2 => 1.5, :k3 => 5.0, :k4 => 5.0, :op_1 => 1.5, :op_2 => 1.5]) oprob = ODEProblem(rn, u0, (0.0, 1000.0), ps) sol = solve(oprob, Tsit5()) - @test sol[:X][1] == u0[:X1]^2 + ps[:op_1]*(u0[:X2] + 2*u0[:X3]) + u0[:X1]*u0[:X4]/ps[:op_2] + ps[:p] + @test sol[:X][1] == u0[:X1]^2 + ps[:op_1]*(u0[:X2] + 2*u0[:X3]) + u0[:X1]*u0[:X4]/ps[:op_2] + ps[:p] end # Checks that ivs are correctly found. @@ -535,7 +535,7 @@ end # Declares observables implicitly/explicitly. # Cannot test `isequal(rn1, rn2)` because the two sets of observables have some obscure Symbolics # substructure that is different. -let +let # Basic case. rn1 = @reaction_network rn_observed begin @observables X ~ X1 + X2 @@ -609,14 +609,14 @@ let (kB,kD), 2X <--> X2 (k1,k2), Y1 <--> Y2 end - + @test isspecies(rn.X) @test !isspecies(rn.Y) @test !isspecies(rn.Z) end # Tests various erroneous declarations throw errors. -let +let # Independent variable in @observables. @test_throws Exception @eval @reaction_network begin @observables X(t) ~ X1 + X2 @@ -699,7 +699,7 @@ let @equations begin X + 5 ~ k*S 3Y + X ~ S + X*d - end + end (p,d), 0 <--> S end @@ -747,7 +747,7 @@ let rn2 = @reaction_network rn begin @parameters k @variables X(t) - @equations begin + @equations begin X + 2 ~ k*S end (p,d), 0 <--> S @@ -759,7 +759,7 @@ end # Tries with interpolating a value into an equation. # Tries using rn.X notation for designating variables. # Tries for empty parameter vector. -let +let c = 6.0 rn = complete(@reaction_network begin @variables X(t) @@ -774,12 +774,12 @@ let end # Checks hierarchical model. -let +let base_rn = @network_component begin @variables V1(t) @equations begin X*3V1 ~ X - 2 - end + end (p,d), 0 <--> X end @unpack X, V1, p, d = base_rn @@ -788,7 +788,7 @@ let @variables V2(t) @equations begin X*4V2 ~ X - 3 - end + end (p,d), 0 <--> X end @@ -813,7 +813,7 @@ let @equations begin X + 5 ~ k*S D(Y) ~ X + S - 5*Y - end + end (p,d), 0 <--> S end @unpack X, Y, S, p, d, k = rn @@ -848,7 +848,7 @@ let end # Tests that various erroneous declarations throw errors. -let +let # Using = instead of ~ (for equation). @test_throws Exception @eval @reaction_network begin @variables X(t) @@ -861,4 +861,38 @@ let @equations X ~ p - S (P,D), 0 <--> S end +end + +# test combinatoric_ratelaws DSL option +let + rn = @reaction_network begin + @combinatoric_ratelaws false + (k1,k2), 2A <--> B + end + combinatoric_ratelaw = Catalyst.get_combinatoric_ratelaws(rn) + @test combinatoric_ratelaw == false + rl = oderatelaw(reactions(rn)[1]; combinatoric_ratelaw) + @unpack k1, A = rn + @test isequal(rl, k1*A^2) + + rn2 = @reaction_network begin + @combinatoric_ratelaws true + (k1,k2), 2A <--> B + end + combinatoric_ratelaw = Catalyst.get_combinatoric_ratelaws(rn2) + @test combinatoric_ratelaw == true + rl = oderatelaw(reactions(rn2)[1]; combinatoric_ratelaw) + @unpack k1, A = rn2 + @test isequal(rl, k1*A^2/2) + + crl = false + rn3 = @reaction_network begin + @combinatoric_ratelaws $crl + (k1,k2), 2A <--> B + end + combinatoric_ratelaw = Catalyst.get_combinatoric_ratelaws(rn3) + @test combinatoric_ratelaw == crl + rl = oderatelaw(reactions(rn3)[1]; combinatoric_ratelaw) + @unpack k1, A = rn3 + @test isequal(rl, k1*A^2) end \ No newline at end of file From a04b5274a7dcd521abd2e82e231392cb199caa3f Mon Sep 17 00:00:00 2001 From: Sam Isaacson Date: Mon, 13 May 2024 11:24:27 -0400 Subject: [PATCH 2/2] remove show --- src/reaction_network.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/src/reaction_network.jl b/src/reaction_network.jl index a8ee74b718..5cafae7a91 100644 --- a/src/reaction_network.jl +++ b/src/reaction_network.jl @@ -370,7 +370,6 @@ function make_reaction_system(ex::Expr; name = :(gensym(:ReactionSystem))) else combinatoric_ratelaws = true end - @show combinatoric_ratelaws # Reads more options. observed_vars, observed_eqs, obs_syms = read_observed_options(options, [species_declared; variables], all_ivs)