From 8431dfafd026a9b96ee3e0f1b4ff2bce4b5f973e Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Tue, 21 Mar 2023 08:56:40 -0400 Subject: [PATCH 01/16] Check Catalyst 13 --- Project.toml | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/Project.toml b/Project.toml index f44b30d..1745673 100644 --- a/Project.toml +++ b/Project.toml @@ -19,17 +19,17 @@ Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7" TupleTools = "9d95972d-f1c8-5527-a6e0-b4b365fa01f6" [compat] -Catalyst = "11, 12" +Catalyst = "13" Combinatorics = "1" Cumulants = "1" DataStructures = "0.18" -Distributions = "0.23, 0.24, 0.25" +Distributions = "0.25" DocStringExtensions = "0.8, 0.9" Latexify = "0.15.14" -ModelingToolkit = "6.5, 7, 8" +ModelingToolkit = "8" SciMLBase = "1.8.3" -SymbolicUtils = "0.19" -Symbolics = "3, 4" +SymbolicUtils = "1" +Symbolics = "5" TupleTools = "1.2.0" julia = "1.6" From 44ccee49405fa48653898025b559bce6e7b168f8 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Wed, 19 Apr 2023 16:36:58 +0200 Subject: [PATCH 02/16] make sure to use BasicSymbolic --- src/symbolic.jl | 32 ++++++++++++++++---------------- 1 file changed, 16 insertions(+), 16 deletions(-) diff --git a/src/symbolic.jl b/src/symbolic.jl index 83c9f15..0113590 100644 --- a/src/symbolic.jl +++ b/src/symbolic.jl @@ -20,13 +20,13 @@ end trim_key(expr) = filter(x -> !(isspace(x) || x == ')' || x== '(' || x==','), string(expr)) # Expand a symbolic expression (no binomial expansion) -expansion_rule_mod = @acrule ~x * +(~~ys) => sum(map(y-> ~x * y, ~~ys)) # apply distribution law +expansion_rule_mod = @acrule ~x * +(~~ys) => sum(map(y-> ~x * y, ~~ys)) # apply distribution law expand_mod = Fixpoint(Prewalk(PassThrough(expansion_rule_mod))) # distributes terms until no longer possible -flatten_rule_mod = @rule(~x::isnotflat(+) => flatten_term(+, ~x)) # -flatten_mod = Fixpoint(PassThrough(flatten_rule_mod)) # +flatten_rule_mod = @rule(~x::isnotflat(+) => flatten_term(+, ~x)) # +flatten_mod = Fixpoint(PassThrough(flatten_rule_mod)) # expand_expr = Fixpoint(PassThrough(Chain([expand_mod, flatten_mod]))) # apply flatten and distribution until no longer possible -function define_μ(iter::AbstractVector, iv::Union{Sym,Num}) +function define_μ(iter::AbstractVector, iv::SymbolicUtils.BasicSymbolic) indices = map(trim_key, iter) @@ -47,15 +47,15 @@ function define_μ(iter::AbstractVector, iv::Union{Sym,Num}) end -define_μ(N::Int, order::Int, iv::Union{Sym,Num}) = define_μ(construct_iter_all(N, order), iv) +define_μ(N::Int, order::Int, iv::SymbolicUtils.BasicSymbolic) = define_μ(construct_iter_all(N, order), iv) -function define_μ(N::Int, order::Int, iter = construct_iter_all(N, order)) +function define_μ(N::Int, order::Int, iter = construct_iter_all(N, order)) @parameters t return define_μ(iter, value(t)) end -function define_M(iter::AbstractVector, iv::Union{Sym,Num}) +function define_M(iter::AbstractVector, iv::SymbolicUtils.BasicSymbolic) indices = map(trim_key, iter) @@ -78,15 +78,15 @@ function define_M(iter::AbstractVector, iv::Union{Sym,Num}) end -define_M(N::Int, order::Int, iv::Union{Sym,Num}) = define_M(construct_iter_all(N, order), iv) +define_M(N::Int, order::Int, iv::SymbolicUtils.BasicSymbolic) = define_M(construct_iter_all(N, order), iv) -function define_M(N::Int, order::Int, iter = construct_iter_all(N, order)) +function define_M(N::Int, order::Int, iter = construct_iter_all(N, order)) @parameters t return define_M(iter, value(t)) end function extract_variables(eqs::Array{Equation, 1}, μ, M=[]) - + vars = vcat(values(μ)..., values(M)...) # extract variables from rhs of each equation eq_vars = unique(vcat(get_variables.(eqs)...)) @@ -158,13 +158,13 @@ end function split_factor(expr::Div, iv, smap, vars) num, denom = arguments(expr) isconstant(denom, vars, iv) || error("The denominator $denom is not constant.") - + factor, powers = split_factor(num, iv, smap, vars) - + factor / denom, powers end -function split_factor(expr, iv, smap, vars) +function split_factor(expr, iv, smap, vars) if isconstant(expr, vars, iv) expr, zeros(length(vars)) elseif isvar(expr, vars) @@ -184,7 +184,7 @@ end function polynomial_propensity(expr::Add, iv, smap, vars) factors = [] powers = Vector{Int}[] - + for term in arguments(expr) factor_term, power_term = try split_factor(term, iv, smap, vars) @@ -208,7 +208,7 @@ function polynomial_propensity(expr, iv, smap, vars) [ factor ], [ powers ] end -function polynomial_propensities(arr::AbstractArray, iv::Sym, smap::AbstractDict) +function polynomial_propensities(arr::AbstractArray, iv::SymbolicUtils.BasicSymbolic, smap::AbstractDict) vars = [ x for (x,_) in Base.sort(collect(smap), by=x->x[2]) ] all_factors = Array{Vector}(undef, size(arr)) @@ -221,4 +221,4 @@ function polynomial_propensities(arr::AbstractArray, iv::Sym, smap::AbstractDict max_power = maximum(sum.(vcat(all_powers...))) all_factors, all_powers, max_power -end \ No newline at end of file +end From 140bde40b94541291ec53d0fa8ac324b31e64ccd Mon Sep 17 00:00:00 2001 From: Augustinas Sukys Date: Tue, 2 May 2023 19:35:38 +0100 Subject: [PATCH 03/16] fix symbolics --- src/MomentClosure.jl | 8 ++-- src/stochastic_stoichiometry.jl | 45 ++++++++++++++----- src/symbolic.jl | 80 ++++++++++++++++++--------------- 3 files changed, 84 insertions(+), 49 deletions(-) diff --git a/src/MomentClosure.jl b/src/MomentClosure.jl index 448d118..9b6f00c 100644 --- a/src/MomentClosure.jl +++ b/src/MomentClosure.jl @@ -7,11 +7,11 @@ using SciMLBase, SciMLBase.EnsembleAnalysis using Random using Distributions: Geometric -using Symbolics: value, var_from_nested_derivative, map_subscripts, hessian, gradient, setmetadata, scalarize +using Symbolics: value, var_from_nested_derivative, map_subscripts, hessian, gradient, setmetadata using SymbolicUtils.Rewriters: Chain, PassThrough, Prewalk, Fixpoint -using SymbolicUtils: Symbolic, Term, Real, expand, simplify, operation, - arguments, @rule, @acrule, isnotflat, flatten_term, - istree, FnType, Mul, Add, Pow, Div +using SymbolicUtils: BasicSymbolic, Real, Term, FnType, expand, simplify, + operation, arguments, @rule, @acrule, isnotflat, flatten_term, + istree, isterm, ismul, isadd, ispow, isdiv using DataStructures: OrderedDict using TupleTools: sort diff --git a/src/stochastic_stoichiometry.jl b/src/stochastic_stoichiometry.jl index 8e7538e..505b800 100644 --- a/src/stochastic_stoichiometry.jl +++ b/src/stochastic_stoichiometry.jl @@ -4,21 +4,31 @@ # TODO: follow changes in Catalyst and extend the compatibility accordingly # TODO: extend for other univariate discrete distributions -expected_coeff(x, n::Int) = x^n # covers Number/Sym/Pow -expected_coeff(x::Mul, n::Int) = prod(expected_coeff(a, n) for a in arguments(x)) +expected_coeff_mul(x) = prod(expected_coeff(a) for a in arguments(x)) +expected_coeff_pow(x) = expected_coeff(arguments(x)...) +expected_coeff_term(x) = expected_coeff(x, 1) -function expected_coeff(x::Div, n::Int) +function expected_coeff(x) + if ismul(x) + expected_coeff_mul(x) + elseif ispow(x) + expected_coeff_pow(x) + elseif isterm(x) + expected_coeff_term(x) + else + x + end +end + +expected_coeff_add(x, n::Int) = sum( expected_coeff(a) for a in arguments(expand(x^n)) ) +expected_coeff_mul(x, n::Int) = prod(expected_coeff(a, n) for a in arguments(x)) + +function expected_coeff_div(x, n::Int) num, denom = arguments(x) (expected_coeff(num, n) / expected_coeff(denom, n)) end -expected_coeff(x::Add, n::Int) = sum( expected_coeff(a) for a in arguments(expand(x^n)) ) -expected_coeff(x::Mul) = prod(expected_coeff(a) for a in arguments(x)) -expected_coeff(x::Pow) = expected_coeff(arguments(x)...) -expected_coeff(x::Term) = expected_coeff(x, 1) -expected_coeff(x) = x - -function expected_coeff(x::Term, n::Int) +function expected_coeff_term(x, n::Int) if isequal(operation(x), rand) args = arguments(x) isone(length(args)) || error("Unexpected arguments in $x") @@ -28,6 +38,21 @@ function expected_coeff(x::Term, n::Int) end end +function expected_coeff(x, n::Int) + if ismul(x) + expected_coeff_mul(x, n) + elseif isdiv(x) + expected_coeff_div(x, n) + elseif isadd(x) + expected_coeff_add(x, n) + elseif isterm(x) + expected_coeff_term(x, n) + else + x^n + end +end + + function resolve_moment(d, args, n) @assert isequal(d, Geometric) "Can only compute moments for the geometric distribution." geometric_raw_moment(args[1], n) diff --git a/src/symbolic.jl b/src/symbolic.jl index 0113590..67c9324 100644 --- a/src/symbolic.jl +++ b/src/symbolic.jl @@ -26,7 +26,7 @@ flatten_rule_mod = @rule(~x::isnotflat(+) => flatten_term(+, ~x)) # flatten_mod = Fixpoint(PassThrough(flatten_rule_mod)) # expand_expr = Fixpoint(PassThrough(Chain([expand_mod, flatten_mod]))) # apply flatten and distribution until no longer possible -function define_μ(iter::AbstractVector, iv::SymbolicUtils.BasicSymbolic) +function define_μ(iter::AbstractVector, iv::BasicSymbolic) indices = map(trim_key, iter) @@ -47,7 +47,7 @@ function define_μ(iter::AbstractVector, iv::SymbolicUtils.BasicSymbolic) end -define_μ(N::Int, order::Int, iv::SymbolicUtils.BasicSymbolic) = define_μ(construct_iter_all(N, order), iv) +define_μ(N::Int, order::Int, iv::BasicSymbolic) = define_μ(construct_iter_all(N, order), iv) function define_μ(N::Int, order::Int, iter = construct_iter_all(N, order)) @parameters t @@ -55,7 +55,7 @@ function define_μ(N::Int, order::Int, iter = construct_iter_all(N, order)) end -function define_M(iter::AbstractVector, iv::SymbolicUtils.BasicSymbolic) +function define_M(iter::AbstractVector, iv::BasicSymbolic) indices = map(trim_key, iter) @@ -78,7 +78,7 @@ function define_M(iter::AbstractVector, iv::SymbolicUtils.BasicSymbolic) end -define_M(N::Int, order::Int, iv::SymbolicUtils.BasicSymbolic) = define_M(construct_iter_all(N, order), iv) +define_M(N::Int, order::Int, iv::BasicSymbolic) = define_M(construct_iter_all(N, order), iv) function define_M(N::Int, order::Int, iter = construct_iter_all(N, order)) @parameters t @@ -103,13 +103,15 @@ end # a simple hack to convert SymbolicUtils Div (fractions) to Pow (multiplication) # more stable in symbolic manipulations that we perform -function div_to_pow(expr::Div) - args = arguments(expr) - args[1] * args[2]^-1 +function div_to_pow(expr) + if isdiv(expr) + args = arguments(expr) + args[1] * args[2]^-1 + else + expr + end end -div_to_pow(expr) = expr - ## Set of functions to deconstruct polynomial propensities ## #= @@ -133,21 +135,20 @@ isvar(x, vars) = any(isequal(x), vars) isconstant(expr, vars, iv) = !istree(expr) || (!isvar(expr, vars) && all(arg -> isconstant(arg, vars, iv), arguments(expr))) -function split_factor(expr::Pow, iv, smap, vars) +function split_factor_pow(expr, iv, vars) base, exp = arguments(expr) - #(exp isa Int && exp >= 0) || error("Unexpected exponent: $expr") exp isa Int || error("Unexpected exponent: $expr") - factor, powers = split_factor(base, iv, smap, vars) + factor, powers = split_factor(base, iv, vars) factor ^ exp, powers .* exp end -function split_factor(expr::Mul, iv, smap, vars) - powers = zeros(Int, length(smap)) +function split_factor_mul(expr, iv, vars) + powers = zeros(Int, length(vars)) factor = 1 for arg in arguments(expr) - factor_arg, power_arg = split_factor(arg, iv, smap, vars) + factor_arg, power_arg = split_factor(arg, iv, vars) factor *= factor_arg powers .+= power_arg end @@ -155,17 +156,24 @@ function split_factor(expr::Mul, iv, smap, vars) factor, powers end -function split_factor(expr::Div, iv, smap, vars) +function split_factor_div(expr, iv, vars) num, denom = arguments(expr) - isconstant(denom, vars, iv) || error("The denominator $denom is not constant.") + isconstant(denom, vars, iv) || error("The denominator $denom in propensity $expr is not constant.") - factor, powers = split_factor(num, iv, smap, vars) + factor, powers = split_factor(num, iv, vars) factor / denom, powers end -function split_factor(expr, iv, smap, vars) - if isconstant(expr, vars, iv) +function split_factor(expr, iv, vars) + + if ispow(expr) + split_factor_pow(expr, iv, vars) + elseif ismul(expr) + split_factor_mul(expr, iv, vars) + elseif isdiv(expr) + split_factor_div(expr, iv, vars) + elseif isconstant(expr, vars, iv) expr, zeros(length(vars)) elseif isvar(expr, vars) 1, map(isequal(expr), vars) @@ -174,20 +182,21 @@ function split_factor(expr, iv, smap, vars) end end -function polynomial_propensity(expr::Div, iv, smap, vars) + +function polynomial_propensity_div(expr, iv, vars) num, denom = arguments(expr) - isconstant(denom, vars, iv) || error("The denominator $denom is not constant.") - factors, powers = polynomial_propensity(num, iv, smap, vars) - factors ./ denom, powers + isconstant(denom, vars, iv) || error("The denominator $denom in propensity $expr is not constant.") + factors, powers = polynomial_propensity(num, iv, vars) + factors ./denom, powers end -function polynomial_propensity(expr::Add, iv, smap, vars) +function polynomial_propensity_add(expr, iv, vars) factors = [] powers = Vector{Int}[] for term in arguments(expr) factor_term, power_term = try - split_factor(term, iv, smap, vars) + split_factor(term, iv, vars) catch e error("Propensity function $term is non-polynomial? \n" * string(e)) end @@ -198,17 +207,18 @@ function polynomial_propensity(expr::Add, iv, smap, vars) factors, powers end -function polynomial_propensity(expr, iv, smap, vars) - factor, powers = try - split_factor(expr, iv, smap, vars) - catch e - error("Propensity function $expr is non-polynomial? \n" * string(e)) +function polynomial_propensity(expr, iv, vars) + if isdiv(expr) + polynomial_propensity_div(expr, iv, vars) + elseif isadd(expr) + polynomial_propensity_add(expr, iv, vars) + else + factor, powers = split_factor(expr, iv, vars) + [factor], [powers] end - - [ factor ], [ powers ] end -function polynomial_propensities(arr::AbstractArray, iv::SymbolicUtils.BasicSymbolic, smap::AbstractDict) +function polynomial_propensities(arr::AbstractArray, iv::BasicSymbolic, smap::AbstractDict) vars = [ x for (x,_) in Base.sort(collect(smap), by=x->x[2]) ] all_factors = Array{Vector}(undef, size(arr)) @@ -216,7 +226,7 @@ function polynomial_propensities(arr::AbstractArray, iv::SymbolicUtils.BasicSymb for (rind, expr) in enumerate(arr) expr = expand(expr) - all_factors[rind], all_powers[rind] = polynomial_propensity(expr, iv, smap, vars) + all_factors[rind], all_powers[rind] = polynomial_propensity(expr, iv, vars) end max_power = maximum(sum.(vcat(all_powers...))) From 993217e800e9cef532bd4c50d2a0a0975765ea66 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Tue, 2 May 2023 15:30:27 -0400 Subject: [PATCH 04/16] Update reaction_system_tests.jl --- test/reaction_system_tests.jl | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/test/reaction_system_tests.jl b/test/reaction_system_tests.jl index 97bd7cb..089dbf0 100644 --- a/test/reaction_system_tests.jl +++ b/test/reaction_system_tests.jl @@ -11,11 +11,12 @@ S_mat = [ 1 -1 1 -1; a = [c₁*X*Y*(X-1)/Ω^2, c₂*X, c₃*Ω, c₄*X] rn = @reaction_network begin + @parameters c₁ c₂ c₃ c₄ Ω (c₁/Ω^2), 2X + Y → 3X (c₂), X → Y (Ω*c₃, c₄), 0 ↔ X -end c₁ c₂ c₃ c₄ Ω +end smap = speciesmap(rn) @test isequal(get_stoichiometry(rn, smap), S_mat) -@test isequal(propensities(rn, combinatoric_ratelaw=false), a) \ No newline at end of file +@test isequal(propensities(rn, combinatoric_ratelaw=false), a) From 9bf8aa46a625522c29df21900c4164d14b3df038 Mon Sep 17 00:00:00 2001 From: Augustinas Sukys Date: Thu, 4 May 2023 23:21:41 +0100 Subject: [PATCH 05/16] simplify less extremely long expressions + simplify can get very slow here --- src/closure_methods/conditional_derivative_matching.jl | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/src/closure_methods/conditional_derivative_matching.jl b/src/closure_methods/conditional_derivative_matching.jl index 6d664bb..0a97673 100644 --- a/src/closure_methods/conditional_derivative_matching.jl +++ b/src/closure_methods/conditional_derivative_matching.jl @@ -163,11 +163,10 @@ function conditional_derivative_matching(sys::MomentEquations, raw_to_central_exp = raw_to_central_moments(N, sys.q_order, μ_M_exp, bernoulli=true) raw_to_central = raw_to_central_moments(N, sys.q_order, μ_M, bernoulli=true) for i in sys.iter_q - closure_exp[sys.M[i]] = simplify(raw_to_central_exp[i]) - closure[sys.M[i]] = simplify(raw_to_central[i]) + closure_exp[sys.M[i]] = raw_to_central_exp[i] + closure[sys.M[i]] = raw_to_central[i] end - else closure_exp = closure_μ_exp closure = closure_μ From 8787ffb898822bc7c57310a6b839e8eeb51433b7 Mon Sep 17 00:00:00 2001 From: Augustinas Sukys Date: Thu, 4 May 2023 23:25:02 +0100 Subject: [PATCH 06/16] expand symbolic div for latexify makes LaTeX expressions a bit easier to read in certain contexts --- src/symbolic.jl | 1 + src/utils.jl | 2 ++ 2 files changed, 3 insertions(+) diff --git a/src/symbolic.jl b/src/symbolic.jl index 67c9324..fe88e9f 100644 --- a/src/symbolic.jl +++ b/src/symbolic.jl @@ -25,6 +25,7 @@ expand_mod = Fixpoint(Prewalk(PassThrough(expansion_rule_mod))) # distributes te flatten_rule_mod = @rule(~x::isnotflat(+) => flatten_term(+, ~x)) # flatten_mod = Fixpoint(PassThrough(flatten_rule_mod)) # expand_expr = Fixpoint(PassThrough(Chain([expand_mod, flatten_mod]))) # apply flatten and distribution until no longer possible +expand_div = PassThrough(@acrule( +(~~xs) / ~a => sum(map(x -> x / ~a, ~~xs)))) function define_μ(iter::AbstractVector, iv::BasicSymbolic) diff --git a/src/utils.jl b/src/utils.jl index 1c611c9..b915655 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -307,6 +307,7 @@ function format_moment_eqs(eqs::MomentEquations) for i in 1:size(odes)[1] key = vars[i] eq = odes[i].rhs + eq = expand_div(eq) expr = "d"*string(key)*"/dt = "*string(eq) expr = replace(expr, "(t)"=>"") expr = replace(expr, ".0"=>"") @@ -341,6 +342,7 @@ function format_closure(eqs::ClosedMomentEquations; format_all::Bool=false) for i in iter eq = closure[i] + eq = expand_div(eq) expr = string(i)*" = "*string(eq) expr = replace(expr, "(t)"=>"") expr = replace(expr, ".0"=>"") From 5ec637b597459730211b259295af55a562877c22 Mon Sep 17 00:00:00 2001 From: Augustinas Sukys Date: Thu, 4 May 2023 23:26:17 +0100 Subject: [PATCH 07/16] add missing imports --- src/MomentClosure.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/MomentClosure.jl b/src/MomentClosure.jl index 9b6f00c..1da2453 100644 --- a/src/MomentClosure.jl +++ b/src/MomentClosure.jl @@ -7,7 +7,8 @@ using SciMLBase, SciMLBase.EnsembleAnalysis using Random using Distributions: Geometric -using Symbolics: value, var_from_nested_derivative, map_subscripts, hessian, gradient, setmetadata +using Symbolics: value, var_from_nested_derivative, map_subscripts, hessian, + gradient, setmetadata, scalarize using SymbolicUtils.Rewriters: Chain, PassThrough, Prewalk, Fixpoint using SymbolicUtils: BasicSymbolic, Real, Term, FnType, expand, simplify, operation, arguments, @rule, @acrule, isnotflat, flatten_term, From ec4fb74f71b8169cd8c0f3ce7022b4295eb2818c Mon Sep 17 00:00:00 2001 From: Augustinas Sukys Date: Thu, 4 May 2023 23:28:03 +0100 Subject: [PATCH 08/16] fix the remaining tests --- src/central_moment_equations.jl | 4 ++-- src/raw_moment_equations.jl | 2 +- test/bernoulli_closure_tests.jl | 10 +++++----- test/closure_tests.jl | 15 ++++++--------- test/conditional_closure_tests.jl | 8 ++++---- test/latexify_tests.jl | 9 ++++++++- test/lma_tests.jl | 10 +++++----- test/moment_equations_SDE_tests.jl | 6 ++++-- test/moment_equations_tests.jl | 24 +++++++++++++----------- test/reaction_system_tests.jl | 2 +- 10 files changed, 49 insertions(+), 41 deletions(-) diff --git a/src/central_moment_equations.jl b/src/central_moment_equations.jl index 57f2809..9dd528f 100644 --- a/src/central_moment_equations.jl +++ b/src/central_moment_equations.jl @@ -128,7 +128,7 @@ function generate_central_moment_eqs(rn::ReactionSystem, m_order::Int, q_order:: else du[i] = S[i, r]*suma + du[i] end - du[i] = expand(du[i]) + du[i] = simplify(du[i]) end end @@ -158,7 +158,7 @@ function generate_central_moment_eqs(rn::ReactionSystem, m_order::Int, q_order:: dM[i] -= i[j]*du[j]*M[i.-iter_1[j]] end end - dM[i] = expand(dM[i]) + dM[i] = simplify(dM[i]) end D = Differential(iv) diff --git a/src/raw_moment_equations.jl b/src/raw_moment_equations.jl index 8390c19..d43ab45 100644 --- a/src/raw_moment_equations.jl +++ b/src/raw_moment_equations.jl @@ -72,7 +72,7 @@ function generate_raw_moment_eqs(rn::ReactionSystem, m_order::Int; dμ[i] += factor_j * suma end end - dμ[i] = expand(dμ[i]) + dμ[i] = simplify(dμ[i]) end D = Differential(iv) diff --git a/test/bernoulli_closure_tests.jl b/test/bernoulli_closure_tests.jl index ae25387..13e2d34 100644 --- a/test/bernoulli_closure_tests.jl +++ b/test/bernoulli_closure_tests.jl @@ -13,7 +13,7 @@ rn = @reaction_network begin k_off*P^2, g --> 0 k_p, g --> g + $m*P γ_p, P --> 0 -end k_on k_off k_p γ_p +end binary_vars = [1] @@ -30,7 +30,7 @@ closed_eqs = moment_closure(sys, "normal", binary_vars) closed_eqs= moment_closure(sys, "log-normal", binary_vars) expr1 = expand(closed_eqs.closure[M[0,3]]) expr2 = M[0,2]^3*μ[0,1]^-3 + 3*M[0,2]^2*μ[0,1]^-1 -@test length(closed_eqs.odes.states) == 4 && isequal(expand(expr1), expr2) +@test length(closed_eqs.odes.states) == 4 && isequal(simplify(expr1), simplify(expr2)) closed_eqs = moment_closure(sys, "poisson", binary_vars) @test length(closed_eqs.odes.states) == 4 && isequal(closed_eqs.closure[M[0,3]], μ[0,1]) @@ -43,7 +43,7 @@ expr2 = 2*M[0,2]^2*μ[0,1]^-1 closed_eqs = moment_closure(sys, "derivative matching", binary_vars) expr1 = expand(closed_eqs.closure[sys.M[0,3]]) expr2 = M[0,2]^3*μ[0,1]^-3 + 3*M[0,2]^2*μ[0,1]^-1 -@test length(closed_eqs.odes.states) == 4 && isequal(expr1, expr2) +@test length(closed_eqs.odes.states) == 4 && isequal(simplify(expr1), simplify(expr2)) sys = generate_raw_moment_eqs(rn, 2) @@ -60,7 +60,7 @@ closed_eqs = moment_closure(sys, "poisson", binary_vars) @test length(closed_eqs.odes.states) == 4 closed_eqs = moment_closure(sys, "gamma", binary_vars) -@test isequal(expr1, expr2) +@test length(closed_eqs.odes.states) == 4 closed_eqs = moment_closure(sys, "derivative matching", binary_vars) -@test isequal(expr1, expr2) +@test length(closed_eqs.odes.states) == 4 \ No newline at end of file diff --git a/test/closure_tests.jl b/test/closure_tests.jl index adf8d73..2b9666d 100644 --- a/test/closure_tests.jl +++ b/test/closure_tests.jl @@ -8,9 +8,8 @@ rn = @reaction_network begin (c₁/Ω^2), 2X + Y → 3X (c₂), X → Y (Ω*c₃, c₄), 0 ↔ X -end c₁ c₂ c₃ c₄ Ω +end -#@parameters c₁, c₂, c₃, c₄, Ω @syms c₁::Real c₂::Real c₃::Real c₄::Real Ω::Real μ = define_μ(2,4) M = define_M(2,4) @@ -24,8 +23,7 @@ closed_eqs = moment_closure(sys, "zero") expr1 = closed_eqs.odes.eqs[1].rhs expr2 = c₃*Ω + M[1,1]*c₁*μ[1,0]*(Ω^-2) + M[1,1]*c₁*(Ω^-2)*(μ[1,0]- 1) + c₁*M[2,0]*μ[0,1]*(Ω^-2) + c₁*μ[0,1]*μ[1,0]*(Ω^-2)*(μ[1,0] - 1) - c₂*μ[1,0] - c₄*μ[1,0] -expr2 = value.(expand(expr2)) -@test isequal(expand(expr1), expand(expr2)) +@test isequal(simplify(expr1), simplify(expr2)) # check that deterministic_IC is working with central moments ic_values = Dict(deterministic_IC([2, 5], closed_eqs)) @@ -39,7 +37,7 @@ closed_eqs= moment_closure(sys, "log-normal") expr1 = closed_eqs.closure[M[1,2]] expr2 = μ[1,0]*μ[0,1]^2*(1.0+M[0,2]*μ[0,1]^-2)*(1.0 + M[1,1]*(μ[0,1]^-1)*(μ[1,0]^-1))^2 - M[0,2]*μ[1,0] - μ[1,0]*μ[0,1]^2 - 2*M[1,1]*μ[0,1] -@test isequal(expand(expr1), expand(expr2)) +@test isequal(expr1, simplify(expand(expr2))) closed_eqs = moment_closure(sys, "poisson") @test isequal(closed_eqs.closure[M[3,0]], μ[1,0]) @@ -54,7 +52,7 @@ closed_eqs = moment_closure(sys, "derivative matching") expr1 = closed_eqs.closure[sys.M[0,4]] expr2 = μ[0,1]^4*(M[0,2]+μ[0,1]^2)^-6*(M[0,3]+μ[0,1]^3+3*M[0,2]*μ[0,1])^4 - μ[0,1]^4 - 6*M[0,2]*μ[0,1]^2 - 4*M[0,3]*μ[0,1] -@test isequal(expr1, expr2) +@test isequal(expand(expr1), expand(simplify(expr2))) # --- Test closures on raw moment equations --- @@ -72,8 +70,7 @@ expr1 = closed_eqs.odes.eqs[5].rhs expr2 = c₂*μ[1,0] + 2*c₂*μ[1,1] + c₁*μ[0,1]*μ[2,0]*Ω^-2 - c₁*μ[1,1]*Ω^-2 - 4*c₁*μ[1,1]^2*Ω^-2 + 4*c₁*μ[0,1]*μ[1,1]*Ω^-2 - 2*c₁*μ[0,1]*μ[1,0]^2*Ω^-2 + 2*c₁*μ[0,2]*μ[1,0]*Ω^-2 -2*c₁*μ[0,2]*μ[2,0]*Ω^-2 + 2*c₁*μ[1,0]*μ[1,1]*Ω^-2 - 4*c₁*μ[1,0]*μ[0,1]^2*Ω^-2 + 4*c₁*μ[0,1]^2*μ[1,0]^2*Ω^-2 -expr2 = simplify(value.(expr2)) -@test isequal(expand(expr1), expr2) +@test isequal(simplify(expr1), simplify(expr2)) # check that deterministic_IC is working with raw moments ic_values = Dict(deterministic_IC([2, 5], closed_eqs)) @@ -93,7 +90,7 @@ expr2 = μ[1,0] + 6 * μ[1,0]^4 + 3*μ[2,0]^2 + 4*μ[1,0]*μ[3,0] -12*μ[2,0]*μ closed_eqs = moment_closure(sys, "gamma") expr1 = closed_eqs.closure[μ[0,3]] expr2 = μ[0,1]^3 + 3*μ[0,1]*(μ[0,2]-μ[0,1]^2) + 2*μ[0,1]^-1*(μ[0,2]-μ[0,1]^2)^2 -@test isequal(expr1, expr2) +@test isequal(simplify(expr1), simplify(expr2)) closed_eqs = moment_closure(sys, "derivative matching") expr1 = closed_eqs.closure[sys.μ[0,3]] diff --git a/test/conditional_closure_tests.jl b/test/conditional_closure_tests.jl index bca71c8..717381a 100644 --- a/test/conditional_closure_tests.jl +++ b/test/conditional_closure_tests.jl @@ -14,7 +14,7 @@ rn = @reaction_network begin k_off*P^2, g --> 0 k_p, g --> g + $m*P γ_p, P --> 0 -end k_on k_off k_p γ_p +end binary_vars = [1] @parameters k_on, k_off, k_p, γ_p @@ -39,10 +39,10 @@ expr2 = simplify(value.(expr2)) closed_eqs = moment_closure(sys, "conditional gaussian", binary_vars) expr1 = closed_eqs.closure[μ[1,4]] expr2 = 4*μ[1,1]*μ[1,3]*μ[1,0]^-1 + 3*μ[1,2]^2*μ[1,0]^-1 - 12*μ[1,2]*μ[1,1]^2*μ[1,0]^-2 + 6*μ[1,1]^4*μ[1,0]^-3 -@test isequal(expr1, expr2) +@test isequal(expr1, simplify(expr2)) expr1 = closed_eqs.closure[μ[1,3]] expr2 = 3*μ[1,2]*μ[1,1]*μ[1,0]^-1 - 2*μ[1,1]^3*μ[1,0]^-2 -@test isequal(expr1, expr2) +@test isequal(expr1, simplify(expr2)) closed_eqs = moment_closure(sys, "conditional derivative matching", binary_vars) expr1 = closed_eqs.closure[μ[1,4]] @@ -74,4 +74,4 @@ closed_eqs = moment_closure(sys, "conditional derivative matching", binary_vars) expr1 = closed_eqs.closure[M[1,3]] expr2 = μ[1,0]*(M[1,1]+μ[0,1]*μ[1,0])^-3*(M[1,2]+M[0,2]*μ[1,0]+μ[1,0]*μ[0,1]^2+2*M[1,1]*μ[0,1])^3 - M[0,3]*μ[1,0] - 3*M[1,1]*μ[0,1]^2 - 3*M[1,2]*μ[0,1] - μ[1,0]*μ[0,1]^3- 3*M[0,2]*μ[0,1]*μ[1,0] -@test isequal(expand(expr1), expand(expr2)) +@test isequal(expr1, simplify(expr2)) \ No newline at end of file diff --git a/test/latexify_tests.jl b/test/latexify_tests.jl index 86687fa..8a6480d 100644 --- a/test/latexify_tests.jl +++ b/test/latexify_tests.jl @@ -13,7 +13,7 @@ rn = @reaction_network begin k_off*P^2, g --> 0 k_p, g --> g + $m*P γ_p, P --> 0 -end k_on k_off k_p γ_p +end binary_vars = [1] raw_eqs = generate_raw_moment_eqs(rn, 2, combinatoric_ratelaw=false) @@ -37,11 +37,18 @@ expr = replace(raw"\begin{align*} \frac{d\mu_{1 0}}{dt} =& k_{on}", "\r\n"=>"\n") @test latexify(clean_eqs)[1:46] == expr +#= expr = replace(raw"\begin{align*} \mu_{1 2} =& \mu_{1 0}^{-1} \mu_{1 1}^{2} \\ \mu_{1 3} =& 3 \mu_{1 0}^{-1} \mu_{1 1} \mu_{1 2} - 2 \mu_{1 0}^{-2} \mu_{1 1}^{3} \end{align*} ", "\r\n"=>"\n") +=# +expr = replace(raw"\begin{align*} +\mu_{1 2} =& \frac{\mu_{1 1}^{2}}{\mu_{1 0}} \\ +\mu_{1 3} =& \frac{-2 \mu_{1 1}^{3}}{\mu_{1 0}^{2}} + \frac{3 \mu_{1 1} \mu_{1 2}}{\mu_{1 0}} +\end{align*} +", "\r\n"=>"\n") @test latexify(closed_raw_eqs, :closure) == expr @test_throws MethodError latexify(raw_eqs, :closure) \ No newline at end of file diff --git a/test/lma_tests.jl b/test/lma_tests.jl index 59c90a3..5399c64 100644 --- a/test/lma_tests.jl +++ b/test/lma_tests.jl @@ -14,7 +14,7 @@ rn_nonlinear = @reaction_network begin ρ_u, g → g + p ρ_b*(1-g), 0 ⇒ p 1, p → 0 -end σ_b σ_u ρ_b ρ_u +end rn_linear = @reaction_network begin σ_b_LMA, g → 0 @@ -22,7 +22,7 @@ rn_linear = @reaction_network begin ρ_u, g → g+p (ρ_b*(1-g)), 0 ⇒ p 1, p → 0 -end σ_b_LMA σ_u ρ_b ρ_u +end binary_vars = [speciesmap(rn_nonlinear)[g]] LMA_eqs, _ = linear_mapping_approximation(rn_nonlinear, rn_linear, binary_vars, combinatoric_ratelaw=false) @@ -36,7 +36,7 @@ expr2 = ρ_b + ρ_u*μ[1,0] - ρ_b*μ[1,0] - μ[0,1] @test isequal(expr1, value.(expr2)) expr1 = LMA_eqs.odes.eqs[3].rhs expr2 = ρ_u*μ[1,0] + σ_u*μ[0,1] - μ[1,1] - σ_u*μ[1,1] - σ_b*μ[1,1]^2*μ[1,0]^-1 -@test isequal(expr1, value.(expr2)) +@test isequal(simplify(expr1), simplify(expr2)) # cooperativity cp=2 rn_nonlinear = @reaction_network begin @@ -45,10 +45,10 @@ rn_nonlinear = @reaction_network begin ρ_u, g → g + p ρ_b*(1-g), 0 ⇒ p 1, p → 0 -end σ_b σ_u ρ_b ρ_u +end binary_vars = [speciesmap(rn_nonlinear)[g]] _, effective_params = linear_mapping_approximation(rn_nonlinear, rn_linear, binary_vars, combinatoric_ratelaw=false) expr1 = effective_params[value(σ_b_LMA)] expr2 = (σ_b*μ[1,2] - σ_b*μ[1,1]) * μ[1,0]^-1 -@test isequal(expr1, value.(expr2)) +@test isequal(expr1, value(expr2)) \ No newline at end of file diff --git a/test/moment_equations_SDE_tests.jl b/test/moment_equations_SDE_tests.jl index 26f9c29..111ccbe 100644 --- a/test/moment_equations_SDE_tests.jl +++ b/test/moment_equations_SDE_tests.jl @@ -24,11 +24,12 @@ cir_central_moments = generate_central_moment_eqs(cir_model, 2) ## Chemical Reaction Networks via Chemical Langevin Equation # unimolecular system schloegl = @reaction_network begin + @parameters k1 k2 1.0, 2X → 3X 1.0, 3X → 2X k1, ∅ → X k2, X → ∅ -end k1 k2 +end k1, k2 = parameters(schloegl) for order in 2:6 @@ -51,10 +52,11 @@ expr = k1 + 9*M[(3,)] + k2*μ[(1,)] + μ[(1,)]^3 + 19*M[(2,)]*μ[(1,)] + μ[(1,) # check if things work for multiple species rn = @reaction_network begin + @parameters c1 c2 c3 c4 (c1), 2X + Y → 3X (c2), X → Y (c3, c4), 0 ↔ X -end c1 c2 c3 c4 +end c1, c2, c3, c4 = parameters(rn) diff --git a/test/moment_equations_tests.jl b/test/moment_equations_tests.jl index 24b0612..9f40e78 100644 --- a/test/moment_equations_tests.jl +++ b/test/moment_equations_tests.jl @@ -1,5 +1,5 @@ using MomentClosure -using MomentClosure: expected_coeff +using MomentClosure: expected_coeff, expand_div using Symbolics: value, simplify, expand using Distributions: Geometric using Catalyst @@ -17,8 +17,8 @@ b = value(b); d = value(d); p = value(p) @test isequal(expected_coeff(sin(b), 2), sin(b)^2) @register_symbolic Geometric(p) m = rand(Geometric(p)) -@test isequal(expected_coeff(m, 1), p^-1 - 1) -@test isequal(expected_coeff(m, 3), -1 + 7p^-1 - 12p^-2 + 6p^-3) +@test isequal(expand_div(expected_coeff(m, 1)), p^-1 - 1) +@test isequal(expand_div(expected_coeff(m, 3)), expand(-1 + 7/p - 12/p^2 + 6/p^3)) @register_symbolic Bernoulli(p) @test_throws AssertionError expected_coeff(rand(Bernoulli(p)), 1) @@ -30,7 +30,7 @@ rn = @reaction_network begin (c₁*Ω^-2), 2X + Y → 3X (c₂), X → Y (Ω*c₃, c₄), 0 ↔ X -end c₁ c₂ c₃ c₄ Ω +end sys = generate_central_moment_eqs(rn, 2, 4, combinatoric_ratelaw=false) expr1 = sys.odes.eqs[2].rhs @@ -41,12 +41,14 @@ M = sys.M expr2 = c₂*μ[1,0] + c₁*M[1,1]*(Ω^-2) + c₁*μ[0,1]*μ[1,0]*Ω^-2 - c₁*M[2,1]*Ω^-2 - 2*c₁*M[1,1]*μ[1,0]*Ω^-2 - c₁*M[2,0]*μ[0,1]*Ω^-2 - c₁*μ[0,1]*Ω^-2*μ[1,0]^2 expr2 = simplify(value.(expr2)) -@test isequal(expand(expr1), expr2) +#@test isequal(expand(expr1), expr2) +@test isequal(expr1, expr2) + sys = generate_central_moment_eqs(rn, 2, combinatoric_ratelaw=false) expr1 = sys.odes.eqs[2].rhs @test isequal(MomentClosure.Differential(t)(sys.μ[1,0]), sys.odes.eqs[1].lhs) -@test isequal(expand(expr1), expr2) +@test isequal(expr1, expr2) sys = generate_raw_moment_eqs(rn, 2, combinatoric_ratelaw=false) μ = sys.μ @@ -54,7 +56,7 @@ expr1 = sys.odes.eqs[4].rhs expr2 = c₂*μ[2,0] + c₁*μ[1,1]*Ω^-2 - c₁*μ[3,1]*Ω^-2 -c₂*(μ[1,0] + μ[1,1]) + c₃*μ[0,1]*Ω - c₄*μ[1,1] - c₁*μ[1,2]*Ω^-2 + c₁*μ[2,2]*Ω^-2 expr2 = simplify(value.(expr2)) -@test isequal(expand(expr1), expand(expr2)) +@test isequal(expr1, expr2) # corner case - a linear propensity with rate coefficient 1 rn = @reaction_network begin @@ -68,21 +70,21 @@ sys = generate_raw_moment_eqs(rn, 2) rn = @reaction_network begin (X^1.5), X ⇒ Y -end c₁ +end @test_throws ErrorException generate_raw_moment_eqs(rn, 2) rn = @reaction_network begin (c₁*Y^3*X^1.5), X ⇒ Y -end c₁ +end @test_throws ErrorException generate_raw_moment_eqs(rn, 2) rn = @reaction_network begin (c₁*sin(exp(X)+Y)), X ⇒ Y -end c₁ +end @test_throws ErrorException generate_raw_moment_eqs(rn, 2) rn = @reaction_network begin (c₁*sin(exp(t))+Y^(3+1)), X ⇒ Y -end c₁ +end eqs = generate_raw_moment_eqs(rn, 2) @test eqs.q_order == 5 \ No newline at end of file diff --git a/test/reaction_system_tests.jl b/test/reaction_system_tests.jl index 089dbf0..489d0fe 100644 --- a/test/reaction_system_tests.jl +++ b/test/reaction_system_tests.jl @@ -3,7 +3,7 @@ using MomentClosure using Catalyst @parameters t, c₁, c₂, c₃, c₄, Ω -@variables X(t), Y(t) +@species X(t), Y(t) S_mat = [ 1 -1 1 -1; -1 1 0 0] From 5befda9705c1aceedef33e1bb1c636d4bf181849 Mon Sep 17 00:00:00 2001 From: Augustinas Sukys Date: Thu, 4 May 2023 23:48:38 +0100 Subject: [PATCH 09/16] combinatoric_ratelaw -> combinatoric_ratelaws now fully consistent with Catalyst --- src/central_moment_equations.jl | 8 ++++---- src/closure_methods/linear_mapping_approximation.jl | 10 +++++----- src/raw_moment_equations.jl | 8 ++++---- src/reaction_systems.jl | 8 ++++---- test/bernoulli_closure_tests.jl | 2 +- test/closure_tests.jl | 4 ++-- test/latexify_tests.jl | 2 +- test/lma_tests.jl | 4 ++-- test/moment_equations_SDE_tests.jl | 10 +++++----- test/moment_equations_tests.jl | 6 +++--- test/reaction_system_tests.jl | 4 ++-- 11 files changed, 33 insertions(+), 33 deletions(-) diff --git a/src/central_moment_equations.jl b/src/central_moment_equations.jl index 9dd528f..0511274 100644 --- a/src/central_moment_equations.jl +++ b/src/central_moment_equations.jl @@ -17,7 +17,7 @@ end """ generate_central_moment_eqs(rn::ReactionSystem, m_order::Int, q_order::Int=0; - langevin::Bool=false, combinatoric_ratelaw::Bool=true, smap=speciesmap(rn)) + langevin::Bool=false, combinatoric_ratelaws::Bool=true, smap=speciesmap(rn)) Given a [`ReactionSystem`](https://catalyst.sciml.ai/stable/api/catalyst_api/#ModelingToolkit.ReactionSystem) return the [`CentralMomentEquations`](@ref) of the system generated up to `m_order`. @@ -32,7 +32,7 @@ Notes: - if `langevin=true`, instead of the Chemical Master Equation the Chemical Langevin Equation (diffusion approximation) is considered, and the moment equations are constructed from the corresponding SDE formulation. -- `combinatoric_ratelaw=true` uses binomials in calculating the propensity functions +- `combinatoric_ratelaws=true` uses binomials in calculating the propensity functions of a `ReactionSystem`, see the notes for [`ModelingToolkit.jumpratelaw`] (https://mtk.sciml.ai/stable/systems/ReactionSystem/#ModelingToolkit.jumpratelaw). *Note* that this field is irrelevant using `ReactionSystemMod` as then the @@ -42,12 +42,12 @@ Notes: accessible with [`Catalyst.speciesmap`](https://catalyst.sciml.ai/stable/api/catalyst_api/#Catalyst.speciesmap). """ function generate_central_moment_eqs(rn::ReactionSystem, m_order::Int, q_order::Int=0; - langevin::Bool=false, combinatoric_ratelaw::Bool=true, smap=speciesmap(rn)) + langevin::Bool=false, combinatoric_ratelaws::Bool=true, smap=speciesmap(rn)) N = numspecies(rn) # no. of molecular species in the network R = numreactions(rn) # no. of reactions in the network iv = get_iv(rn) # independent variable (usually time) - a = propensities(rn; combinatoric_ratelaw) # propensity functions of all reactions in the network + a = propensities(rn; combinatoric_ratelaws) # propensity functions of all reactions in the network a = div_to_pow.(a) # convert Div to Pow (more stable at the moment than dealing with symbolic fractions) S = get_stoichiometry(rn, smap) # net stoichiometric matrix diff --git a/src/closure_methods/linear_mapping_approximation.jl b/src/closure_methods/linear_mapping_approximation.jl index 0a5859c..43d32a3 100644 --- a/src/closure_methods/linear_mapping_approximation.jl +++ b/src/closure_methods/linear_mapping_approximation.jl @@ -1,6 +1,6 @@ """ linear_mapping_approximation(rn_nonlinear::T, rn_linear::T, binary_vars::Array{Int,1}=Int[], m_order::Int=0; - combinatoric_ratelaw = true) where T <: ReactionSystem + combinatoric_ratelaws = true) where T <: ReactionSystem Given a *nonlinear* [`ReactionSystem`](https://catalyst.sciml.ai/stable/api/catalyst_api/#ModelingToolkit.ReactionSystem) and an equivalent *linear* `ReactionSystem`, perform the Linear Mapping Approximation (LMA) @@ -21,12 +21,12 @@ Notes: - By default the moment equations will be generated up to the order determined by the degree of nonlinearity of the nonlinear system's reactions. However, if higher order moment information is required, the optional `m_order` argument may be provided to increase the expansion order manually. -- `combinatoric_ratelaw=true` uses binomials in calculating the propensity functions +- `combinatoric_ratelaws=true` uses binomials in calculating the propensity functions of a `ReactionSystem`, see the notes for [`ModelingToolkit.jumpratelaw`] (https://mtk.sciml.ai/stable/systems/ReactionSystem/#ModelingToolkit.jumpratelaw). """ function linear_mapping_approximation(rn_nonlinear::T, rn_linear::T, binary_vars::Array{Int,1}=Int[], m_order::Int=0; - combinatoric_ratelaw = true) where T <: ReactionSystem + combinatoric_ratelaws = true) where T <: ReactionSystem # check that all necessary information is provided and that LMA is applicable on the given nonlinear system @@ -47,11 +47,11 @@ function linear_mapping_approximation(rn_nonlinear::T, rn_linear::T, binary_vars # apply LMA: substitute reaction parameters in the linear network moment equations to reflect # the nonlinear interactions according to the LMA methodology - term_factors, term_powers, poly_order = polynomial_propensities(propensities(rn_nonlinear; combinatoric_ratelaw)[nonlinear_rs_inds], + term_factors, term_powers, poly_order = polynomial_propensities(propensities(rn_nonlinear; combinatoric_ratelaws)[nonlinear_rs_inds], get_iv(rn_nonlinear), speciesmap(rn_nonlinear)) order = iszero(m_order) ? poly_order : max(poly_order, m_order) - try sys = generate_raw_moment_eqs(rn_linear, order; combinatoric_ratelaw, smap=speciesmap(rn_nonlinear)) + try sys = generate_raw_moment_eqs(rn_linear, order; combinatoric_ratelaws, smap=speciesmap(rn_nonlinear)) catch e error("LMA cannot handle reactions with non-polynomial rates\n $e") end diff --git a/src/raw_moment_equations.jl b/src/raw_moment_equations.jl index d43ab45..e43972a 100644 --- a/src/raw_moment_equations.jl +++ b/src/raw_moment_equations.jl @@ -1,6 +1,6 @@ """ generate_raw_moment_eqs(rn::ReactionSystem, m_order::Int; - langevin::Bool=false, combinatoric_ratelaw::Bool=true, smap=speciesmap(rn)) + langevin::Bool=false, combinatoric_ratelaws::Bool=true, smap=speciesmap(rn)) Given a [`ReactionSystem`](https://catalyst.sciml.ai/stable/api/catalyst_api/#ModelingToolkit.ReactionSystem) return the [`RawMomentEquations`](@ref) of the system generated up to `m_order`. @@ -13,7 +13,7 @@ Notes: - if `langevin=true`, instead of the Chemical Master Equation the Chemical Langevin Equation (diffusion approximation) is considered, and the moment equations are constructed from the corresponding SDE formulation. -- `combinatoric_ratelaw=true` uses binomials in calculating the propensity functions +- `combinatoric_ratelaws=true` uses binomials in calculating the propensity functions of a `ReactionSystem`, see the notes for [`ModelingToolkit.jumpratelaw`] (https://mtk.sciml.ai/stable/systems/ReactionSystem/#ModelingToolkit.jumpratelaw). *Note* that this field is irrelevant using `ReactionSystemMod` as then the @@ -23,12 +23,12 @@ Notes: accessible with [`Catalyst.speciesmap`](https://catalyst.sciml.ai/stable/api/catalyst_api/#Catalyst.speciesmap). """ function generate_raw_moment_eqs(rn::ReactionSystem, m_order::Int; - langevin::Bool=false, combinatoric_ratelaw::Bool=true, smap=speciesmap(rn)) + langevin::Bool=false, combinatoric_ratelaws::Bool=true, smap=speciesmap(rn)) iv = get_iv(rn) N = numspecies(rn) S = get_stoichiometry(rn, smap) - a = propensities(rn; combinatoric_ratelaw) + a = propensities(rn; combinatoric_ratelaws) a = div_to_pow.(a) # NOTE: more stable at the moment than dealing with symbolic fractions if langevin diff --git a/src/reaction_systems.jl b/src/reaction_systems.jl index d7fc320..11d2de2 100644 --- a/src/reaction_systems.jl +++ b/src/reaction_systems.jl @@ -1,15 +1,15 @@ """ - propensities(rn::Union{ReactionSystem, ReactionSystemMod}; combinatoric_ratelaw=true) + propensities(rn::Union{ReactionSystem, ReactionSystemMod}; combinatoric_ratelaws=true) Return a vector of propensity functions of all reactions in the given [`ReactionSystem`] (https://catalyst.sciml.ai/stable/api/catalyst_api/#ModelingToolkit.ReactionSystem). Notes: -- `combinatoric_ratelaw=true` uses binomials in calculating the propensity functions +- `combinatoric_ratelaws=true` uses binomials in calculating the propensity functions of a `ReactionSystem`, see the notes for [`ModelingToolkit.jumpratelaw`] (https://mtk.sciml.ai/stable/systems/ReactionSystem/#ModelingToolkit.jumpratelaw). """ -function propensities(rn::ReactionSystem; combinatoric_ratelaw=true) - simplify.(jumpratelaw.(reactions(rn); combinatoric_ratelaw)) +function propensities(rn::ReactionSystem; combinatoric_ratelaws::Bool=true) + simplify.(jumpratelaw.(reactions(rn); combinatoric_ratelaw=combinatoric_ratelaws)) end """ diff --git a/test/bernoulli_closure_tests.jl b/test/bernoulli_closure_tests.jl index 13e2d34..ac75972 100644 --- a/test/bernoulli_closure_tests.jl +++ b/test/bernoulli_closure_tests.jl @@ -19,7 +19,7 @@ binary_vars = [1] μ = define_μ(2,4) M = define_M(2,4) -sys = generate_central_moment_eqs(rn, 2, combinatoric_ratelaw=false) +sys = generate_central_moment_eqs(rn, 2, combinatoric_ratelaws=false) closed_eqs = moment_closure(sys, "zero", binary_vars) @test length(closed_eqs.odes.states) == 4 && isequal(closed_eqs.closure[M[0,4]], 0) diff --git a/test/closure_tests.jl b/test/closure_tests.jl index 2b9666d..432d0b9 100644 --- a/test/closure_tests.jl +++ b/test/closure_tests.jl @@ -16,7 +16,7 @@ M = define_M(2,4) # --- Test closures on central moment equations --- -sys = generate_central_moment_eqs(rn, 2, 4, combinatoric_ratelaw=false) +sys = generate_central_moment_eqs(rn, 2, 4, combinatoric_ratelaws=false) expr1 = sys.odes.eqs[1].rhs closed_eqs = moment_closure(sys, "zero") @test closed_eqs.closure[M[2,2]] == 0 @@ -56,7 +56,7 @@ expr2 = μ[0,1]^4*(M[0,2]+μ[0,1]^2)^-6*(M[0,3]+μ[0,1]^3+3*M[0,2]*μ[0,1])^4 - # --- Test closures on raw moment equations --- -sys = generate_raw_moment_eqs(rn, 2, combinatoric_ratelaw=false) +sys = generate_raw_moment_eqs(rn, 2, combinatoric_ratelaws=false) closed_eqs = moment_closure(sys, "zero") expr1 = closed_eqs.closure[μ[3,0]] expr2 = -2*μ[1,0]^3 + 3*μ[1,0]*μ[2,0] diff --git a/test/latexify_tests.jl b/test/latexify_tests.jl index 8a6480d..002892f 100644 --- a/test/latexify_tests.jl +++ b/test/latexify_tests.jl @@ -16,7 +16,7 @@ rn = @reaction_network begin end binary_vars = [1] -raw_eqs = generate_raw_moment_eqs(rn, 2, combinatoric_ratelaw=false) +raw_eqs = generate_raw_moment_eqs(rn, 2, combinatoric_ratelaws=false) clean_eqs = bernoulli_moment_eqs(raw_eqs, binary_vars) closed_raw_eqs = moment_closure(raw_eqs, "conditional gaussian", binary_vars) diff --git a/test/lma_tests.jl b/test/lma_tests.jl index 5399c64..1aa27e8 100644 --- a/test/lma_tests.jl +++ b/test/lma_tests.jl @@ -25,7 +25,7 @@ rn_linear = @reaction_network begin end binary_vars = [speciesmap(rn_nonlinear)[g]] -LMA_eqs, _ = linear_mapping_approximation(rn_nonlinear, rn_linear, binary_vars, combinatoric_ratelaw=false) +LMA_eqs, _ = linear_mapping_approximation(rn_nonlinear, rn_linear, binary_vars, combinatoric_ratelaws=false) μ = define_μ(2,3) expr1 = LMA_eqs.odes.eqs[1].rhs @@ -48,7 +48,7 @@ rn_nonlinear = @reaction_network begin end binary_vars = [speciesmap(rn_nonlinear)[g]] -_, effective_params = linear_mapping_approximation(rn_nonlinear, rn_linear, binary_vars, combinatoric_ratelaw=false) +_, effective_params = linear_mapping_approximation(rn_nonlinear, rn_linear, binary_vars, combinatoric_ratelaws=false) expr1 = effective_params[value(σ_b_LMA)] expr2 = (σ_b*μ[1,2] - σ_b*μ[1,1]) * μ[1,0]^-1 @test isequal(expr1, value(expr2)) \ No newline at end of file diff --git a/test/moment_equations_SDE_tests.jl b/test/moment_equations_SDE_tests.jl index 111ccbe..a1a85b2 100644 --- a/test/moment_equations_SDE_tests.jl +++ b/test/moment_equations_SDE_tests.jl @@ -33,7 +33,7 @@ end k1, k2 = parameters(schloegl) for order in 2:6 - schloegl_moments = generate_raw_moment_eqs(schloegl, order; langevin = true, combinatoric_ratelaw = false) + schloegl_moments = generate_raw_moment_eqs(schloegl, order; langevin = true, combinatoric_ratelaws = false) local μ = schloegl_moments.μ rhs(j) = j > 1 ? expand( j * ( k1*μ[(j-1,)] - k2*μ[(j,)] + (μ[(j+1,)] - μ[(j,)]) - (μ[(j+2,)] - 3*μ[(j+1,)] + 2*μ[(j,)]) ) + j*(j-1)/2 * ( k1*μ[(j-2,)] + k2*μ[(j-1,)] + (μ[(j,)] - μ[(j-1,)]) + (μ[(j+1,)] - 3*μ[(j,)] + 2*μ[(j-1,)]) ) ) : @@ -43,7 +43,7 @@ for order in 2:6 end μ = define_μ(1, 4); M = define_M(1, 4) -schloegl_moments = generate_central_moment_eqs(schloegl, 2, 4, langevin=true, combinatoric_ratelaw = false) +schloegl_moments = generate_central_moment_eqs(schloegl, 2, 4, langevin=true, combinatoric_ratelaws = false) expr = k1 + 4*μ[(1,)]^2 + 4*M[(2,)] - μ[(1,)]^3 - M[(3,)] - 3*μ[(1,)] - k2*μ[(1,)] - 3*M[(2,)]*μ[(1,)] @test isequal(schloegl_moments.odes.eqs[1].rhs, expr) expr = k1 + 9*M[(3,)] + k2*μ[(1,)] + μ[(1,)]^3 + 19*M[(2,)]*μ[(1,)] + μ[(1,)] - 2*μ[(1,)]^2 - @@ -61,7 +61,7 @@ end c1, c2, c3, c4 = parameters(rn) for order in 2:10 - rn_moments = generate_raw_moment_eqs(rn, order; langevin = true, combinatoric_ratelaw = false) + rn_moments = generate_raw_moment_eqs(rn, order; langevin = true, combinatoric_ratelaws = false) local μ = rn_moments.μ rhs(i,j) = (i >= 1 ? i*( c1*(μ[(i+1,j+1)]-μ[(i,j+1)]) - (c2 + c4)*μ[(i,j)] + c3*μ[(i-1,j)] ) : 0) + (j >= 1 ? j*( -c1*(μ[(i+2,j)]-μ[(i+1,j)]) + c2*μ[(i+1,j-1)] ) : 0) + @@ -75,8 +75,8 @@ for order in 2:10 @test isequal(analytic_moment_eqs, rn_moments.odes.eqs) end -raw_eqs = generate_raw_moment_eqs(rn, 2, langevin = true, combinatoric_ratelaw=false) -central_eqs = generate_central_moment_eqs(rn, 2, langevin=true, combinatoric_ratelaw=false) +raw_eqs = generate_raw_moment_eqs(rn, 2, langevin = true, combinatoric_ratelaws=false) +central_eqs = generate_central_moment_eqs(rn, 2, langevin=true, combinatoric_ratelaws=false) μ = define_μ(2, 4); M = define_μ(2, 4) central_to_raw = central_to_raw_moments(2, 4) subdict = Dict( μ[iter] => central_to_raw[iter] for iter in filter(x -> sum(x) > 1, raw_eqs.iter_all)) diff --git a/test/moment_equations_tests.jl b/test/moment_equations_tests.jl index 9f40e78..bc7ec28 100644 --- a/test/moment_equations_tests.jl +++ b/test/moment_equations_tests.jl @@ -32,7 +32,7 @@ rn = @reaction_network begin (Ω*c₃, c₄), 0 ↔ X end -sys = generate_central_moment_eqs(rn, 2, 4, combinatoric_ratelaw=false) +sys = generate_central_moment_eqs(rn, 2, 4, combinatoric_ratelaws=false) expr1 = sys.odes.eqs[2].rhs @test isequal(MomentClosure.Differential(t)(sys.μ[1,0]), sys.odes.eqs[1].lhs) μ = sys.μ @@ -45,12 +45,12 @@ expr2 = simplify(value.(expr2)) @test isequal(expr1, expr2) -sys = generate_central_moment_eqs(rn, 2, combinatoric_ratelaw=false) +sys = generate_central_moment_eqs(rn, 2, combinatoric_ratelaws=false) expr1 = sys.odes.eqs[2].rhs @test isequal(MomentClosure.Differential(t)(sys.μ[1,0]), sys.odes.eqs[1].lhs) @test isequal(expr1, expr2) -sys = generate_raw_moment_eqs(rn, 2, combinatoric_ratelaw=false) +sys = generate_raw_moment_eqs(rn, 2, combinatoric_ratelaws=false) μ = sys.μ expr1 = sys.odes.eqs[4].rhs expr2 = c₂*μ[2,0] + c₁*μ[1,1]*Ω^-2 - c₁*μ[3,1]*Ω^-2 -c₂*(μ[1,0] + μ[1,1]) + diff --git a/test/reaction_system_tests.jl b/test/reaction_system_tests.jl index 489d0fe..0c8fd2f 100644 --- a/test/reaction_system_tests.jl +++ b/test/reaction_system_tests.jl @@ -7,7 +7,7 @@ using Catalyst S_mat = [ 1 -1 1 -1; -1 1 0 0] -# NOTE: only holds if combinatoric_ratelaw=false +# NOTE: only holds if combinatoric_ratelaws=false a = [c₁*X*Y*(X-1)/Ω^2, c₂*X, c₃*Ω, c₄*X] rn = @reaction_network begin @@ -19,4 +19,4 @@ end smap = speciesmap(rn) @test isequal(get_stoichiometry(rn, smap), S_mat) -@test isequal(propensities(rn, combinatoric_ratelaw=false), a) +@test isequal(propensities(rn, combinatoric_ratelaws=false), a) From 14dd10f4676b84f3fd0aad277c744292f2607c00 Mon Sep 17 00:00:00 2001 From: Augustinas Sukys Date: Fri, 5 May 2023 00:04:07 +0100 Subject: [PATCH 10/16] remove redundant code not needed after SymbolicUtils update --- src/central_moment_equations.jl | 1 - src/raw_moment_equations.jl | 1 - src/symbolic.jl | 11 ----------- 3 files changed, 13 deletions(-) diff --git a/src/central_moment_equations.jl b/src/central_moment_equations.jl index 0511274..cd4b38e 100644 --- a/src/central_moment_equations.jl +++ b/src/central_moment_equations.jl @@ -48,7 +48,6 @@ function generate_central_moment_eqs(rn::ReactionSystem, m_order::Int, q_order:: R = numreactions(rn) # no. of reactions in the network iv = get_iv(rn) # independent variable (usually time) a = propensities(rn; combinatoric_ratelaws) # propensity functions of all reactions in the network - a = div_to_pow.(a) # convert Div to Pow (more stable at the moment than dealing with symbolic fractions) S = get_stoichiometry(rn, smap) # net stoichiometric matrix if langevin diff --git a/src/raw_moment_equations.jl b/src/raw_moment_equations.jl index e43972a..63dd322 100644 --- a/src/raw_moment_equations.jl +++ b/src/raw_moment_equations.jl @@ -29,7 +29,6 @@ function generate_raw_moment_eqs(rn::ReactionSystem, m_order::Int; N = numspecies(rn) S = get_stoichiometry(rn, smap) a = propensities(rn; combinatoric_ratelaws) - a = div_to_pow.(a) # NOTE: more stable at the moment than dealing with symbolic fractions if langevin drift = S*a diff --git a/src/symbolic.jl b/src/symbolic.jl index fe88e9f..dadb82b 100644 --- a/src/symbolic.jl +++ b/src/symbolic.jl @@ -102,17 +102,6 @@ function extract_variables(eqs::Array{Equation, 1}, μ, M=[]) end -# a simple hack to convert SymbolicUtils Div (fractions) to Pow (multiplication) -# more stable in symbolic manipulations that we perform -function div_to_pow(expr) - if isdiv(expr) - args = arguments(expr) - args[1] * args[2]^-1 - else - expr - end -end - ## Set of functions to deconstruct polynomial propensities ## #= From f46e7e162042689f9cf42ca0afd65d2a85b4a639 Mon Sep 17 00:00:00 2001 From: Augustinas Sukys Date: Fri, 5 May 2023 00:19:43 +0100 Subject: [PATCH 11/16] update github actions --- .github/workflows/Documentation.yml | 4 ++-- .github/workflows/ci.yml | 6 +++--- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/.github/workflows/Documentation.yml b/.github/workflows/Documentation.yml index 01e71aa..245daf2 100644 --- a/.github/workflows/Documentation.yml +++ b/.github/workflows/Documentation.yml @@ -11,7 +11,7 @@ jobs: build: runs-on: ubuntu-latest steps: - - uses: actions/checkout@v2 + - uses: actions/checkout@v3 - uses: julia-actions/setup-julia@latest with: version: '1' @@ -20,5 +20,5 @@ jobs: - name: Build and deploy env: GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }} # For authentication with GitHub Actions token - #DOCUMENTER_KEY: ${{ secrets.DOCUMENTER_KEY }} # For authentication with SSH deploy key + DOCUMENTER_KEY: ${{ secrets.DOCUMENTER_KEY }} # For authentication with SSH deploy key run: julia --project=docs/ docs/make.jl diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 0308eb3..7544a1b 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -16,12 +16,12 @@ jobs: arch: - x64 steps: - - uses: actions/checkout@v2 + - uses: actions/checkout@v3 - uses: julia-actions/setup-julia@v1 with: version: ${{ matrix.version }} arch: ${{ matrix.arch }} - - uses: actions/cache@v1 + - uses: actions/cache@v3 env: cache-name: cache-artifacts with: @@ -34,6 +34,6 @@ jobs: - uses: julia-actions/julia-buildpkg@v1 - uses: julia-actions/julia-runtest@v1 - uses: julia-actions/julia-processcoverage@v1 - - uses: codecov/codecov-action@v1 + - uses: codecov/codecov-action@v3 with: file: lcov.info From 653268086de951e7036ace8503031cd7829574df Mon Sep 17 00:00:00 2001 From: Augustinas Sukys Date: Fri, 5 May 2023 00:20:15 +0100 Subject: [PATCH 12/16] enable dependabot for github actions --- .github/dependabot.yml | 7 +++++++ 1 file changed, 7 insertions(+) create mode 100644 .github/dependabot.yml diff --git a/.github/dependabot.yml b/.github/dependabot.yml new file mode 100644 index 0000000..ff6499d --- /dev/null +++ b/.github/dependabot.yml @@ -0,0 +1,7 @@ +# https://docs.github.com/github/administering-a-repository/configuration-options-for-dependency-updates +version: 2 +updates: + - package-ecosystem: "github-actions" + directory: "/" # Location of package manifests + schedule: + interval: "weekly" \ No newline at end of file From 92393d99e12557599f123422f4a70f0a054841e6 Mon Sep 17 00:00:00 2001 From: Augustinas Sukys Date: Fri, 5 May 2023 00:45:13 +0100 Subject: [PATCH 13/16] bump Latexify compat --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 1745673..37c9a79 100644 --- a/Project.toml +++ b/Project.toml @@ -25,7 +25,7 @@ Cumulants = "1" DataStructures = "0.18" Distributions = "0.25" DocStringExtensions = "0.8, 0.9" -Latexify = "0.15.14" +Latexify = "0.15.14, 0.16" ModelingToolkit = "8" SciMLBase = "1.8.3" SymbolicUtils = "1" From 57e83274cae779a4c0c573ca7e05f0f4066006a7 Mon Sep 17 00:00:00 2001 From: Augustinas Sukys Date: Sat, 6 May 2023 01:32:46 +0100 Subject: [PATCH 14/16] fine-tune symbolics --- src/central_moment_equations.jl | 6 ++---- src/closure_methods/closure.jl | 3 +-- .../conditional_derivative_matching.jl | 4 ++-- src/closure_methods/derivative_matching.jl | 2 +- src/closure_methods/gamma_closure.jl | 6 +++--- src/closure_methods/log_normal_closure.jl | 6 +++--- src/moment_convert.jl | 4 ++-- src/raw_moment_equations.jl | 2 +- src/symbolic.jl | 12 +++++------- src/utils.jl | 4 ++-- test/closure_tests.jl | 5 +++-- test/conditional_closure_tests.jl | 4 ++-- test/moment_equations_tests.jl | 7 +++---- 13 files changed, 30 insertions(+), 35 deletions(-) diff --git a/src/central_moment_equations.jl b/src/central_moment_equations.jl index cd4b38e..7567445 100644 --- a/src/central_moment_equations.jl +++ b/src/central_moment_equations.jl @@ -116,9 +116,7 @@ function generate_central_moment_eqs(rn::ReactionSystem, m_order::Int, q_order:: suma = 0 for j in iter_all suma = suma + Da[r][j]*M[j]*1//fact(j) - #suma = simplify(suma) # here // gives fractions (otherwise it's a float number) - # only issue is that it does not latexify properly end for i in 1:N @@ -127,7 +125,7 @@ function generate_central_moment_eqs(rn::ReactionSystem, m_order::Int, q_order:: else du[i] = S[i, r]*suma + du[i] end - du[i] = simplify(du[i]) + du[i] = expand_expr(du[i]) end end @@ -157,7 +155,7 @@ function generate_central_moment_eqs(rn::ReactionSystem, m_order::Int, q_order:: dM[i] -= i[j]*du[j]*M[i.-iter_1[j]] end end - dM[i] = simplify(dM[i]) + dM[i] = expand_expr(dM[i]) end D = Differential(iv) diff --git a/src/closure_methods/closure.jl b/src/closure_methods/closure.jl index ceefc86..4d37ca4 100644 --- a/src/closure_methods/closure.jl +++ b/src/closure_methods/closure.jl @@ -1,7 +1,6 @@ function close_eqs(sys::MomentEquations, closure_exp::OrderedDict, closure::OrderedDict, polynorm::Bool) - # TODO: improve performance closed_eqs = Equation[] for eq in sys.odes.eqs @@ -10,7 +9,7 @@ function close_eqs(sys::MomentEquations, closure_exp::OrderedDict, if polynorm # depending on the functional form closed_rhs = expand(closed_rhs) else - closed_rhs = expand_expr(closed_rhs) + closed_rhs = expand_mod(closed_rhs) end push!(closed_eqs, Equation(eq.lhs, closed_rhs)) end diff --git a/src/closure_methods/conditional_derivative_matching.jl b/src/closure_methods/conditional_derivative_matching.jl index 0a97673..09f3b9a 100644 --- a/src/closure_methods/conditional_derivative_matching.jl +++ b/src/closure_methods/conditional_derivative_matching.jl @@ -163,8 +163,8 @@ function conditional_derivative_matching(sys::MomentEquations, raw_to_central_exp = raw_to_central_moments(N, sys.q_order, μ_M_exp, bernoulli=true) raw_to_central = raw_to_central_moments(N, sys.q_order, μ_M, bernoulli=true) for i in sys.iter_q - closure_exp[sys.M[i]] = raw_to_central_exp[i] - closure[sys.M[i]] = raw_to_central[i] + closure_exp[sys.M[i]] = expand_mod(raw_to_central_exp[i]) + closure[sys.M[i]] = expand_mod(raw_to_central[i]) end else diff --git a/src/closure_methods/derivative_matching.jl b/src/closure_methods/derivative_matching.jl index 15d6152..0f1c460 100644 --- a/src/closure_methods/derivative_matching.jl +++ b/src/closure_methods/derivative_matching.jl @@ -5,7 +5,7 @@ function multi_binomial(a, b) if length(a) != length(b) throw(ErrorException("vector lengths are different for some reason")) else - C = prod([binomial(a[i], b[i]) for i in 1:length(a)]) + C = prod([binomial(a[i], b[i]) for i in eachindex(a)]) return C end end diff --git a/src/closure_methods/gamma_closure.jl b/src/closure_methods/gamma_closure.jl index b2d9d09..545fc44 100644 --- a/src/closure_methods/gamma_closure.jl +++ b/src/closure_methods/gamma_closure.jl @@ -185,12 +185,12 @@ function gamma_closure(sys::MomentEquations, binary_vars::Array{Int,1}=Int[]) for i in sys.iter_q closure_exp[sys.μ[i]] = substitute(μ[i], M_to_μ) expr = substitute(closure[sys.μ[i]], M_to_μ) - closure[sys.μ[i]] = simplify(expr) + closure[sys.μ[i]] = simplify(expr, simplify_fractions=false) end end if isempty(binary_vars) - return close_eqs(sys, closure_exp, closure, true) + return close_eqs(sys, closure_exp, closure, false) else # have to perform the bernoulli simplify in the end for gamma closure # as otherwise some symbolic terms do not cancel out and break it @@ -204,7 +204,7 @@ function gamma_closure(sys::MomentEquations, binary_vars::Array{Int,1}=Int[]) closed_rhs = substitute(eq.rhs, closure_exp) closed_rhs = expand(closed_rhs) closed_rhs = substitute(closed_rhs, iter_sub) - closed_rhs = simplify(closed_rhs) + closed_rhs = simplify(closed_rhs, simplify_fractions=false) push!(closed_eqs, Equation(eq.lhs, closed_rhs)) diff --git a/src/closure_methods/log_normal_closure.jl b/src/closure_methods/log_normal_closure.jl index 0174eb3..e1ee438 100644 --- a/src/closure_methods/log_normal_closure.jl +++ b/src/closure_methods/log_normal_closure.jl @@ -1,7 +1,7 @@ function log_normal_closure(sys::MomentEquations, binary_vars::Array{Int,1}=Int[]) closure = OrderedDict() - closure_exp = OrderedDict() # here it does not play a role + closure_exp = OrderedDict() N = sys.N if isempty(binary_vars) @@ -76,9 +76,9 @@ function log_normal_closure(sys::MomentEquations, binary_vars::Array{Int,1}=Int[ central_to_raw = central_to_raw_moments(N, sys.q_order) closure_M = OrderedDict() for i in sys.iter_q - closure_exp[M[i]] = raw_to_central[i] + closure_exp[M[i]] = expand_mod(raw_to_central[i]) closure_M[M[i]] = closure[μ_symbolic[i]]-(central_to_raw[i]-M[i]) - closure_M[M[i]] = simplify(closure_M[M[i]]) + closure_M[M[i]] = expand_mod(closure_M[M[i]]) end closure = closure_M else diff --git a/src/moment_convert.jl b/src/moment_convert.jl index 6016606..dc4035d 100644 --- a/src/moment_convert.jl +++ b/src/moment_convert.jl @@ -183,7 +183,7 @@ function raw_to_central_moments(N::Int, order::Int, μ=nothing; bernoulli=false) end suma += term end - raw_to_central[i] = simplify(suma) + raw_to_central[i] = simplify(suma, simplify_fractions=false) end raw_to_central @@ -217,7 +217,7 @@ function central_to_raw_moments(N::Int, order::Int) end suma += term end - central_to_raw[i] = simplify(suma) + central_to_raw[i] = simplify(suma, simplify_fractions=false) end central_to_raw diff --git a/src/raw_moment_equations.jl b/src/raw_moment_equations.jl index 63dd322..2742dfa 100644 --- a/src/raw_moment_equations.jl +++ b/src/raw_moment_equations.jl @@ -71,7 +71,7 @@ function generate_raw_moment_eqs(rn::ReactionSystem, m_order::Int; dμ[i] += factor_j * suma end end - dμ[i] = simplify(dμ[i]) + dμ[i] = expand_expr(dμ[i]) end D = Differential(iv) diff --git a/src/symbolic.jl b/src/symbolic.jl index dadb82b..270089c 100644 --- a/src/symbolic.jl +++ b/src/symbolic.jl @@ -19,13 +19,11 @@ end # Trim a string of form "(a, b, c, d, ...)" to "abcd..." trim_key(expr) = filter(x -> !(isspace(x) || x == ')' || x== '(' || x==','), string(expr)) -# Expand a symbolic expression (no binomial expansion) -expansion_rule_mod = @acrule ~x * +(~~ys) => sum(map(y-> ~x * y, ~~ys)) # apply distribution law -expand_mod = Fixpoint(Prewalk(PassThrough(expansion_rule_mod))) # distributes terms until no longer possible -flatten_rule_mod = @rule(~x::isnotflat(+) => flatten_term(+, ~x)) # -flatten_mod = Fixpoint(PassThrough(flatten_rule_mod)) # -expand_expr = Fixpoint(PassThrough(Chain([expand_mod, flatten_mod]))) # apply flatten and distribution until no longer possible -expand_div = PassThrough(@acrule( +(~~xs) / ~a => sum(map(x -> x / ~a, ~~xs)))) +# Symbolic rules to expand expressions avoiding binomial expansion and fraction simplification +expansion_rule_mod = @acrule ~x * +(~~ys) => sum(map(y-> ~x * y, ~~ys)) +expand_expr = Fixpoint(Prewalk(PassThrough(expansion_rule_mod))) +expand_div = Fixpoint(Prewalk(PassThrough(@rule( +(~~xs) / ~a => sum(map(x -> x / ~a, ~~xs) ))))) +expand_mod = Fixpoint(Chain([expand_div, expand_expr])) function define_μ(iter::AbstractVector, iv::BasicSymbolic) diff --git a/src/utils.jl b/src/utils.jl index b915655..99053fc 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -307,7 +307,7 @@ function format_moment_eqs(eqs::MomentEquations) for i in 1:size(odes)[1] key = vars[i] eq = odes[i].rhs - eq = expand_div(eq) + eq = expand_mod(eq) expr = "d"*string(key)*"/dt = "*string(eq) expr = replace(expr, "(t)"=>"") expr = replace(expr, ".0"=>"") @@ -342,7 +342,7 @@ function format_closure(eqs::ClosedMomentEquations; format_all::Bool=false) for i in iter eq = closure[i] - eq = expand_div(eq) + eq = expand_mod(eq) expr = string(i)*" = "*string(eq) expr = replace(expr, "(t)"=>"") expr = replace(expr, ".0"=>"") diff --git a/test/closure_tests.jl b/test/closure_tests.jl index 432d0b9..0f014b6 100644 --- a/test/closure_tests.jl +++ b/test/closure_tests.jl @@ -1,6 +1,7 @@ using MomentClosure using MomentClosure: define_M, define_μ using Symbolics: value, expand +using SymbolicUtils: Fixpoint using Test using Catalyst @@ -37,7 +38,7 @@ closed_eqs= moment_closure(sys, "log-normal") expr1 = closed_eqs.closure[M[1,2]] expr2 = μ[1,0]*μ[0,1]^2*(1.0+M[0,2]*μ[0,1]^-2)*(1.0 + M[1,1]*(μ[0,1]^-1)*(μ[1,0]^-1))^2 - M[0,2]*μ[1,0] - μ[1,0]*μ[0,1]^2 - 2*M[1,1]*μ[0,1] -@test isequal(expr1, simplify(expand(expr2))) +@test isequal(Fixpoint(simplify)(expr1), Fixpoint(simplify)(expr2)) closed_eqs = moment_closure(sys, "poisson") @test isequal(closed_eqs.closure[M[3,0]], μ[1,0]) @@ -52,7 +53,7 @@ closed_eqs = moment_closure(sys, "derivative matching") expr1 = closed_eqs.closure[sys.M[0,4]] expr2 = μ[0,1]^4*(M[0,2]+μ[0,1]^2)^-6*(M[0,3]+μ[0,1]^3+3*M[0,2]*μ[0,1])^4 - μ[0,1]^4 - 6*M[0,2]*μ[0,1]^2 - 4*M[0,3]*μ[0,1] -@test isequal(expand(expr1), expand(simplify(expr2))) +@test isequal(simplify(expr1), simplify(expr2)) # --- Test closures on raw moment equations --- diff --git a/test/conditional_closure_tests.jl b/test/conditional_closure_tests.jl index 717381a..3134689 100644 --- a/test/conditional_closure_tests.jl +++ b/test/conditional_closure_tests.jl @@ -70,8 +70,8 @@ expr1 = expand(closed_eqs.closure[M[0,5]]) expr2 = 10*M[0,2]*M[0,3] + 5*M[0,4]*μ[0,1] - 15*μ[0,1]*M[0,2]^2 @test isequal(expr1, expr2) -closed_eqs = moment_closure(sys, "conditional derivative matching", binary_vars) +@time closed_eqs = moment_closure(sys, "conditional derivative matching", binary_vars) expr1 = closed_eqs.closure[M[1,3]] expr2 = μ[1,0]*(M[1,1]+μ[0,1]*μ[1,0])^-3*(M[1,2]+M[0,2]*μ[1,0]+μ[1,0]*μ[0,1]^2+2*M[1,1]*μ[0,1])^3 - M[0,3]*μ[1,0] - 3*M[1,1]*μ[0,1]^2 - 3*M[1,2]*μ[0,1] - μ[1,0]*μ[0,1]^3- 3*M[0,2]*μ[0,1]*μ[1,0] -@test isequal(expr1, simplify(expr2)) \ No newline at end of file +@test isequal(simplify(expr1), simplify(expr2)) \ No newline at end of file diff --git a/test/moment_equations_tests.jl b/test/moment_equations_tests.jl index bc7ec28..cee2033 100644 --- a/test/moment_equations_tests.jl +++ b/test/moment_equations_tests.jl @@ -41,14 +41,13 @@ M = sys.M expr2 = c₂*μ[1,0] + c₁*M[1,1]*(Ω^-2) + c₁*μ[0,1]*μ[1,0]*Ω^-2 - c₁*M[2,1]*Ω^-2 - 2*c₁*M[1,1]*μ[1,0]*Ω^-2 - c₁*M[2,0]*μ[0,1]*Ω^-2 - c₁*μ[0,1]*Ω^-2*μ[1,0]^2 expr2 = simplify(value.(expr2)) -#@test isequal(expand(expr1), expr2) -@test isequal(expr1, expr2) +@test isequal(simplify(expr1), expr2) sys = generate_central_moment_eqs(rn, 2, combinatoric_ratelaws=false) expr1 = sys.odes.eqs[2].rhs @test isequal(MomentClosure.Differential(t)(sys.μ[1,0]), sys.odes.eqs[1].lhs) -@test isequal(expr1, expr2) +@test isequal(simplify(expr1), expr2) sys = generate_raw_moment_eqs(rn, 2, combinatoric_ratelaws=false) μ = sys.μ @@ -56,7 +55,7 @@ expr1 = sys.odes.eqs[4].rhs expr2 = c₂*μ[2,0] + c₁*μ[1,1]*Ω^-2 - c₁*μ[3,1]*Ω^-2 -c₂*(μ[1,0] + μ[1,1]) + c₃*μ[0,1]*Ω - c₄*μ[1,1] - c₁*μ[1,2]*Ω^-2 + c₁*μ[2,2]*Ω^-2 expr2 = simplify(value.(expr2)) -@test isequal(expr1, expr2) +@test isequal(simplify(expr1), expr2) # corner case - a linear propensity with rate coefficient 1 rn = @reaction_network begin From 50e428ea5ea0039670c3b691bee36478fea9cc2a Mon Sep 17 00:00:00 2001 From: Augustinas Sukys Date: Sat, 6 May 2023 01:34:06 +0100 Subject: [PATCH 15/16] update docs --- .../assets/derivative_matching_cumulant.svg | 169 +++++------------- docs/src/tutorials/LMA_example.md | 19 +- docs/src/tutorials/P53_system_example.md | 31 ++-- docs/src/tutorials/common_issues.md | 29 +-- .../tutorials/derivative_matching_example.md | 5 +- ...eometric_reactions+conditional_closures.md | 97 +++------- .../src/tutorials/parameter_estimation_SDE.md | 2 +- .../tutorials/time-dependent_propensities.md | 9 +- docs/src/tutorials/using_momentclosure.md | 11 +- 9 files changed, 125 insertions(+), 247 deletions(-) diff --git a/docs/src/assets/derivative_matching_cumulant.svg b/docs/src/assets/derivative_matching_cumulant.svg index 04e74fc..5ed423d 100644 --- a/docs/src/assets/derivative_matching_cumulant.svg +++ b/docs/src/assets/derivative_matching_cumulant.svg @@ -1,141 +1,54 @@ - + - + - + - - + + - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/docs/src/tutorials/LMA_example.md b/docs/src/tutorials/LMA_example.md index 57a0f10..75c4aa4 100644 --- a/docs/src/tutorials/LMA_example.md +++ b/docs/src/tutorials/LMA_example.md @@ -20,20 +20,22 @@ using Catalyst # how the nonlinear reactions are to be transformed using LMA rn_nonlinear = @reaction_network begin + @parameters σ_b σ_u ρ_b ρ_u σ_b, g + p → 0 σ_u*(1-g), 0 ⇒ g + p ρ_u, g → g + p ρ_b*(1-g), 0 ⇒ p 1, p → 0 -end σ_b σ_u ρ_b ρ_u +end rn_linear = @reaction_network begin + @parameters σ_b_LMA σ_u ρ_b ρ_u σ_b_LMA, g → 0 # typing ̄σ_b is not allowed it seems σ_u*(1-g), 0 ⇒ g ρ_u, g → g+p (ρ_b*(1-g)), 0 ⇒ p 1, p → 0 -end σ_b_LMA σ_u ρ_b ρ_u +end ``` We can now apply the LMA to find the effective parameter $\bar{σ}_b$ and generate the corresponding moment equations of the linear GRN using MomentClosure's [`linear_mapping_approximation`](@ref): ```julia @@ -45,7 +47,7 @@ using MomentClosure @variables g(t) binary_vars = [speciesmap(rn_nonlinear)[g]] -LMA_eqs, effective_params = linear_mapping_approximation(rn_nonlinear, rn_linear, binary_vars, combinatoric_ratelaw=false) +LMA_eqs, effective_params = linear_mapping_approximation(rn_nonlinear, rn_linear, binary_vars, combinatoric_ratelaws=false) display(effective_params) ``` ```julia @@ -79,7 +81,7 @@ u₀map = deterministic_IC(u₀, LMA_eqs) oprob_LMA = ODEProblem(LMA_eqs, u₀map, tspan, p) sol_LMA = solve(oprob_LMA, CVODE_BDF(), saveat=dt) -plot(sol_LMA, vars=(0, [2]), label="LMA", ylabel="⟨p⟩", xlabel="time", fmt="svg") +plot(sol_LMA, idxs=[2], label="LMA", ylabel="⟨p⟩", xlabel="time", fmt="svg") ``` ![LMA feedback loop mean protein number](../assets/LMA_feedback_loop_mean_protein_number.svg) @@ -141,11 +143,12 @@ As pointed out by the authors in [1], a more efficient way is to expand the gene However, TaylorSeries only supports elementary function operations at the time and hence evaluating the Kummer's function $M(\cdot,\cdot,\cdot)$ requires some more work (these specialised numerics are readily available in more established scientific computing frameworks such as Mathematica but there's no fun in that). We can extend the TaylorSeries framework by constructing a function `t_pFq` that implements a recurrence relation between the Taylor coefficients for the generalized hypergeometric function `pFq` as defined in [HypergeometricFunctions.jl](https://github.com/JuliaMath/HypergeometricFunctions.jl). This can be done as follows (note that our construction is valid only for a single-variable Taylor series [`Taylor1`](https://juliadiff.org/TaylorSeries.jl/stable/api/#TaylorSeries.Taylor1)): ```julia using TaylorSeries, HypergeometricFunctions +using HypergeometricFunctions: pFqweniger # please let me know if a simpler and more efficient way to do this exists! function t_pFq(α::AbstractVector, β::AbstractVector, a::Taylor1) order = a.order - aux = pFq(α, β, constant_term(a)) + aux = pFqweniger(α, β, constant_term(a)) c = Taylor1(aux, order) iszero(order) && return c @@ -170,13 +173,13 @@ oprob_LMA = remake(oprob_LMA, tspan=tspan) sol_LMA = solve(oprob_LMA, CVODE_BDF(), saveat=dt) # rebuild the symbolic expression for the effective parameter as a function of raw moments -μ_sym = LMA_eqs.odes.states -p_sub = Pair.(LMA_eqs.odes.ps, p) +μ_sym = [LMA_eqs.odes.states...] +p_sub = Dict.(Pair.(LMA_eqs.odes.ps, p)...) avg_σ_b_sym = collect(values(effective_params))[1] fn = build_function(substitute(avg_σ_b_sym, p_sub), μ_sym) avg_σ_b = eval(fn) # evaluate the time-averaged value of the effective parameter -σ_b_avg = sum(avg_σ_b.(sol_LMA[:])) * dt / T +@time σ_b_avg = sum(avg_σ_b.(sol_LMA.u)) * dt / T ``` We proceed with the very last steps of the LMA to obtain the probability distribution: ```julia diff --git a/docs/src/tutorials/P53_system_example.md b/docs/src/tutorials/P53_system_example.md index 120bd69..2d15979 100644 --- a/docs/src/tutorials/P53_system_example.md +++ b/docs/src/tutorials/P53_system_example.md @@ -40,7 +40,7 @@ with parameters We begin by loading all the packages we will need ```julia -using Catalyst, MomentClosure, OrdinaryDiffEq, DiffEqJump, +using Catalyst, MomentClosure, OrdinaryDiffEq, JumpProcesses, DiffEqBase.EnsembleAnalysis, Plots, Plots.PlotMeasures ``` and then build the model using Catalyst and set its parameters as follows: @@ -48,13 +48,14 @@ and then build the model using Catalyst and set its parameters as follows: # → for mass-actions rate # ⇒ for non mass-actions rate rn = @reaction_network begin + @parameters k₁ k₂ k₃ k₄ k₅ k₆ k₇ (k₁), 0 → x (k₂), x → 0 (k₃*x*y/(x+k₇)), x ⇒ 0 (k₄*x), 0 ⇒ y₀ (k₅), y₀ → y (k₆), y → 0 -end k₁ k₂ k₃ k₄ k₅ k₆ k₇ +end # parameters [k₁, k₂, k₃, k₄, k₅, k₆, k₇] p = [90, 0.002, 1.7, 1.1, 0.93, 0.96, 0.01] @@ -106,7 +107,7 @@ h2 = histogram(data[2], normalize=true, xlabel="y₀", ylabel="P(y₀)") h3 = histogram(data[3], normalize=true, xlabel="y", ylabel="P(y)") using Plots.PlotMeasures plot(h1, h2, h3, legend=false, layout=(1,3), size = (1050, 250), - left_margin = 5mm, bottom_margin = 7mm, guidefontsize=10) + left_margin = 5PlotMeasures.mm, bottom_margin = 7PlotMeasures.mm, guidefontsize=10) ``` ![P53-Mdm2 distribution](../assets/p53-Mdm2_distribution.svg) @@ -120,7 +121,7 @@ closures = ["normal", "log-normal", "gamma"] plts = [plot() for i in 1:length(closures)] for q in 3:6 - eqs = generate_central_moment_eqs(rn, 2, q, combinatoric_ratelaw=false) + eqs = generate_central_moment_eqs(rn, 2, q, combinatoric_ratelaws=false) for (closure, plt) in zip(closures, plts) closed_eqs = moment_closure(eqs, closure) @@ -128,7 +129,7 @@ for q in 3:6 oprob = ODEProblem(closed_eqs, u₀map, tspan, p) sol = solve(oprob, Tsit5(), saveat=0.1) - plt = plot!(plt, sol, vars=(0, 1), lw=3, label = "q = "*string(q)) + plt = plot!(plt, sol, idxs=[1], lw=3, label = "q = "*string(q)) end end @@ -139,13 +140,13 @@ end ``` Normal closure: ```julia -plot(plts[1], size=(750, 450), leftmargin=2mm) +plot(plts[1], size=(750, 450), leftmargin=2PlotMeasures.mm) ``` ![P53-Mdm2 normal means 2nd order expansion](../assets/p53-Mdm2_normal_2nd_order.svg) Log-normal closure: ```julia -plot(plts[2], size=(750, 450), leftmargin=2mm) +plot(plts[2], size=(750, 450), leftmargin=2PlotMeasures.mm) ``` ![P53-Mdm2 log-normal means 2nd order expansion](../assets/p53-Mdm2_log-normal_2nd_order.svg) @@ -157,7 +158,7 @@ plot(plts[2], xlims=(0., 50.), lw=3) Gamma closure: ```julia -plot(plts[3], size=(750, 450), leftmargin=2mm) +plot(plts[3], size=(750, 450), leftmargin=2PlotMeasures.mm) ``` ![P53-Mdm2 gamma means 2nd order expansion](../assets/p53-Mdm2_gamma_2nd_order.svg) @@ -168,7 +169,7 @@ We can also plot the variance predictions: # rerunning the same calculations as they are reasonably fast plt = plot() for q in [4,6] - eqs = generate_central_moment_eqs(rn, 2, q, combinatoric_ratelaw=false) + eqs = generate_central_moment_eqs(rn, 2, q, combinatoric_ratelaws=false) for closure in closures closed_eqs = moment_closure(eqs, closure) @@ -177,7 +178,7 @@ for q in [4,6] sol = solve(oprob, Tsit5(), saveat=0.1) # index of M₂₀₀ can be checked with `u₀map` or `closed_eqs.odes.states` - plt = plot!(plt, sol, vars=(0, 4), lw=3, label = closure*" q = "*string(q)) + plt = plot!(plt, sol, idxs=[4], lw=3, label = closure*" q = "*string(q)) end end @@ -200,7 +201,7 @@ q_vals = [4, 6] for (q, plt_m, plt_v) in zip(q_vals, plt_means, plt_vars) - eqs = generate_central_moment_eqs(rn, 3, q, combinatoric_ratelaw=false) + eqs = generate_central_moment_eqs(rn, 3, q, combinatoric_ratelaws=false) for closure in closures closed_eqs = moment_closure(eqs, closure) @@ -209,8 +210,8 @@ for (q, plt_m, plt_v) in zip(q_vals, plt_means, plt_vars) oprob = ODEProblem(closed_eqs, u₀map, tspan, p) sol = solve(oprob, Tsit5(), saveat=0.1) - plt_m = plot!(plt_m, sol, vars=(0, 1), label = closure) - plt_v = plot!(plt_v, sol, vars=(0, 4), label = closure) + plt_m = plot!(plt_m, sol, idxs=[1], label = closure) + plt_v = plot!(plt_v, sol, idxs=[4], label = closure) end @@ -270,7 +271,7 @@ Finally, we can check whether better estimates can be obtained using even higher plt = plot() closures = ["zero", "normal", "log-normal", "gamma"] -eqs = generate_central_moment_eqs(rn, 5, 6, combinatoric_ratelaw=false) +eqs = generate_central_moment_eqs(rn, 5, 6, combinatoric_ratelaws=false) # faster to store than recompute in case we want to try different solvers/params oprobs = Dict() @@ -281,7 +282,7 @@ for closure in closures oprobs[closure] = ODEProblem(closed_eqs, u₀map, tspan, p) sol = solve(oprobs[closure], Tsit5(), saveat=0.1) - plt = plot!(plt, sol, vars=(0, 1), label = closure) + plt = plot!(plt, sol, idxs=[1], label = closure) end plt = plot!(plt, xlabel = "Time [h]", ylabel = "Mean p53 molecule number") diff --git a/docs/src/tutorials/common_issues.md b/docs/src/tutorials/common_issues.md index ef7452e..1811970 100644 --- a/docs/src/tutorials/common_issues.md +++ b/docs/src/tutorials/common_issues.md @@ -7,16 +7,17 @@ We first redefine the system and its parameters for completeness: using MomentClosure, Catalyst, OrdinaryDiffEq, Plots rn = @reaction_network begin + @parameters c₁ c₂ c₃ c₄ Ω (c₁/Ω^2), 2X + Y → 3X (c₂), X → Y (c₃*Ω, c₄), 0 ↔ X -end c₁ c₂ c₃ c₄ Ω +end p = [0.9, 2, 1, 1, 100] u₀ = [1, 1] tspan = (0., 100.) -raw_eqs = generate_raw_moment_eqs(rn, 2, combinatoric_ratelaw=false) +raw_eqs = generate_raw_moment_eqs(rn, 2, combinatoric_ratelaws=false) ``` As we have seen earlier, second-order moment expansion using normal closure approximates the true system dynamics sufficiently accurately but it's interesting to see how other closures compare. Let's try applying zero closure: ```julia @@ -26,7 +27,7 @@ u₀map = deterministic_IC(u₀, closed_raw_eqs) oprob = ODEProblem(closed_raw_eqs, u₀map, tspan, p) sol = solve(oprob, Tsit5(), saveat=0.1) -plot(sol, vars=(0, [1,2]), lw=2) +plot(sol, idxs=[1,2], lw=2) ``` ![Brusselator issue 1](../assets/brusselator_issue_1.svg) @@ -40,7 +41,7 @@ u₀map = deterministic_IC(u₀, closed_raw_eqs) oprob = ODEProblem(closed_raw_eqs, u₀map, tspan, p) sol = solve(oprob, Tsit5(), saveat=0.1) -plot(sol, vars=(0, [1,2]), lw=2, legend=:bottomright) +plot(sol, idxs=[1,2], lw=2, legend=:bottomright) ``` ![Brusselator issue 2](../assets/brusselator_issue_2.svg) @@ -48,66 +49,66 @@ We observe sustained oscillatory behaviour instead of the expected damped oscill Normal closure is also quite fragile. This can be seen by simply including the combinatorial scaling of the mass-action propensity functions with `combinatoric_ratelaw=true` which leads to unphysical sustained oscillatory trajectories: ```julia -raw_eqs = generate_raw_moment_eqs(rn, 2, combinatoric_ratelaw=true) +raw_eqs = generate_raw_moment_eqs(rn, 2, combinatoric_ratelaws=true) closed_raw_eqs = moment_closure(raw_eqs, "normal") u₀map = deterministic_IC(u₀, closed_raw_eqs) oprob = ODEProblem(closed_raw_eqs, u₀map, tspan, p) sol = solve(oprob, Tsit5(), saveat=0.1) -plot(sol, vars=(0, [1,2]), lw=2) +plot(sol, idxs=[1,2], lw=2) ``` ![Brusselator issue 3](../assets/brusselator_issue_3.svg) Nevertheless, this can be improved upon by increasing the order of moment expansion: ```julia -raw_eqs = generate_raw_moment_eqs(rn, 3, combinatoric_ratelaw=true) +raw_eqs = generate_raw_moment_eqs(rn, 3, combinatoric_ratelaws=true) closed_raw_eqs = moment_closure(raw_eqs, "normal") u₀map = deterministic_IC(u₀, closed_raw_eqs) oprob = ODEProblem(closed_raw_eqs, u₀map, tspan, p) sol = solve(oprob, Tsit5(), saveat=0.1) -plot(sol, vars=(0, [1,2]), lw=2, legend=:bottomright) +plot(sol, idxs=[1,2], lw=2, legend=:bottomright) ``` ![Brusselator issue 4](../assets/brusselator_issue_4.svg) Some dampening in the system is now visible. Increasing the expansion order to `4` finally leads to physically sensible results: ```julia -raw_eqs = generate_raw_moment_eqs(rn, 4, combinatoric_ratelaw=true) +raw_eqs = generate_raw_moment_eqs(rn, 4, combinatoric_ratelaws=true) closed_raw_eqs = moment_closure(raw_eqs, "normal") u₀map = deterministic_IC(u₀, closed_raw_eqs) oprob = ODEProblem(closed_raw_eqs, u₀map, tspan, p) sol = solve(oprob, Tsit5(), saveat=0.1) -plot(sol, vars=(0, [1,2]), lw=2) +plot(sol, idxs=[1,2], lw=2) ``` ![Brusselator issue 5](../assets/brusselator_issue_5.svg) For dessert, we consider unphysical divergent trajectories—a frequent problem with moment equations [3]. A good example is the second-order moment expansion including the combinatorial scaling of propensities with log-normal closure applied: ```julia -raw_eqs = generate_raw_moment_eqs(rn, 2, combinatoric_ratelaw=true) +raw_eqs = generate_raw_moment_eqs(rn, 2, combinatoric_ratelaws=true) closed_raw_eqs = moment_closure(raw_eqs, "log-normal") u₀map = deterministic_IC(u₀, closed_raw_eqs) oprob = ODEProblem(closed_raw_eqs, u₀map, tspan, p) sol = solve(oprob, Rodas4P(), saveat=0.1) -plot(sol, vars=(0, [1,2]), lw=2) +plot(sol, idxs=[1,2], lw=2) ``` ![Brusselator issue 6](../assets/brusselator_issue_6.svg) In contrast to normal closure, increasing the expansion order makes the problem worse: ```julia -raw_eqs = generate_raw_moment_eqs(rn, 3, combinatoric_ratelaw=true) +raw_eqs = generate_raw_moment_eqs(rn, 3, combinatoric_ratelaws=true) closed_raw_eqs = moment_closure(raw_eqs, "log-normal") u₀map = deterministic_IC(u₀, closed_raw_eqs) oprob = ODEProblem(closed_raw_eqs, u₀map, tspan, p) sol = solve(oprob, Rodas4P(), saveat=0.1) -plot(sol, vars=(0, [1,2]), lw=2) +plot(sol, idxs=[1,2], lw=2) ``` ```julia ┌ Warning: Interrupted. Larger maxiters is needed. diff --git a/docs/src/tutorials/derivative_matching_example.md b/docs/src/tutorials/derivative_matching_example.md index cb5704c..57eace8 100644 --- a/docs/src/tutorials/derivative_matching_example.md +++ b/docs/src/tutorials/derivative_matching_example.md @@ -12,9 +12,10 @@ The reaction network and its parameters can be defined as follows: using Catalyst rn = @reaction_network begin + @parameters c₁ c₂ (c₁), x₁ → 2x₁+x₂ (c₂), x₁+x₂ → x₂ -end c₁ c₂ +end # parameter values p = [1.0, 1.0] @@ -133,7 +134,7 @@ true ``` The last ingredient we need for a proper comparison between the second and third order moment expansions is a reference value predicted by the SSA. We can simulate $10^5$ SSA trajectories as follows: ```julia -using DiffEqJump +using JumpProcesses dprob = DiscreteProblem(rn, u0, tspan, p) jprob = JumpProblem(rn, dprob, Direct(), save_positions=(false, false)) diff --git a/docs/src/tutorials/geometric_reactions+conditional_closures.md b/docs/src/tutorials/geometric_reactions+conditional_closures.md index 31f2134..3e26e7a 100644 --- a/docs/src/tutorials/geometric_reactions+conditional_closures.md +++ b/docs/src/tutorials/geometric_reactions+conditional_closures.md @@ -28,11 +28,12 @@ using MomentClosure, Catalyst, Distributions, JumpProcesses, DiffEqBase, Ordinar m = rand(Distributions.Geometric(1/(1+b))) rn = @reaction_network begin + @parameters k_on k_off k_p γ_p k_on*(1-g), 0 --> g # G* -> G k_off*P^2, g --> 0 # G -> G* k_p, g --> g + $m*P # G -> G + mP, m ~ Geometric(p) γ_p, P --> 0 # P -> ∅ -end k_on k_off k_p γ_p +end ``` We can now generate the raw moment equations up to third order: ```julia @@ -40,62 +41,18 @@ eqs = generate_raw_moment_eqs(rn, 3) latexify(eqs) ``` ```math - \begin{align*} -\frac{d\mu_{1 0}}{dt} =& k_{on} - k_{on} \mu_{1 0} - k_{off} \mu_{1 2} \\ +\frac{d\mu_{1 0}}{dt} =& k_{on} - k_{off} \mu_{1 2} - k_{on} \mu_{1 0} \\ \frac{d\mu_{0 1}}{dt} =& b k_{p} \mu_{1 0} - \gamma_{p} \mu_{0 1} \\ -\frac{d\mu_{2 0}}{dt} =& k_{on} + k_{on} \mu_{1 0} + k_{off} \mu_{1 2} - 2 k_{on} \mu_{2 0} - 2 k_{off} \mu_{2 2} \\ -\frac{d\mu_{1 1}}{dt} =& k_{on} \mu_{0 1} + b k_{p} \mu_{2 0} - k_{off} \mu_{1 3} - k_{on} \mu_{1 1} - \gamma_{p} \mu_{1 1} \\ -\frac{d\mu_{0 2}}{dt} =& \gamma_{p} \mu_{0 1} + \frac{k_{p} \left( \frac{-1}{1 + b} \right)^{2} \mu_{1 0} + 2 b k_{p} \mu_{1 0} + b k_{p} \left( \frac{-1}{1 + b} \right)^{2} \mu_{1 0} - k_{p} \mu_{1 0}}{b \left( \frac{1}{1 + b} \right)^{2} + \left( \frac{1}{1 + b} \right)^{2}} + 2 b k_{p} \mu_{1 1} - 2 \gamma_{p} \mu_{0 2} \\ -\frac{d\mu_{3 0}}{dt} =& k_{on} + 2 k_{on} \mu_{1 0} + 3 k_{off} \mu_{2 2} - k_{off} \mu_{1 2} - 3 k_{on} \mu_{3 0} - 3 k_{off} \mu_{3 2} \\ -\frac{d\mu_{2 1}}{dt} =& k_{off} \mu_{1 3} + k_{on} \mu_{0 1} + k_{on} \mu_{1 1} + b k_{p} \mu_{3 0} - 2 k_{off} \mu_{2 3} - \gamma_{p} \mu_{2 1} - 2 k_{on} \mu_{2 1} \\ -\frac{d\mu_{1 2}}{dt} =& k_{on} \mu_{0 2} + \gamma_{p} \mu_{1 1} + \frac{k_{p} \left( \frac{-1}{1 + b} \right)^{2} \mu_{2 0} + 2 b k_{p} \mu_{2 0} + b k_{p} \left( \frac{-1}{1 + b} \right)^{2} \mu_{2 0} - k_{p} \mu_{2 0}}{b \left( \frac{1}{1 + b} \right)^{2} + \left( \frac{1}{1 + b} \right)^{2}} + 2 b k_{p} \mu_{2 1} - k_{off} \mu_{1 4} - k_{on} \mu_{1 2} - 2 \gamma_{p} \mu_{1 2} \\ -\frac{d\mu_{0 3}}{dt} =& \frac{3 k_{p} \left( \frac{-1}{1 + b} \right)^{2} \mu_{1 1} + 6 b k_{p} \mu_{1 1} + 3 b k_{p} \left( \frac{-1}{1 + b} \right)^{2} \mu_{1 1} - 3 k_{p} \mu_{1 1}}{b \left( \frac{1}{1 + b} \right)^{2} + \left( \frac{1}{1 + b} \right)^{2}} + \frac{k_{p} \left( \frac{-1}{1 + b} \right)^{3} \mu_{1 0} + 6 b k_{p} \mu_{1 0} + 7 k_{p} \left( \frac{-1}{1 + b} \right)^{2} \mu_{1 0} + b k_{p} \left( \frac{-1}{1 + b} \right)^{3} \mu_{1 0} + 7 b k_{p} \left( \frac{-1}{1 + b} \right)^{2} \mu_{1 0} - 6 k_{p} \mu_{1 0}}{b \left( \frac{1}{1 + b} \right)^{3} + \left( \frac{1}{1 + b} \right)^{3}} + 3 \gamma_{p} \mu_{0 2} + 3 b k_{p} \mu_{1 2} - \gamma_{p} \mu_{0 1} - 3 \gamma_{p} \mu_{0 3} -\end{align*} -``` -Note that certain symbolic expressions have not been fully simplified: this happens as [Symbolics.jl](https://github.com/JuliaSymbolics/Symbolics.jl) used internally is not yet able to simplify more complicated fractions to the extent that is required here. For reference, the simplified moment equations take the form: -```math -\begin{align*} -\frac{d\mu{_{10}}}{dt} =& k_{on} - k_{off} \mu{_{12}} - k_{on} \mu{_{10}} \\ -\frac{d\mu{_{01}}}{dt} =& b k_{p} \mu{_{10}} - \gamma_{p} \mu{_{01}} \\ -\frac{d\mu{_{20}}}{dt} =& k_{on} + k_{off} \mu{_{12}} + k_{on} \mu{_{10}} - 2 k_{off} \mu{_{22}} - 2 k_{on} \mu{_{20}} \\ -\frac{d\mu{_{11}}}{dt} =& k_{on} \mu{_{01}} + b k_{p} \mu{_{20}} - k_{off} \mu{_{13}} - k_{on} \mu{_{11}} - \gamma_{p} \mu{_{11}} \\ -\frac{d\mu{_{02}}}{dt} =& \gamma_{p} \mu{_{01}} + b k_{p} \mu{_{10}} + 2 b k_{p} \mu{_{11}} + 2 k_{p} \mu{_{10}} b^{2} - 2 \gamma_{p} \mu{_{02}} \\ -\frac{d\mu{_{30}}}{dt} =& k_{on} + 3 k_{off} \mu{_{22}} + 2 k_{on} \mu{_{10}} - k_{off} \mu{_{12}} - 3 k_{off} \mu{_{32}} - 3 k_{on} \mu{_{30}} \\ -\frac{d\mu{_{21}}}{dt} =& k_{off} \mu{_{13}} + k_{on} \mu{_{01}} + k_{on} \mu{_{11}} + b k_{p} \mu{_{30}} - 2 k_{off} \mu{_{23}} - 2 k_{on} \mu{_{21}} - \gamma_{p} \mu{_{21}} \\ -\frac{d\mu{_{12}}}{dt} =& k_{on} \mu{_{02}} + \gamma_{p} \mu{_{11}} + b k_{p} \mu{_{20}} + 2 b k_{p} \mu{_{21}} + 2 k_{p} \mu{_{20}} b^{2} - k_{off} \mu{_{14}} - k_{on} \mu{_{12}} - 2 \gamma_{p} \mu{_{12}} \\ -\frac{d\mu{_{03}}}{dt} =& b k_{p} \mu{_{10}} + 3 \gamma_{p} \mu{_{02}} + 3 b k_{p} \mu{_{11}} + 3 b k_{p} \mu{_{12}} + 6 k_{p} \mu{_{10}} b^{2} + 6 k_{p} \mu{_{10}} b^{3} + 6 k_{p} \mu{_{11}} b^{2} - \gamma_{p} \mu{_{01}} - 3 \gamma_{p} \mu{_{03}} -\end{align*} -``` -The fractions appeared in our equations as we defined the success probability $p$ of the geometric distribution as $p = \frac{1}{1+b}$. This could have been avoided by leaving $p$ in the symbolic expressions and substituting its numerical value expressed in terms of the mean $b$ later on. We proceed to do that and redefine the reaction system accordingly, making the moment equations easier to handle: -```julia -@parameters p -m = rand(Distributions.Geometric(p)) - -rn = @reaction_network begin - k_on*(1-g), 0 --> g - k_off*P^2, g --> 0* - k_p, g --> g + $m*P - γ_p, P --> 0 -end k_on k_off k_p γ_p - -eqs = generate_raw_moment_eqs(rn, 3) -latexify(eqs) -``` -```math -\begin{align*} -\frac{d\mu_{1 0}}{dt} =& k_{on} - k_{on} \mu_{1 0} - k_{off} \mu_{1 2} \\ -\frac{d\mu_{0 1}}{dt} =& k_{p} p^{-1} \mu_{1 0} - k_{p} \mu_{1 0} - \gamma_{p} \mu_{0 1} \\ -\frac{d\mu_{2 0}}{dt} =& k_{on} + k_{on} \mu_{1 0} + k_{off} \mu_{1 2} - 2 k_{on} \mu_{2 0} - 2 k_{off} \mu_{2 2} \\ -\frac{d\mu_{1 1}}{dt} =& k_{on} \mu_{0 1} + k_{p} p^{-1} \mu_{2 0} - k_{off} \mu_{1 3} - k_{on} \mu_{1 1} - \gamma_{p} \mu_{1 1} - k_{p} \mu_{2 0} \\ -\frac{d\mu_{0 2}}{dt} =& k_{p} \mu_{1 0} + \gamma_{p} \mu_{0 1} + 2 k_{p} p^{-2} \mu_{1 0} + 2 k_{p} p^{-1} \mu_{1 1} - 2 k_{p} \mu_{1 1} - 2 \gamma_{p} \mu_{0 2} - 3 k_{p} p^{-1} \mu_{1 0} \\ -\frac{d\mu_{3 0}}{dt} =& k_{on} + 2 k_{on} \mu_{1 0} + 3 k_{off} \mu_{2 2} - k_{off} \mu_{1 2} - 3 k_{on} \mu_{3 0} - 3 k_{off} \mu_{3 2} \\ -\frac{d\mu_{2 1}}{dt} =& k_{off} \mu_{1 3} + k_{on} \mu_{0 1} + k_{on} \mu_{1 1} + k_{p} p^{-1} \mu_{3 0} - 2 k_{off} \mu_{2 3} - \gamma_{p} \mu_{2 1} - k_{p} \mu_{3 0} - 2 k_{on} \mu_{2 1} \\ -\frac{d\mu_{1 2}}{dt} =& k_{on} \mu_{0 2} + \gamma_{p} \mu_{1 1} + k_{p} \mu_{2 0} + 2 k_{p} p^{-1} \mu_{2 1} + 2 k_{p} p^{-2} \mu_{2 0} - k_{off} \mu_{1 4} - k_{on} \mu_{1 2} - 2 \gamma_{p} \mu_{1 2} - 2 k_{p} \mu_{2 1} - 3 k_{p} p^{-1} \mu_{2 0} \\ -\frac{d\mu_{0 3}}{dt} =& 3 k_{p} \mu_{1 1} + 3 \gamma_{p} \mu_{0 2} + 7 k_{p} p^{-1} \mu_{1 0} + 3 k_{p} p^{-1} \mu_{1 2} + 6 k_{p} p^{-3} \mu_{1 0} + 6 k_{p} p^{-2} \mu_{1 1} - \gamma_{p} \mu_{0 1} - k_{p} \mu_{1 0} - 3 k_{p} \mu_{1 2} - 3 \gamma_{p} \mu_{0 3} - 9 k_{p} p^{-1} \mu_{1 1} - 12 k_{p} p^{-2} \mu_{1 0} +\frac{d\mu_{2 0}}{dt} =& k_{on} + k_{off} \mu_{1 2} + k_{on} \mu_{1 0} - 2 k_{on} \mu_{2 0} - 2 k_{off} \mu_{2 2} \\ +\frac{d\mu_{1 1}}{dt} =& k_{on} \mu_{0 1} + b k_{p} \mu_{2 0} - k_{on} \mu_{1 1} - \gamma_{p} \mu_{1 1} - k_{off} \mu_{1 3} \\ +\frac{d\mu_{0 2}}{dt} =& \gamma_{p} \mu_{0 1} + b k_{p} \mu_{1 0} + 2 b k_{p} \mu_{1 1} + 2 k_{p} b^{2} \mu_{1 0} - 2 \gamma_{p} \mu_{0 2} \\ +\frac{d\mu_{3 0}}{dt} =& k_{on} + 3 k_{off} \mu_{2 2} + 2 k_{on} \mu_{1 0} - k_{off} \mu_{1 2} - 3 k_{off} \mu_{3 2} - 3 k_{on} \mu_{3 0} \\ +\frac{d\mu_{2 1}}{dt} =& k_{on} \mu_{0 1} + k_{on} \mu_{1 1} + k_{off} \mu_{1 3} + b k_{p} \mu_{3 0} - 2 k_{on} \mu_{2 1} - \gamma_{p} \mu_{2 1} - 2 k_{off} \mu_{2 3} \\ +\frac{d\mu_{1 2}}{dt} =& k_{on} \mu_{0 2} + \gamma_{p} \mu_{1 1} + b k_{p} \mu_{2 0} + 2 k_{p} b^{2} \mu_{2 0} + 2 b k_{p} \mu_{2 1} - k_{off} \mu_{1 4} - k_{on} \mu_{1 2} - 2 \gamma_{p} \mu_{1 2} \\ +\frac{d\mu_{0 3}}{dt} =& b k_{p} \mu_{1 0} + 3 \gamma_{p} \mu_{0 2} + 3 b k_{p} \mu_{1 1} + 3 b k_{p} \mu_{1 2} + 6 k_{p} b^{2} \mu_{1 0} + 6 k_{p} b^{3} \mu_{1 0} + 6 k_{p} b^{2} \mu_{1 1} - \gamma_{p} \mu_{0 1} - 3 \gamma_{p} \mu_{0 3} \end{align*} ``` - A lot of information in this system of ODEs is redundant as the gene state is a Bernoulli variable that (in our case) has the following properties: ```math \begin{align*} @@ -113,12 +70,12 @@ latexify(clean_eqs) ``` ```math \begin{align*} -\frac{d\mu_{1 0}}{dt} =& k_{on} - k_{on} \mu_{1 0} - k_{off} \mu_{1 2} \\ -\frac{d\mu_{0 1}}{dt} =& k_{p} p^{-1} \mu_{1 0} - \gamma_{p} \mu_{0 1} - k_{p} \mu_{1 0} \\ -\frac{d\mu_{1 1}}{dt} =& k_{on} \mu_{0 1} + k_{p} p^{-1} \mu_{1 0} - k_{off} \mu_{1 3} - k_{p} \mu_{1 0} - k_{on} \mu_{1 1} - \gamma_{p} \mu_{1 1} \\ -\frac{d\mu_{0 2}}{dt} =& \gamma_{p} \mu_{0 1} + k_{p} \mu_{1 0} + 2 k_{p} p^{-2} \mu_{1 0} + 2 k_{p} p^{-1} \mu_{1 1} - 2 k_{p} \mu_{1 1} - 2 \gamma_{p} \mu_{0 2} - 3 k_{p} p^{-1} \mu_{1 0} \\ -\frac{d\mu_{1 2}}{dt} =& k_{on} \mu_{0 2} + k_{p} \mu_{1 0} + \gamma_{p} \mu_{1 1} + 2 k_{p} p^{-2} \mu_{1 0} + 2 k_{p} p^{-1} \mu_{1 1} - 2 k_{p} \mu_{1 1} - k_{off} \mu_{1 4} - k_{on} \mu_{1 2} - 2 \gamma_{p} \mu_{1 2} - 3 k_{p} p^{-1} \mu_{1 0} \\ -\frac{d\mu_{0 3}}{dt} =& 3 k_{p} \mu_{1 1} + 3 \gamma_{p} \mu_{0 2} + 7 k_{p} p^{-1} \mu_{1 0} + 3 k_{p} p^{-1} \mu_{1 2} + 6 k_{p} p^{-3} \mu_{1 0} + 6 k_{p} p^{-2} \mu_{1 1} - \gamma_{p} \mu_{0 1} - k_{p} \mu_{1 0} - 3 k_{p} \mu_{1 2} - 3 \gamma_{p} \mu_{0 3} - 9 k_{p} p^{-1} \mu_{1 1} - 12 k_{p} p^{-2} \mu_{1 0} +\frac{d\mu_{1 0}}{dt} =& k_{on} - k_{off} \mu_{1 2} - k_{on} \mu_{1 0} \\ +\frac{d\mu_{0 1}}{dt} =& b k_{p} \mu_{1 0} - \gamma_{p} \mu_{0 1} \\ +\frac{d\mu_{1 1}}{dt} =& k_{on} \mu_{0 1} + b k_{p} \mu_{1 0} - k_{on} \mu_{1 1} - \gamma_{p} \mu_{1 1} - k_{off} \mu_{1 3} \\ +\frac{d\mu_{0 2}}{dt} =& \gamma_{p} \mu_{0 1} + b k_{p} \mu_{1 0} + 2 b k_{p} \mu_{1 1} + 2 k_{p} b^{2} \mu_{1 0} - 2 \gamma_{p} \mu_{0 2} \\ +\frac{d\mu_{1 2}}{dt} =& k_{on} \mu_{0 2} + \gamma_{p} \mu_{1 1} + b k_{p} \mu_{1 0} + 2 b k_{p} \mu_{1 1} + 2 k_{p} b^{2} \mu_{1 0} - k_{off} \mu_{1 4} - k_{on} \mu_{1 2} - 2 \gamma_{p} \mu_{1 2} \\ +\frac{d\mu_{0 3}}{dt} =& b k_{p} \mu_{1 0} + 3 \gamma_{p} \mu_{0 2} + 3 b k_{p} \mu_{1 1} + 3 b k_{p} \mu_{1 2} + 6 k_{p} b^{2} \mu_{1 0} + 6 k_{p} b^{3} \mu_{1 0} + 6 k_{p} b^{2} \mu_{1 1} - \gamma_{p} \mu_{0 1} - 3 \gamma_{p} \mu_{0 3} \end{align*} ``` The system of ODEs is now much simpler and we can see that there are two higher-order moments we need to truncate: $\mu_{13}$ and $μ_{14}$. We consider normal, derivative matching, conditional gaussian and conditional derivative matching closures to see how well they compare. First we apply different closures and print out the corresponding higher-order moment expressions in order to check that our results are consistent with those published in [1]. @@ -143,8 +100,8 @@ latexify(dm_eqs, :closure) ``` ```math \begin{align*} -\mu{_{13}} =& \mu{_{03}} \mu{_{10}} \mu{_{01}}^{3} \mu{_{02}}^{-3} \mu{_{11}}^{-3} \mu{_{12}}^{3} \\ -\mu{_{14}} =& \mu{_{04}} \mu{_{01}}^{-4} \mu{_{02}}^{6} \mu{_{03}}^{-4} \mu{_{10}}^{-1} \mu{_{11}}^{4} \mu{_{12}}^{-6} \mu{_{13}}^{4} +\mu_{1 3} =& \frac{\mu_{0 1}^{3} \mu_{1 2}^{3} \mu_{0 3} \mu_{1 0}}{\mu_{0 2}^{3} \mu_{1 1}^{3}} \\ +\mu_{1 4} =& \frac{\mu_{0 2}^{6} \mu_{1 1}^{4} \mu_{1 3}^{4} \mu_{0 4}}{\mu_{0 1}^{4} \mu_{0 3}^{4} \mu_{1 2}^{6} \mu_{1 0}} \end{align*} ``` For conditional gaussian closure: @@ -154,8 +111,8 @@ latexify(cond_gaussian_eqs, :closure) ``` ```math \begin{align*} -\mu{_{13}} =& 3 \mu{_{11}} \mu{_{12}} \mu{_{10}}^{-1} - 2 \mu{_{10}}^{-2} \mu{_{11}}^{3} \\ -\mu{_{14}} =& 3 \mu{_{10}}^{-1} \mu{_{12}}^{2} + 6 \mu{_{10}}^{-3} \mu{_{11}}^{4} + 4 \mu{_{11}} \mu{_{13}} \mu{_{10}}^{-1} - 12 \mu{_{12}} \mu{_{10}}^{-2} \mu{_{11}}^{2} +\mu_{1 3} =& \frac{-2 \mu_{1 1}^{3}}{\mu_{1 0}^{2}} + \frac{3 \mu_{1 1} \mu_{1 2}}{\mu_{1 0}} \\ +\mu_{1 4} =& \frac{3 \mu_{1 2}^{2}}{\mu_{1 0}} + \frac{6 \mu_{1 1}^{4}}{\mu_{1 0}^{3}} + \frac{-12 \mu_{1 1}^{2} \mu_{1 2}}{\mu_{1 0}^{2}} + \frac{4 \mu_{1 1} \mu_{1 3}}{\mu_{1 0}} \end{align*} ``` And, finally, for conditional derivative matching: @@ -165,15 +122,14 @@ latexify(cond_dm_eqs, :closure) ``` ```math \begin{align*} -\mu{_{13}} =& \mu{_{10}} \mu{_{11}}^{-3} \mu{_{12}}^{3} \\ -\mu{_{14}} =& \mu{_{10}}^{-1} \mu{_{11}}^{4} \mu{_{12}}^{-6} \mu{_{13}}^{4} +\mu_{1 3} =& \frac{\mu_{0 1}^{3} \mu_{1 2}^{3} \mu_{0 3} \mu_{1 0}}{\mu_{0 2}^{3} \mu_{1 1}^{3}} \\ +\mu_{1 4} =& \frac{\mu_{0 2}^{6} \mu_{1 1}^{4} \mu_{1 3}^{4} \mu_{0 4}}{\mu_{0 1}^{4} \mu_{0 3}^{4} \mu_{1 2}^{6} \mu_{1 0}} \end{align*} ``` All these results are consistent with [1], reassuring that the model and closures are implemented as intended. Finally, we can proceed to solve the resulting ODEs and compare the resulting means and standard deviations to the SSA. Following Soltani et al. [1], we define all model parameters as: ```julia mean_p = 200 mean_b = 70 -p_val = 1/(1+mean_b) # note we use p = 1/(1+b) γ_p_val = 1 k_off_val = 0.001 k_on_val = 0.05 @@ -184,7 +140,7 @@ symmap = [:k_on => k_on_val, :k_off => k_off_val, :k_p => k_p_val, :γ_p => γ_p_val, - :p => p_val] + :b => mean_b] pmap = symmap_to_varmap(rn, symmap) # initial gene state and protein number, order [g, p] @@ -234,7 +190,7 @@ for closure in ["normal", "derivative matching", # μ₀₁ is 2nd and μ₀₂ is 4th element in sol # can check the order with `closed_eqs.odes.states` - plt_m = plot!(plt_m, sol, vars=(0, 2), label=closure) + plt_m = plot!(plt_m, sol, idxs=[2], label=closure) plt_std = plot!(plt_std, sol.t, sqrt.(sol[4, :] .- sol[2, :].^2), label=closure) end @@ -267,6 +223,7 @@ l = rand(Distributions.Geometric(1/(1+b_y))) # g_y - gene state of Y protein producing gene # x, y - proteins X and Y rn = @reaction_network begin + @parameters kx_on kx_off ky_on ky_off k_x γ_x k_y γ_y kx_on*(1-g_x)*y, 0 --> g_x # 0 -> g_x kx_off, g_x --> 0 # g_x -> 0 ky_on*(1-g_y), 0 --> g_y # 0 -> g_y @@ -275,7 +232,7 @@ rn = @reaction_network begin γ_x, x --> 0 # x -> 0 k_y*g_y, 0 --> $l*y # 0 -> ly, l ~ Geometric(mean_b_y) γ_y, y --> 0 # y -> 0 -end kx_on kx_off ky_on ky_off k_x γ_x k_y γ_y +end # both g_x and g_y are Bernoulli random variables binary_vars = [1, 2] @@ -351,7 +308,7 @@ for closure in ["derivative matching", "conditional derivative matching"] sol = solve(oprob, Tsit5(), saveat=0.1) # μ₀₀₀₁ is the 4th and μ₀₀₀₂ is the 12th element in sol (can check with closed_eqs.odes.states) - plt_m = plot!(plt_m, sol, vars=(0, 4), label=closure) + plt_m = plot!(plt_m, sol, idxs=[4], label=closure) plt_std = plot!(plt_std, sol.t, sqrt.(sol[12, :] .- sol[4, :].^2), label=closure) end diff --git a/docs/src/tutorials/parameter_estimation_SDE.md b/docs/src/tutorials/parameter_estimation_SDE.md index 1c00b3e..1988cbe 100644 --- a/docs/src/tutorials/parameter_estimation_SDE.md +++ b/docs/src/tutorials/parameter_estimation_SDE.md @@ -111,7 +111,7 @@ closed_moment_prob = ODEProblem(LV_moments, u0map, (0.0, Tf), zeros(5)) function obj_MCA(p) prob = remake(closed_moment_prob; p=p) sol = solve(prob, Tsit5(), saveat = t_data) - if sol.retcode == :Success + if sol.retcode == ReturnCode.Success obj = sum(norm(sol.u[i][1:2] - means[i])^2 for i in 1:length(t_data)) obj += 1e4*sum((sol.u[i][3] - sol.u[i][1]^2 - vars[i][1])^2 for i in 1:length(t_data)) obj += 1e4*sum((sol.u[i][5] - sol.u[i][2]^2 - vars[i][2])^2 for i in 1:length(t_data)) diff --git a/docs/src/tutorials/time-dependent_propensities.md b/docs/src/tutorials/time-dependent_propensities.md index 58f8ee0..5e360e5 100644 --- a/docs/src/tutorials/time-dependent_propensities.md +++ b/docs/src/tutorials/time-dependent_propensities.md @@ -20,12 +20,13 @@ where $c_2^0$ is a fixed value and the sinusoidal modulation with frequency $\om using Catalyst, MomentClosure, Latexify rn = @reaction_network begin + @parameters c₁ c₂ c₃ c₄ Ω ω τ (c₁/Ω^2), 2X + Y → 3X (c₂*(1+0.5*sin(ω*(t<τ)*t))), X → Y (c₃*Ω, c₄), 0 ↔ X -end c₁ c₂ c₃ c₄ Ω ω τ +end -raw_eqs = generate_raw_moment_eqs(rn, 2, combinatoric_ratelaw=false) +raw_eqs = generate_raw_moment_eqs(rn, 2, combinatoric_ratelaws=false) latexify(raw_eqs) ``` ```math @@ -60,13 +61,13 @@ oprob = ODEProblem(closed_raw_eqs, u₀map, tspan, p) # solve using Tsit5 solver sol = solve(oprob, Tsit5(), saveat=0.2) -plot(sol, vars=(0, [1,2]), lw=2) +plot(sol, idxs=[1,2], lw=2) ``` ![Time-dependent Brusselator normal](../assets/Brusselator_time-dependent_normal.svg) It would be great to compare our results to the true dynamics. Using DifferentialEquations, we can run a modified SSA taking into account the time-dependent propensity functions ([VariableRateJumps](https://diffeq.sciml.ai/stable/tutorials/discrete_stochastic_example/#Adding-a-VariableRateJump)). This requires some care and is done by combining an `ODEProblem` with a `JumpProblem`: ```julia -using DiffEqJump +using JumpProcesses # convert ReactionSystem into JumpSystem jsys = convert(JumpSystem, rn, combinatoric_ratelaws=false) diff --git a/docs/src/tutorials/using_momentclosure.md b/docs/src/tutorials/using_momentclosure.md index 9341534..390e8bf 100644 --- a/docs/src/tutorials/using_momentclosure.md +++ b/docs/src/tutorials/using_momentclosure.md @@ -19,10 +19,11 @@ The terminology and the notation used throughout is consistent with the **Theory using Catalyst rn = @reaction_network begin # including system-size parameter Ω + @parameters c₁ c₂ c₃ c₄ Ω (c₁/Ω^2), 2X + Y → 3X (c₂), X → Y (c₃*Ω, c₄), 0 ↔ X -end c₁ c₂ c₃ c₄ Ω +end ``` The returned `rn` is an instance of [`ModelingToolkit.ReactionSystem`](https://catalyst.sciml.ai/stable/api/catalyst_api/#ModelingToolkit.ReactionSystem). The net stoichiometry matrix and an array of the corresponding propensities, if needed, can be extracted directly from the model using [`Catalyst.netstoichmat`](https://catalyst.sciml.ai/dev/api/catalyst_api/#Catalyst.netstoichmat) and MomentClosure function [`propensities`](@ref) respectively. @@ -37,7 +38,7 @@ We can now obtain the moment equations. The system follows the law of mass actio Let's start with the raw moment equations which we choose to generate up to second order ($m=2$): ```julia using MomentClosure -raw_eqs = generate_raw_moment_eqs(rn, 2, combinatoric_ratelaw=false) +raw_eqs = generate_raw_moment_eqs(rn, 2, combinatoric_ratelaws=false) ``` Note that we have set `combinatoric_ratelaw=false` in order to ignore the factorial scaling factors which [Catalyst adds to mass-action reactions](https://catalyst.sciml.ai/stable/tutorials/models/#Reaction-rate-laws-used-in-simulations). The function [`generate_raw_moment_eqs`](@ref) returns an instance of [`RawMomentEquations`](@ref MomentClosure.RawMomentEquations) that contains a [`ModelingToolkit.ODESystem`](https://mtk.sciml.ai/stable/systems/ODESystem/) composed of all the moment equations (accessed by `raw_eqs.odes`). @@ -75,7 +76,7 @@ Coming back to the generated moment equations, we observe that they depend on hi The corresponding central moment equations can also be easily generated: ```julia -central_eqs = generate_central_moment_eqs(rn, 2, combinatoric_ratelaw=false) +central_eqs = generate_central_moment_eqs(rn, 2, combinatoric_ratelaws=false) ``` Note that in case of non-polynomial propensity functions the [Taylor expansion order $q$](@ref central_moment_eqs) must also be specified, see the [P53 system example](P53_system_example.md) for more details. Luckily, the Brusselator contains only mass-action reactions and hence $q$ is automatically determined by the highest order (polynomial) propensity. The function [`generate_central_moment_eqs`](@ref) returns an instance of [`CentralMomentEquations`](@ref MomentClosure.CentralMomentEquations). As before, we can visualise the central moment equations: ```julia @@ -189,7 +190,7 @@ using OrdinaryDiffEq sol = solve(oprob, Tsit5(), saveat=0.1) using Plots -plot(sol, vars=(0, [1,2]), lw=2) +plot(sol, idxs=[1,2], lw=2) ``` ![Brusselator means 1](../assets/brusselator_means_1.svg) @@ -213,7 +214,7 @@ dprob = DiscreteProblem(jsys, u₀, tspan, p) # same parameters as defined earli jprob = JumpProblem(jsys, dprob, Direct(), save_positions=(false, false)) # define an EnsembleProblem to simulate multiple trajectories -ensembleprob = EnsembleProblem(jprob) +ensembleprob = EnsembleProblem(jprob) # simulate 10000 SSA trajectories @time sol_SSA = solve(ensembleprob, SSAStepper(), saveat=0.1, trajectories=10000) From c3ca652419311a44b7007a8c5f5bb80fe461d552 Mon Sep 17 00:00:00 2001 From: Augustinas Sukys Date: Sat, 6 May 2023 01:35:10 +0100 Subject: [PATCH 16/16] update Project.toml --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 37c9a79..0759ce6 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "MomentClosure" uuid = "01a1b25a-ecf0-48c5-ae58-55bfd5393600" authors = ["Augustinas Sukys"] -version = "0.2.0" +version = "0.3.0" [deps] Catalyst = "479239e8-5488-4da2-87a7-35f2df7eef83"