From c15f84e265ce4e881b2cc63ece2445d05bfc2ce8 Mon Sep 17 00:00:00 2001 From: Torkel Date: Sun, 17 Sep 2023 15:28:48 -0400 Subject: [PATCH 01/41] start noise scaling remake --- src/reactionsystem.jl | 20 ++++++++++++++++++-- 1 file changed, 18 insertions(+), 2 deletions(-) diff --git a/src/reactionsystem.jl b/src/reactionsystem.jl index 770eb938ba..7f8691b391 100644 --- a/src/reactionsystem.jl +++ b/src/reactionsystem.jl @@ -54,6 +54,17 @@ end # variables drop_dynamics(s) = isconstant(s) || isbc(s) || (!isspecies(s)) +### Other Metadata Fields ### + +# Denotes that a parameter controls the scaling of noise in the CLE. +struct NoiseScalingParameter end +Symbolics.option_to_metadata_type(::Val{:isnoisescalingparameter}) = NoiseScalingParameter + +isnoisescalingparameter(s::Num) = isnoisescalingparameter(MT.value(s)) +function isnoisescalingparameter(s) + MT.getmetadata(s, NoiseScalingParameter, false) +end + """ $(TYPEDEF) @@ -1489,7 +1500,7 @@ Notes: differential equations. """ function Base.convert(::Type{<:SDESystem}, rs::ReactionSystem; - noise_scaling = nothing, name = nameof(rs), + noise_scaling = get_noise_scaling(rs), name = nameof(rs), combinatoric_ratelaws = get_combinatoric_ratelaws(rs), include_zero_odes = true, checks = false, remove_conserved = false, default_u0 = Dict(), default_p = Dict(), defaults = _merge(Dict(default_u0), Dict(default_p)), @@ -1533,6 +1544,11 @@ function Base.convert(::Type{<:SDESystem}, rs::ReactionSystem; kwargs...) end +# Extracts any noise scaling parameters from a reaction system. +function get_noise_scaling(rs::ReactionSystem) + return nothing +end + """ ```julia Base.convert(::Type{<:JumpSystem},rs::ReactionSystem; combinatoric_ratelaws=true) @@ -1616,7 +1632,7 @@ end # SDEProblem from AbstractReactionNetwork function DiffEqBase.SDEProblem(rs::ReactionSystem, u0, tspan, p = DiffEqBase.NullParameters(), args...; - noise_scaling = nothing, name = nameof(rs), + noise_scaling = get_noise_scaling(rs), name = nameof(rs), combinatoric_ratelaws = get_combinatoric_ratelaws(rs), include_zero_odes = true, checks = false, check_length = false, From 8b1bdb1b151d1e1813781cbaeccc4803bef82449 Mon Sep 17 00:00:00 2001 From: Torkel Date: Sun, 17 Sep 2023 17:28:00 -0400 Subject: [PATCH 02/41] improved noise scaling --- src/networkapi.jl | 26 +++++++++++++++++++++----- src/reaction_network.jl | 35 ++++++++++++++++++++++++++++++++++- src/reactionsystem.jl | 13 +++++++++++-- 3 files changed, 66 insertions(+), 8 deletions(-) diff --git a/src/networkapi.jl b/src/networkapi.jl index 1926fd34e2..b75e4ac30d 100644 --- a/src/networkapi.jl +++ b/src/networkapi.jl @@ -504,7 +504,9 @@ end # convert symbol of the form :sys.a.b.c to a symbolic a.b.c function _symbol_to_var(sys, sym) - if hasproperty(sys, sym) + if sym isa Num + return sym + elseif hasproperty(sys, sym) var = getproperty(sys, sym, namespace = false) else strs = split(String(sym), "₊") # need to check if this should be split of not!!! @@ -563,25 +565,39 @@ pmap = symmap_to_varmap(sys, [:β => 1.0, :ν => 1.0, :subsys₊k => 1.0]) Notes: - Any `Symbol`, `sym`, within `symmap` must be a valid field of `sys`. i.e. `sys.sym` must be defined. +- The function works even if symbols and symbolics are mixed, e.g. +```julia +sir = @reaction_network sir begin + β, S + I --> 2I + ν, I --> R +end +@unpack S = sir +mixed_u0map = [S => 1.0, :I => 1.0, :R => 1.0] +u0map = symmap_to_varmap(sir, mixed_u0map) +``` + """ function symmap_to_varmap(sys, symmap::Tuple) if all(p -> p isa Pair{Symbol}, symmap) return ((_symbol_to_var(sys, sym) => val for (sym, val) in symmap)...,) else # if not all entries map a symbol to value pass through - return symmap + return symmapAny end end -function symmap_to_varmap(sys, symmap::AbstractArray{Pair{Symbol, T}}) where {T} +function symmap_to_varmap(sys, symmap::AbstractArray{Pair{Any, T}}) where {T} [_symbol_to_var(sys, sym) => val for (sym, val) in symmap] end -function symmap_to_varmap(sys, symmap::Dict{Symbol, T}) where {T} +function symmap_to_varmap(sys, symmap::Dict{Any, T}) where {T} Dict(_symbol_to_var(sys, sym) => val for (sym, val) in symmap) end -# don't permute any other types and let varmap_to_vars handle erroring +# don't permute any other types and let varmap_to_vars handle erroring. +# If all keys are `Num`, conversion not needed. symmap_to_varmap(sys, symmap) = symmap +symmap_to_varmap(sys, symmap::AbstractArray{Pair{Num, T}}) where {T} = symmap +symmap_to_varmap(sys, symmap::Dict{Num, T}) where {T} = symmap #error("symmap_to_varmap requires a Dict, AbstractArray or Tuple to map Symbols to values.") ######################## reaction complexes and reaction rates ############################### diff --git a/src/reaction_network.jl b/src/reaction_network.jl index 753847a6fb..9995bb2c55 100644 --- a/src/reaction_network.jl +++ b/src/reaction_network.jl @@ -369,6 +369,9 @@ function make_reaction_system(ex::Expr; name = :(gensym(:ReactionSystem))) parameters_declared = extract_syms(options, :parameters) variables = extract_syms(options, :variables) + # Handles a (potential) noise scaling parameter. + #nosie_scaling_p_args = handle_noise_scaling_ps!(parameters_declared, options) + # handle independent variables if haskey(options, :ivs) ivs = Tuple(extract_syms(options, :ivs)) @@ -391,7 +394,7 @@ function make_reaction_system(ex::Expr; name = :(gensym(:ReactionSystem))) declared_syms) species = vcat(species_declared, species_extracted) parameters = vcat(parameters_declared, parameters_extracted) - + # Checks for input errors. (sum(length.([reaction_lines, option_lines])) != length(ex.args)) && error("@reaction_network input contain $(length(ex.args) - sum(length.([reaction_lines,option_lines]))) malformed lines.") @@ -414,6 +417,7 @@ function make_reaction_system(ex::Expr; name = :(gensym(:ReactionSystem))) push!(rxexprs.args, get_rxexprs(reaction)) end + # Returns the rephrased expression. quote $ps $ivexpr @@ -556,8 +560,22 @@ function get_pexpr(parameters_extracted, options) pexprs = (haskey(options, :parameters) ? options[:parameters] : (isempty(parameters_extracted) ? :() : :(@parameters))) foreach(p -> push!(pexprs.args, p), parameters_extracted) + append!(pexprs.args, get_noise_scaling_pexpr(options)) pexprs end +# Extracts any decalred nosie scaling parameters. +function get_noise_scaling_pexpr(options) + haskey(options, :noise_scaling_parameters) || return [] + ns_expr = options[:noise_scaling_parameters] + for idx = length(ns_expr.args):-1:3 + if ns_expr.args[idx] isa Symbol + insert!(ns_expr.args, idx+1, :([noisescalingparameter=true])) + elseif (ns_expr.args[idx] isa Expr) && (ns_expr.args[idx].head == :ref) + insert!(ns_expr.args, idx+1, :([noisescalingparameter=true])) + end + end + return ns_expr.args[3:end] +end # Creates the reaction vector declaration statement. function get_rxexprs(rxstruct) @@ -811,6 +829,21 @@ function recursive_expand_functions!(expr::ExprValues) expr end +### Option Handling ### + +# Extracts any potential nosie scaling parameters and add them to teh decalred parameters. +function add_noise_scaling_ps!(parameters_declared, options) + haskey(options, :noise_scaling_parameters) || return + for arg in options[:noise_scaling_parameters].args[end:-1:3] + if arg isa Symbol + push!(parameters_declared, arg) + elseif arg isa Expr + push!(parameters_declared, arg.args[1]) + end + end +end + + # ### Old functions (for deleting). # function get_rx_species_deletethis(rxs, ps) diff --git a/src/reactionsystem.jl b/src/reactionsystem.jl index 7f8691b391..f62b376d1b 100644 --- a/src/reactionsystem.jl +++ b/src/reactionsystem.jl @@ -58,7 +58,7 @@ drop_dynamics(s) = isconstant(s) || isbc(s) || (!isspecies(s)) # Denotes that a parameter controls the scaling of noise in the CLE. struct NoiseScalingParameter end -Symbolics.option_to_metadata_type(::Val{:isnoisescalingparameter}) = NoiseScalingParameter +Symbolics.option_to_metadata_type(::Val{:noisescalingparameter}) = NoiseScalingParameter isnoisescalingparameter(s::Num) = isnoisescalingparameter(MT.value(s)) function isnoisescalingparameter(s) @@ -1546,7 +1546,16 @@ end # Extracts any noise scaling parameters from a reaction system. function get_noise_scaling(rs::ReactionSystem) - return nothing + ns_params = filter(p -> isnoisescalingparameter(p), parameters(rs)) + if isempty(ns_params) + return nothing + elseif length(ns_params) == 1 + return Num(ns_params[1]) + elseif length(ns_params) == length(reactions(rs)) + return ns_params + else + error("The system have $(length(ns_params)) noise scaling parameters. This number should be equal to 0, 1, or the number of reactions ($(length(reactions(rs)))).") + end end """ From 4ed3a25546d55aabf37221f53cef9bbd6f985c0d Mon Sep 17 00:00:00 2001 From: Torkel Date: Mon, 18 Sep 2023 09:56:23 -0400 Subject: [PATCH 03/41] add tests --- test/model_simulation/simulate_SDEs.jl | 55 +++++++++++++++++++++++++- 1 file changed, 54 insertions(+), 1 deletion(-) diff --git a/test/model_simulation/simulate_SDEs.jl b/test/model_simulation/simulate_SDEs.jl index aa1ff5fe80..a2f033702d 100644 --- a/test/model_simulation/simulate_SDEs.jl +++ b/test/model_simulation/simulate_SDEs.jl @@ -121,7 +121,7 @@ end ### Noise Scaling ### -# Tests with a single noise scaling parameter. +# Tests where a single noise scaling parameter is given directly to the SDEProblem. let noise_scaling_network = @reaction_network begin (k1, k2), X1 ↔ X2 end for repeat in 1:5 @@ -144,6 +144,59 @@ let end end +# Tests with multiple noise scaling parameters directly in the macro. +let + noise_scaling_network_1 = @reaction_network begin + @noise_scaling_parameters η1 η2 + (k1, k2), X1 ↔ X2 + end + u0 = [:X1 => 1000.0, :X2 => 3000.0] + sol_1_1 = solve(SDEProblem(noise_scaling_network_1, u0, (0.0, 1000.0), [:k1 => 2.0, :k2 => 0.66, :η1 => 2.0, :η2 => 2.0]), ImplicitEM(); saveat=1.0) + sol_1_2 = solve(SDEProblem(noise_scaling_network_1, u0, (0.0, 1000.0), [:k1 => 2.0, :k2 => 0.66, :η1 => 2.0, :η2 => 0.2]), ImplicitEM(); saveat=1.0) + sol_1_3 = solve(SDEProblem(noise_scaling_network_1, u0, (0.0, 1000.0), [:k1 => 2.0, :k2 => 0.66, :η1 => 0.2, :η2 => 0.2]), ImplicitEM(); saveat=1.0) + @test var(sol_1_1[1,:]) > var(sol_1_2[1,:]) > var(sol_1_3[1,:]) + + noise_scaling_network_2 = @reaction_network begin + @noise_scaling_parameters η[1:2] + (k1, k2), X1 ↔ X2 + end + @unpack η = noise_scaling_network_2 + sol_2_1 = solve(SDEProblem(noise_scaling_network_2, u0, (0.0, 1000.0), [:k1 => 2.0, :k2 => 0.66, η[1] => 2.0, η[2] => 2.0]), ImplicitEM(); saveat=1.0) + sol_2_2 = solve(SDEProblem(noise_scaling_network_2, u0, (0.0, 1000.0), [:k1 => 2.0, :k2 => 0.66, η[1] => 2.0, η[2] => 0.2]), ImplicitEM(); saveat=1.0) + sol_2_3 = solve(SDEProblem(noise_scaling_network_2, u0, (0.0, 1000.0), [:k1 => 2.0, :k2 => 0.66, η[1] => 0.2, η[2] => 0.2]), ImplicitEM(); saveat=1.0) + @test var(sol_2_1[1,:]) > var(sol_2_2[1,:]) > var(sol_2_3[1,:]) +end + +# Tests using default values for nosie scaling. +let + noise_scaling_network = @reaction_network begin + @noise_scaling_parameters η=0.0 + (k1, k2), X1 ↔ X2 + end + u0 = [:X1 => 1100.0, :X2 => 2900.0] + p = [:k1 => 2.0, :k2 => 0.66] + ss = solve(SDEProblem(noise_scaling_network, u0, (0.0, 1000.0), p), ImplicitEM())[end] + @test ss ≈ [1000.0, 3000.0] +end + +# Complicated test with many combinations of options. +let + noise_scaling_network = @reaction_network begin + @parameters k1 par1 [description="Parameter par1"] par2 η1 [noisescalingparameter=true] + @noise_scaling_parameters η2=0.0 [description="Parameter η2"] η3=1.0 η4 [description="Parameter η4"] + (p, d), 0 ↔ X1 + (k1, k2), X1 ↔ X2 + end + @unpack X1, η4 = noise_scaling_network + u0 = [X1 => 500.0, :X2 => 500.0] + p = [p => 20.0, :d => 0.1, :η1 => 0.0, :η3 => 0.0, η4 => 0.0, :k1 => 2.0, :k2 => 2.0, :par1 => 1000.0, :par2 => 1000.0] + + ss = solve(SDEProblem(noise_scaling_network, u0, (0.0, 1000.0), p), ImplicitEM())[end] + ss ≈ [200.0, 200.0] +end + +# Tests that nosie scaling wor + ### Checks Simulations Don't Error ### #Tries to create a large number of problem, ensuring there are no errors (cannot solve as solution likely to go into negatives). From 031f60e00cf55d67fac15599ac777780a8f8527b Mon Sep 17 00:00:00 2001 From: Torkel Date: Mon, 18 Sep 2023 13:43:51 -0400 Subject: [PATCH 04/41] Handle noise parameter default values --- src/networkapi.jl | 8 +++++--- src/reaction_network.jl | 12 +++++++---- test/model_simulation/simulate_SDEs.jl | 17 +++++++++------- .../model_simulation/u0_n_parameter_inputs.jl | 20 +++++++++++++++++++ 4 files changed, 43 insertions(+), 14 deletions(-) diff --git a/src/networkapi.jl b/src/networkapi.jl index b75e4ac30d..ebf2e65ee4 100644 --- a/src/networkapi.jl +++ b/src/networkapi.jl @@ -581,22 +581,24 @@ function symmap_to_varmap(sys, symmap::Tuple) if all(p -> p isa Pair{Symbol}, symmap) return ((_symbol_to_var(sys, sym) => val for (sym, val) in symmap)...,) else # if not all entries map a symbol to value pass through - return symmapAny + return symmap end end -function symmap_to_varmap(sys, symmap::AbstractArray{Pair{Any, T}}) where {T} +function symmap_to_varmap(sys, symmap::AbstractArray{Pair{T, V}}) where {T, V} [_symbol_to_var(sys, sym) => val for (sym, val) in symmap] end -function symmap_to_varmap(sys, symmap::Dict{Any, T}) where {T} +function symmap_to_varmap(sys, symmap::Dict{T, V}) where {T, V} Dict(_symbol_to_var(sys, sym) => val for (sym, val) in symmap) end # don't permute any other types and let varmap_to_vars handle erroring. # If all keys are `Num`, conversion not needed. symmap_to_varmap(sys, symmap) = symmap +symmap_to_varmap(sys, symmap::AbstractArray{Pair{SymbolicUtils.BasicSymbolic{Real}, T}}) where {T} = symmap symmap_to_varmap(sys, symmap::AbstractArray{Pair{Num, T}}) where {T} = symmap +symmap_to_varmap(sys, symmap::Dict{SymbolicUtils.BasicSymbolic{Real}, T}) where {T} = symmap symmap_to_varmap(sys, symmap::Dict{Num, T}) where {T} = symmap #error("symmap_to_varmap requires a Dict, AbstractArray or Tuple to map Symbols to values.") diff --git a/src/reaction_network.jl b/src/reaction_network.jl index 9995bb2c55..8d1869d8aa 100644 --- a/src/reaction_network.jl +++ b/src/reaction_network.jl @@ -568,10 +568,14 @@ function get_noise_scaling_pexpr(options) haskey(options, :noise_scaling_parameters) || return [] ns_expr = options[:noise_scaling_parameters] for idx = length(ns_expr.args):-1:3 - if ns_expr.args[idx] isa Symbol - insert!(ns_expr.args, idx+1, :([noisescalingparameter=true])) - elseif (ns_expr.args[idx] isa Expr) && (ns_expr.args[idx].head == :ref) - insert!(ns_expr.args, idx+1, :([noisescalingparameter=true])) + if (ns_expr.args[idx] isa Symbol) || # Parameter on form η. + (ns_expr.args[idx] isa Expr) && (ns_expr.args[idx].head == :ref) || # Parameter on form η[1:3]. + (ns_expr.args[idx] isa Expr) && (ns_expr.args[idx].head == :(=)) # Parameter on form η=0.1. + if idx < length(ns_expr.args) && (ns_expr.args[idx+1] isa Expr) && (ns_expr.args[idx+1].head == :vect) + push!(ns_expr.args[idx+1].args,:(noisescalingparameter=true)) + else + insert!(ns_expr.args, idx+1, :([noisescalingparameter=true])) + end end end return ns_expr.args[3:end] diff --git a/test/model_simulation/simulate_SDEs.jl b/test/model_simulation/simulate_SDEs.jl index a2f033702d..3cad57f147 100644 --- a/test/model_simulation/simulate_SDEs.jl +++ b/test/model_simulation/simulate_SDEs.jl @@ -167,32 +167,35 @@ let @test var(sol_2_1[1,:]) > var(sol_2_2[1,:]) > var(sol_2_3[1,:]) end -# Tests using default values for nosie scaling. +# Tests using default values for noise scaling. let noise_scaling_network = @reaction_network begin @noise_scaling_parameters η=0.0 (k1, k2), X1 ↔ X2 end - u0 = [:X1 => 1100.0, :X2 => 2900.0] - p = [:k1 => 2.0, :k2 => 0.66] + u0 = [:X1 => 1100.0, :X2 => 3900.0] + p = [:k1 => 2.0, :k2 => 0.5, :η=>0.0] ss = solve(SDEProblem(noise_scaling_network, u0, (0.0, 1000.0), p), ImplicitEM())[end] - @test ss ≈ [1000.0, 3000.0] + @test ss ≈ [1000.0, 4000.0] end # Complicated test with many combinations of options. let noise_scaling_network = @reaction_network begin @parameters k1 par1 [description="Parameter par1"] par2 η1 [noisescalingparameter=true] - @noise_scaling_parameters η2=0.0 [description="Parameter η2"] η3=1.0 η4 [description="Parameter η4"] + @noise_scaling_parameters η2=0.0 [description="Parameter η2"] η3=1.0 η4 (p, d), 0 ↔ X1 (k1, k2), X1 ↔ X2 end - @unpack X1, η4 = noise_scaling_network + @unpack X1, η4, p = noise_scaling_network u0 = [X1 => 500.0, :X2 => 500.0] p = [p => 20.0, :d => 0.1, :η1 => 0.0, :η3 => 0.0, η4 => 0.0, :k1 => 2.0, :k2 => 2.0, :par1 => 1000.0, :par2 => 1000.0] + @test getdescription(parameters(noise_scaling_network)[2]) == "Parameter par1" + @test getdescription(parameters(noise_scaling_network)[8]) == "Parameter η2" + ss = solve(SDEProblem(noise_scaling_network, u0, (0.0, 1000.0), p), ImplicitEM())[end] - ss ≈ [200.0, 200.0] + @test ss ≈ [200.0, 200.0] end # Tests that nosie scaling wor diff --git a/test/model_simulation/u0_n_parameter_inputs.jl b/test/model_simulation/u0_n_parameter_inputs.jl index 1ac44fdca7..8ef750c67c 100644 --- a/test/model_simulation/u0_n_parameter_inputs.jl +++ b/test/model_simulation/u0_n_parameter_inputs.jl @@ -72,3 +72,23 @@ let end end end + +# Tests uding mix of symbols and symbolics in input. +let + test_network = @reaction_network begin + (p1, d1), 0 ↔ X1 + (p2, d2), 0 ↔ X2 + end + @unpack p1, d1, p2, d2, X1, X2 = test_network + u0_1 = [X1 => 0.7, X2 => 3.6] + u0_2 = [:X1 => 0.7, X2 => 3.6] + u0_3 = [:X1 => 0.7, :X2 => 3.6] + p_1 = [p1 => 1.2, d1 => 4.0, p2 => 2.5, d2 =>0.1] + p_2 = [:p1 => 1.2, d1 => 4.0, :p2 => 2.5, d2 =>0.1] + p_3 = [:p1 => 1.2, :d1 => 4.0, :p2 => 2.5, :d2 =>0.1] + + ss_base = solve(ODEProblem(test_network, u0_1, (0.0, 10.0), p_1), Tsit5())[end] + for u0 in [u0_1, u0_2, u0_3], p in [p_1, p_2, p_3] + @test ss_base == solve(ODEProblem(test_network, u0, (0.0, 10.0), p), Tsit5())[end] + end +end \ No newline at end of file From e172c47f7983977e34d004aebd1025b1512346d3 Mon Sep 17 00:00:00 2001 From: Torkel Date: Mon, 18 Sep 2023 16:26:03 -0400 Subject: [PATCH 05/41] clean --- src/reaction_network.jl | 14 -------------- 1 file changed, 14 deletions(-) diff --git a/src/reaction_network.jl b/src/reaction_network.jl index 8d1869d8aa..e4f07509a4 100644 --- a/src/reaction_network.jl +++ b/src/reaction_network.jl @@ -833,20 +833,6 @@ function recursive_expand_functions!(expr::ExprValues) expr end -### Option Handling ### - -# Extracts any potential nosie scaling parameters and add them to teh decalred parameters. -function add_noise_scaling_ps!(parameters_declared, options) - haskey(options, :noise_scaling_parameters) || return - for arg in options[:noise_scaling_parameters].args[end:-1:3] - if arg isa Symbol - push!(parameters_declared, arg) - elseif arg isa Expr - push!(parameters_declared, arg.args[1]) - end - end -end - # ### Old functions (for deleting). From a5f6cf8ffd9e48b47176c0f5b29eb989d19291be Mon Sep 17 00:00:00 2001 From: Torkel Date: Mon, 18 Sep 2023 18:45:39 -0400 Subject: [PATCH 06/41] update --- src/networkapi.jl | 11 ++--- src/reaction_network.jl | 7 ++- src/reactionsystem.jl | 7 ++- test/model_simulation/simulate_SDEs.jl | 44 ++++--------------- .../model_simulation/u0_n_parameter_inputs.jl | 20 --------- 5 files changed, 23 insertions(+), 66 deletions(-) diff --git a/src/networkapi.jl b/src/networkapi.jl index ebf2e65ee4..ea50cbee06 100644 --- a/src/networkapi.jl +++ b/src/networkapi.jl @@ -585,21 +585,16 @@ function symmap_to_varmap(sys, symmap::Tuple) end end -function symmap_to_varmap(sys, symmap::AbstractArray{Pair{T, V}}) where {T, V} +function symmap_to_varmap(sys, symmap::AbstractArray{Pair{Symbol, T}}) where {T} [_symbol_to_var(sys, sym) => val for (sym, val) in symmap] end -function symmap_to_varmap(sys, symmap::Dict{T, V}) where {T, V} +function symmap_to_varmap(sys, symmap::Dict{Symbol, T}) where {T} Dict(_symbol_to_var(sys, sym) => val for (sym, val) in symmap) end -# don't permute any other types and let varmap_to_vars handle erroring. -# If all keys are `Num`, conversion not needed. +# don't permute any other types and let varmap_to_vars handle erroring symmap_to_varmap(sys, symmap) = symmap -symmap_to_varmap(sys, symmap::AbstractArray{Pair{SymbolicUtils.BasicSymbolic{Real}, T}}) where {T} = symmap -symmap_to_varmap(sys, symmap::AbstractArray{Pair{Num, T}}) where {T} = symmap -symmap_to_varmap(sys, symmap::Dict{SymbolicUtils.BasicSymbolic{Real}, T}) where {T} = symmap -symmap_to_varmap(sys, symmap::Dict{Num, T}) where {T} = symmap #error("symmap_to_varmap requires a Dict, AbstractArray or Tuple to map Symbols to values.") ######################## reaction complexes and reaction rates ############################### diff --git a/src/reaction_network.jl b/src/reaction_network.jl index e4f07509a4..8be05fed33 100644 --- a/src/reaction_network.jl +++ b/src/reaction_network.jl @@ -370,7 +370,7 @@ function make_reaction_system(ex::Expr; name = :(gensym(:ReactionSystem))) variables = extract_syms(options, :variables) # Handles a (potential) noise scaling parameter. - #nosie_scaling_p_args = handle_noise_scaling_ps!(parameters_declared, options) + #noise_scaling_p_args = handle_noise_scaling_ps!(parameters_declared, options) # handle independent variables if haskey(options, :ivs) @@ -417,6 +417,9 @@ function make_reaction_system(ex::Expr; name = :(gensym(:ReactionSystem))) push!(rxexprs.args, get_rxexprs(reaction)) end + # Adds potentital noise scaling parameters. + #append!(ps.args, noise_scaling_p_args) + # Returns the rephrased expression. quote $ps @@ -563,7 +566,7 @@ function get_pexpr(parameters_extracted, options) append!(pexprs.args, get_noise_scaling_pexpr(options)) pexprs end -# Extracts any decalred nosie scaling parameters. +# Extracts any decalred noise scaling parameters. function get_noise_scaling_pexpr(options) haskey(options, :noise_scaling_parameters) || return [] ns_expr = options[:noise_scaling_parameters] diff --git a/src/reactionsystem.jl b/src/reactionsystem.jl index f62b376d1b..10c098d55f 100644 --- a/src/reactionsystem.jl +++ b/src/reactionsystem.jl @@ -1506,6 +1506,7 @@ function Base.convert(::Type{<:SDESystem}, rs::ReactionSystem; default_u0 = Dict(), default_p = Dict(), defaults = _merge(Dict(default_u0), Dict(default_p)), kwargs...) spatial_convert_err(rs::ReactionSystem, SDESystem) + isnothing(noise_scaling) && (noise_scaling = get_noise_scaling(rs)) # Required before final deprication of previous noise scaling functionality. flatrs = Catalyst.flatten(rs) error_if_constraints(SDESystem, flatrs) @@ -1641,11 +1642,15 @@ end # SDEProblem from AbstractReactionNetwork function DiffEqBase.SDEProblem(rs::ReactionSystem, u0, tspan, p = DiffEqBase.NullParameters(), args...; - noise_scaling = get_noise_scaling(rs), name = nameof(rs), + noise_scaling = nothing, name = nameof(rs), combinatoric_ratelaws = get_combinatoric_ratelaws(rs), include_zero_odes = true, checks = false, check_length = false, remove_conserved = false, kwargs...) + if !isnothing(noise_scaling) + any(isnoisescalingparameter.(parameters(rs))) && error("You have declared some paraemters as noise scaling parameters, and also given a \"noise_scaling\" argument to SDEProblem. Please remove the \"noise_scaling\", as this way of scaling CLE noise is being depricated.") + @warn "Passing noise scaling input into SDEProblem will be depricated. New standard is to declare one (or several) paraemter as noise scaling parameters when the ReactionSystem is created. Please read https://docs.sciml.ai/Catalyst/stable/catalyst_applications/advanced_simulations/#Scaling-the-noise-magnitude-in-the-chemical-Langevin-equations." + end u0map = symmap_to_varmap(rs, u0) pmap = symmap_to_varmap(rs, p) sde_sys = convert(SDESystem, rs; noise_scaling, name, combinatoric_ratelaws, diff --git a/test/model_simulation/simulate_SDEs.jl b/test/model_simulation/simulate_SDEs.jl index 3cad57f147..b51ee1cab5 100644 --- a/test/model_simulation/simulate_SDEs.jl +++ b/test/model_simulation/simulate_SDEs.jl @@ -121,29 +121,6 @@ end ### Noise Scaling ### -# Tests where a single noise scaling parameter is given directly to the SDEProblem. -let - noise_scaling_network = @reaction_network begin (k1, k2), X1 ↔ X2 end - for repeat in 1:5 - p = 1.0 .+ rand(rng, 2) - u0 = 10000 * (1.0 .+ rand(rng, 2)) - sol001 = solve(SDEProblem(noise_scaling_network, u0, (0.0, 1000.0), vcat(p, 0.01), - noise_scaling = (@variables η1)[1]), ImplicitEM()) - sol01 = solve(SDEProblem(noise_scaling_network, u0, (0.0, 1000.0), vcat(p, 0.1), - noise_scaling = (@variables η1)[1]), ImplicitEM()) - sol1 = solve(SDEProblem(noise_scaling_network, u0, (0.0, 1000.0), vcat(p, 1.0), - noise_scaling = (@variables η2)[1]), ImplicitEM()) - sol10 = solve(SDEProblem(noise_scaling_network, u0, (0.0, 1000.0), vcat(p, 10.0), - noise_scaling = (@variables η3)[1]), ImplicitEM()) - @test 2 * std(first.(sol001.u)[100:end]) < std(first.(sol01.u)[100:end]) - @test 2 * std(last.(sol001.u)[100:end]) < std(last.(sol01.u)[100:end]) - @test 2 * std(first.(sol01.u)[100:end]) < std(first.(sol1.u)[100:end]) - @test 2 * std(last.(sol01.u)[100:end]) < std(last.(sol1.u)[100:end]) - @test 2 * std(first.(sol1.u)[100:end]) < std(first.(sol10.u)[100:end]) - @test 2 * std(last.(sol1.u)[100:end]) < std(last.(sol10.u)[100:end]) - end -end - # Tests with multiple noise scaling parameters directly in the macro. let noise_scaling_network_1 = @reaction_network begin @@ -160,10 +137,10 @@ let @noise_scaling_parameters η[1:2] (k1, k2), X1 ↔ X2 end - @unpack η = noise_scaling_network_2 - sol_2_1 = solve(SDEProblem(noise_scaling_network_2, u0, (0.0, 1000.0), [:k1 => 2.0, :k2 => 0.66, η[1] => 2.0, η[2] => 2.0]), ImplicitEM(); saveat=1.0) - sol_2_2 = solve(SDEProblem(noise_scaling_network_2, u0, (0.0, 1000.0), [:k1 => 2.0, :k2 => 0.66, η[1] => 2.0, η[2] => 0.2]), ImplicitEM(); saveat=1.0) - sol_2_3 = solve(SDEProblem(noise_scaling_network_2, u0, (0.0, 1000.0), [:k1 => 2.0, :k2 => 0.66, η[1] => 0.2, η[2] => 0.2]), ImplicitEM(); saveat=1.0) + @unpack k1, k2, η = noise_scaling_network_2 + sol_2_1 = solve(SDEProblem(noise_scaling_network_2, u0, (0.0, 1000.0), [k1 => 2.0, k2 => 0.66, η[1] => 2.0, η[2] => 2.0]), ImplicitEM(); saveat=1.0) + sol_2_2 = solve(SDEProblem(noise_scaling_network_2, u0, (0.0, 1000.0), [k1 => 2.0, k2 => 0.66, η[1] => 2.0, η[2] => 0.2]), ImplicitEM(); saveat=1.0) + sol_2_3 = solve(SDEProblem(noise_scaling_network_2, u0, (0.0, 1000.0), [k1 => 2.0, k2 => 0.66, η[1] => 0.2, η[2] => 0.2]), ImplicitEM(); saveat=1.0) @test var(sol_2_1[1,:]) > var(sol_2_2[1,:]) > var(sol_2_3[1,:]) end @@ -175,8 +152,7 @@ let end u0 = [:X1 => 1100.0, :X2 => 3900.0] p = [:k1 => 2.0, :k2 => 0.5, :η=>0.0] - ss = solve(SDEProblem(noise_scaling_network, u0, (0.0, 1000.0), p), ImplicitEM())[end] - @test ss ≈ [1000.0, 4000.0] + @test SDEProblem(noise_scaling_network, u0, (0.0, 1000.0), p)[:η] == 0.0 end # Complicated test with many combinations of options. @@ -187,18 +163,16 @@ let (p, d), 0 ↔ X1 (k1, k2), X1 ↔ X2 end - @unpack X1, η4, p = noise_scaling_network - u0 = [X1 => 500.0, :X2 => 500.0] - p = [p => 20.0, :d => 0.1, :η1 => 0.0, :η3 => 0.0, η4 => 0.0, :k1 => 2.0, :k2 => 2.0, :par1 => 1000.0, :par2 => 1000.0] + u0 = [:X1 => 500.0, :X2 => 500.0] + p = [:p => 20.0, :d => 0.1, :η1 => 0.0, :η3 => 0.0, :η4 => 0.0, :k1 => 2.0, :k2 => 2.0, :par1 => 1000.0, :par2 => 1000.0] @test getdescription(parameters(noise_scaling_network)[2]) == "Parameter par1" @test getdescription(parameters(noise_scaling_network)[8]) == "Parameter η2" - ss = solve(SDEProblem(noise_scaling_network, u0, (0.0, 1000.0), p), ImplicitEM())[end] - @test ss ≈ [200.0, 200.0] + sprob = SDEProblem(noise_scaling_network, u0, (0.0, 1000.0), p) + @test sprob[:η1] == sprob[:η2] == sprob[:η3] == sprob[:η4] == 0.0 end -# Tests that nosie scaling wor ### Checks Simulations Don't Error ### diff --git a/test/model_simulation/u0_n_parameter_inputs.jl b/test/model_simulation/u0_n_parameter_inputs.jl index 8ef750c67c..6b3a0adf03 100644 --- a/test/model_simulation/u0_n_parameter_inputs.jl +++ b/test/model_simulation/u0_n_parameter_inputs.jl @@ -71,24 +71,4 @@ let @test discrete_probs[1].u0 == discrete_probs[i].u0 end end -end - -# Tests uding mix of symbols and symbolics in input. -let - test_network = @reaction_network begin - (p1, d1), 0 ↔ X1 - (p2, d2), 0 ↔ X2 - end - @unpack p1, d1, p2, d2, X1, X2 = test_network - u0_1 = [X1 => 0.7, X2 => 3.6] - u0_2 = [:X1 => 0.7, X2 => 3.6] - u0_3 = [:X1 => 0.7, :X2 => 3.6] - p_1 = [p1 => 1.2, d1 => 4.0, p2 => 2.5, d2 =>0.1] - p_2 = [:p1 => 1.2, d1 => 4.0, :p2 => 2.5, d2 =>0.1] - p_3 = [:p1 => 1.2, :d1 => 4.0, :p2 => 2.5, :d2 =>0.1] - - ss_base = solve(ODEProblem(test_network, u0_1, (0.0, 10.0), p_1), Tsit5())[end] - for u0 in [u0_1, u0_2, u0_3], p in [p_1, p_2, p_3] - @test ss_base == solve(ODEProblem(test_network, u0, (0.0, 10.0), p), Tsit5())[end] - end end \ No newline at end of file From cd3a3aa29dfd4dc0c8408419cb2fe6538167fd30 Mon Sep 17 00:00:00 2001 From: Torkel Date: Mon, 18 Sep 2023 19:59:21 -0400 Subject: [PATCH 07/41] Update src/networkapi.jl Co-authored-by: Sam Isaacson --- src/networkapi.jl | 11 ----------- 1 file changed, 11 deletions(-) diff --git a/src/networkapi.jl b/src/networkapi.jl index ea50cbee06..4868b5d055 100644 --- a/src/networkapi.jl +++ b/src/networkapi.jl @@ -565,17 +565,6 @@ pmap = symmap_to_varmap(sys, [:β => 1.0, :ν => 1.0, :subsys₊k => 1.0]) Notes: - Any `Symbol`, `sym`, within `symmap` must be a valid field of `sys`. i.e. `sys.sym` must be defined. -- The function works even if symbols and symbolics are mixed, e.g. -```julia -sir = @reaction_network sir begin - β, S + I --> 2I - ν, I --> R -end -@unpack S = sir -mixed_u0map = [S => 1.0, :I => 1.0, :R => 1.0] -u0map = symmap_to_varmap(sir, mixed_u0map) -``` - """ function symmap_to_varmap(sys, symmap::Tuple) if all(p -> p isa Pair{Symbol}, symmap) From 516e7cce3ac155b95d57a065c6aece2408af075b Mon Sep 17 00:00:00 2001 From: Torkel Date: Mon, 18 Sep 2023 19:59:33 -0400 Subject: [PATCH 08/41] Update src/networkapi.jl Co-authored-by: Sam Isaacson --- src/networkapi.jl | 3 --- 1 file changed, 3 deletions(-) diff --git a/src/networkapi.jl b/src/networkapi.jl index 4868b5d055..a12c0ed0d7 100644 --- a/src/networkapi.jl +++ b/src/networkapi.jl @@ -504,9 +504,6 @@ end # convert symbol of the form :sys.a.b.c to a symbolic a.b.c function _symbol_to_var(sys, sym) - if sym isa Num - return sym - elseif hasproperty(sys, sym) var = getproperty(sys, sym, namespace = false) else strs = split(String(sym), "₊") # need to check if this should be split of not!!! From a32c88596623957522d9f05d0965505044a4e4db Mon Sep 17 00:00:00 2001 From: Torkel Date: Mon, 18 Sep 2023 19:59:47 -0400 Subject: [PATCH 09/41] Update src/reaction_network.jl Co-authored-by: Sam Isaacson --- src/reaction_network.jl | 3 --- 1 file changed, 3 deletions(-) diff --git a/src/reaction_network.jl b/src/reaction_network.jl index 8be05fed33..5ecf567ca5 100644 --- a/src/reaction_network.jl +++ b/src/reaction_network.jl @@ -369,9 +369,6 @@ function make_reaction_system(ex::Expr; name = :(gensym(:ReactionSystem))) parameters_declared = extract_syms(options, :parameters) variables = extract_syms(options, :variables) - # Handles a (potential) noise scaling parameter. - #noise_scaling_p_args = handle_noise_scaling_ps!(parameters_declared, options) - # handle independent variables if haskey(options, :ivs) ivs = Tuple(extract_syms(options, :ivs)) From 82036876af11881ccd449207d7d41fdbd6e2e7de Mon Sep 17 00:00:00 2001 From: Torkel Date: Mon, 18 Sep 2023 19:59:59 -0400 Subject: [PATCH 10/41] Update src/reaction_network.jl Co-authored-by: Sam Isaacson --- src/reaction_network.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/src/reaction_network.jl b/src/reaction_network.jl index 5ecf567ca5..6a5908bac4 100644 --- a/src/reaction_network.jl +++ b/src/reaction_network.jl @@ -414,7 +414,6 @@ function make_reaction_system(ex::Expr; name = :(gensym(:ReactionSystem))) push!(rxexprs.args, get_rxexprs(reaction)) end - # Adds potentital noise scaling parameters. #append!(ps.args, noise_scaling_p_args) # Returns the rephrased expression. From 131c7ad629a17e3eaddb50ba6089a33df2bf75fc Mon Sep 17 00:00:00 2001 From: Torkel Date: Mon, 18 Sep 2023 22:29:50 -0400 Subject: [PATCH 11/41] update --- src/Catalyst.jl | 2 +- src/networkapi.jl | 1 + src/reaction_network.jl | 30 +++++++++----------------- src/reactionsystem.jl | 17 +++++++-------- test/model_simulation/simulate_SDEs.jl | 2 +- 5 files changed, 21 insertions(+), 31 deletions(-) diff --git a/src/Catalyst.jl b/src/Catalyst.jl index 6bd26b35c9..3f0573fd6a 100644 --- a/src/Catalyst.jl +++ b/src/Catalyst.jl @@ -77,7 +77,7 @@ export ODEProblem, const ExprValues = Union{Expr, Symbol, Float64, Int, Bool} include("expression_utils.jl") include("reaction_network.jl") -export @reaction_network, @add_reactions, @reaction, @species +export @reaction_network, @add_reactions, @reaction, @species, @noise_scaling_parameters # registers CRN specific functions using Symbolics.jl include("registered_functions.jl") diff --git a/src/networkapi.jl b/src/networkapi.jl index a12c0ed0d7..1926fd34e2 100644 --- a/src/networkapi.jl +++ b/src/networkapi.jl @@ -504,6 +504,7 @@ end # convert symbol of the form :sys.a.b.c to a symbolic a.b.c function _symbol_to_var(sys, sym) + if hasproperty(sys, sym) var = getproperty(sys, sym, namespace = false) else strs = split(String(sym), "₊") # need to check if this should be split of not!!! diff --git a/src/reaction_network.jl b/src/reaction_network.jl index 6a5908bac4..aa2c2638b7 100644 --- a/src/reaction_network.jl +++ b/src/reaction_network.jl @@ -153,6 +153,7 @@ macro species(ex...) # put back the vector of the new species symbols vars.args[idx] = lastarg + println(vars) esc(vars) end @@ -369,6 +370,9 @@ function make_reaction_system(ex::Expr; name = :(gensym(:ReactionSystem))) parameters_declared = extract_syms(options, :parameters) variables = extract_syms(options, :variables) + # Handles noise scaling parameters. + noise_scaling_pexpr = handle_noise_scaling_parameters!(parameters_declared, options) + # handle independent variables if haskey(options, :ivs) ivs = Tuple(extract_syms(options, :ivs)) @@ -405,6 +409,7 @@ function make_reaction_system(ex::Expr; name = :(gensym(:ReactionSystem))) sexprs = get_sexpr(species_extracted, options; iv_symbols = ivs) vexprs = haskey(options, :variables) ? options[:variables] : :() pexprs = get_pexpr(parameters_extracted, options) + ns_ps, ns_pssym = scalarize_macro(haskey(options, :noise_scaling_parameters), noise_scaling_pexpr, "ns_ps") ps, pssym = scalarize_macro(!isempty(parameters), pexprs, "ps") vars, varssym = scalarize_macro(!isempty(variables), vexprs, "vars") sps, spssym = scalarize_macro(!isempty(species), sexprs, "specs") @@ -414,11 +419,14 @@ function make_reaction_system(ex::Expr; name = :(gensym(:ReactionSystem))) push!(rxexprs.args, get_rxexprs(reaction)) end - #append!(ps.args, noise_scaling_p_args) - + println(ps) + println(ns_ps) + println(sps) + # Returns the rephrased expression. quote $ps + $ns_ps $ivexpr $vars $sps @@ -559,26 +567,8 @@ function get_pexpr(parameters_extracted, options) pexprs = (haskey(options, :parameters) ? options[:parameters] : (isempty(parameters_extracted) ? :() : :(@parameters))) foreach(p -> push!(pexprs.args, p), parameters_extracted) - append!(pexprs.args, get_noise_scaling_pexpr(options)) pexprs end -# Extracts any decalred noise scaling parameters. -function get_noise_scaling_pexpr(options) - haskey(options, :noise_scaling_parameters) || return [] - ns_expr = options[:noise_scaling_parameters] - for idx = length(ns_expr.args):-1:3 - if (ns_expr.args[idx] isa Symbol) || # Parameter on form η. - (ns_expr.args[idx] isa Expr) && (ns_expr.args[idx].head == :ref) || # Parameter on form η[1:3]. - (ns_expr.args[idx] isa Expr) && (ns_expr.args[idx].head == :(=)) # Parameter on form η=0.1. - if idx < length(ns_expr.args) && (ns_expr.args[idx+1] isa Expr) && (ns_expr.args[idx+1].head == :vect) - push!(ns_expr.args[idx+1].args,:(noisescalingparameter=true)) - else - insert!(ns_expr.args, idx+1, :([noisescalingparameter=true])) - end - end - end - return ns_expr.args[3:end] -end # Creates the reaction vector declaration statement. function get_rxexprs(rxstruct) diff --git a/src/reactionsystem.jl b/src/reactionsystem.jl index 10c098d55f..efc92d02eb 100644 --- a/src/reactionsystem.jl +++ b/src/reactionsystem.jl @@ -58,10 +58,10 @@ drop_dynamics(s) = isconstant(s) || isbc(s) || (!isspecies(s)) # Denotes that a parameter controls the scaling of noise in the CLE. struct NoiseScalingParameter end -Symbolics.option_to_metadata_type(::Val{:noisescalingparameter}) = NoiseScalingParameter +Symbolics.option_to_metadata_type(::Val{:noise_scaling_parameter}) = NoiseScalingParameter -isnoisescalingparameter(s::Num) = isnoisescalingparameter(MT.value(s)) -function isnoisescalingparameter(s) +is_noise_scaling_parameter(s::Num) = is_noise_scaling_parameter(MT.value(s)) +function is_noise_scaling_parameter(s) MT.getmetadata(s, NoiseScalingParameter, false) end @@ -1506,10 +1506,13 @@ function Base.convert(::Type{<:SDESystem}, rs::ReactionSystem; default_u0 = Dict(), default_p = Dict(), defaults = _merge(Dict(default_u0), Dict(default_p)), kwargs...) spatial_convert_err(rs::ReactionSystem, SDESystem) - isnothing(noise_scaling) && (noise_scaling = get_noise_scaling(rs)) # Required before final deprication of previous noise scaling functionality. flatrs = Catalyst.flatten(rs) error_if_constraints(SDESystem, flatrs) + if any(isnoisescalingparameter, get_ps(flatrs)) + any(is_noise_scaling_parameter.(parameters(rs))) && error("You have declared some paraemters as noise scaling parameters, and also given a \"noise_scaling\" argument to SDEProblem. Please remove the \"noise_scaling\", as this way of scaling CLE noise is being depricated.") + @warn "Passing noise scaling input into SDEProblem will be deprecated. New standard is to declare one (or several) paraemter as noise scaling parameters when the ReactionSystem is created. Please read https://docs.sciml.ai/Catalyst/stable/catalyst_applications/advanced_simulations/#Scaling-the-noise-magnitude-in-the-chemical-Langevin-equations." + end if noise_scaling isa AbstractArray (length(noise_scaling) != numreactions(flatrs)) && @@ -1547,7 +1550,7 @@ end # Extracts any noise scaling parameters from a reaction system. function get_noise_scaling(rs::ReactionSystem) - ns_params = filter(p -> isnoisescalingparameter(p), parameters(rs)) + ns_params = filter(p -> is_noise_scaling_parameter(p), parameters(rs)) if isempty(ns_params) return nothing elseif length(ns_params) == 1 @@ -1647,10 +1650,6 @@ function DiffEqBase.SDEProblem(rs::ReactionSystem, u0, tspan, include_zero_odes = true, checks = false, check_length = false, remove_conserved = false, kwargs...) - if !isnothing(noise_scaling) - any(isnoisescalingparameter.(parameters(rs))) && error("You have declared some paraemters as noise scaling parameters, and also given a \"noise_scaling\" argument to SDEProblem. Please remove the \"noise_scaling\", as this way of scaling CLE noise is being depricated.") - @warn "Passing noise scaling input into SDEProblem will be depricated. New standard is to declare one (or several) paraemter as noise scaling parameters when the ReactionSystem is created. Please read https://docs.sciml.ai/Catalyst/stable/catalyst_applications/advanced_simulations/#Scaling-the-noise-magnitude-in-the-chemical-Langevin-equations." - end u0map = symmap_to_varmap(rs, u0) pmap = symmap_to_varmap(rs, p) sde_sys = convert(SDESystem, rs; noise_scaling, name, combinatoric_ratelaws, diff --git a/test/model_simulation/simulate_SDEs.jl b/test/model_simulation/simulate_SDEs.jl index b51ee1cab5..6480382371 100644 --- a/test/model_simulation/simulate_SDEs.jl +++ b/test/model_simulation/simulate_SDEs.jl @@ -158,7 +158,7 @@ end # Complicated test with many combinations of options. let noise_scaling_network = @reaction_network begin - @parameters k1 par1 [description="Parameter par1"] par2 η1 [noisescalingparameter=true] + @parameters k1 par1 [description="Parameter par1"] par2 η1 [noise_scaling_parameter=true] @noise_scaling_parameters η2=0.0 [description="Parameter η2"] η3=1.0 η4 (p, d), 0 ↔ X1 (k1, k2), X1 ↔ X2 From 228eb05c3ea50c1972e056858e8186986325f294 Mon Sep 17 00:00:00 2001 From: Torkel Date: Mon, 18 Sep 2023 22:32:44 -0400 Subject: [PATCH 12/41] change --- src/reaction_network.jl | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/reaction_network.jl b/src/reaction_network.jl index aa2c2638b7..c4c0aaf71b 100644 --- a/src/reaction_network.jl +++ b/src/reaction_network.jl @@ -727,8 +727,9 @@ function extract_metadata(metadata_line::Expr) return metadata end -### DSL Options Handling ### -# Most options handled in previous sections, when code re-organised, these should ideally be moved to the same place. +# macro, similar to @parameters, but all paraemters become noise scaling parameters. +macro noise_scaling_parameters(ex...) + vars = Symbolics._parse_vars(:parameters, Real, ex) # Reads the observables options. Outputs an expression ofr creating the obervable variables, and a vector of observable equations. function read_observed_options(options, species_n_vars_declared, ivs_sorted) From 98a051181fa80905456979d382b7e1050a02ea17 Mon Sep 17 00:00:00 2001 From: Torkel Date: Mon, 18 Sep 2023 22:36:51 -0400 Subject: [PATCH 13/41] update --- src/reaction_network.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/src/reaction_network.jl b/src/reaction_network.jl index c4c0aaf71b..a5a133e839 100644 --- a/src/reaction_network.jl +++ b/src/reaction_network.jl @@ -153,7 +153,6 @@ macro species(ex...) # put back the vector of the new species symbols vars.args[idx] = lastarg - println(vars) esc(vars) end From 4c85f2c4076dd71c3e2a5e822c57a215c8851428 Mon Sep 17 00:00:00 2001 From: Torkel Date: Mon, 25 Sep 2023 17:43:52 -0400 Subject: [PATCH 14/41] fix error --- src/reaction_network.jl | 6 +----- src/reactionsystem.jl | 11 +++++++---- test/model_simulation/simulate_SDEs.jl | 13 +++++++++++++ 3 files changed, 21 insertions(+), 9 deletions(-) diff --git a/src/reaction_network.jl b/src/reaction_network.jl index a5a133e839..9aa2c10482 100644 --- a/src/reaction_network.jl +++ b/src/reaction_network.jl @@ -418,10 +418,6 @@ function make_reaction_system(ex::Expr; name = :(gensym(:ReactionSystem))) push!(rxexprs.args, get_rxexprs(reaction)) end - println(ps) - println(ns_ps) - println(sps) - # Returns the rephrased expression. quote $ps @@ -728,7 +724,7 @@ end # macro, similar to @parameters, but all paraemters become noise scaling parameters. macro noise_scaling_parameters(ex...) - vars = Symbolics._parse_vars(:parameters, Real, ex) + vars = Symbolics._parse_vars(:parameters, Real, ex, toparam) # Reads the observables options. Outputs an expression ofr creating the obervable variables, and a vector of observable equations. function read_observed_options(options, species_n_vars_declared, ivs_sorted) diff --git a/src/reactionsystem.jl b/src/reactionsystem.jl index efc92d02eb..bfd52502a3 100644 --- a/src/reactionsystem.jl +++ b/src/reactionsystem.jl @@ -1509,11 +1509,8 @@ function Base.convert(::Type{<:SDESystem}, rs::ReactionSystem; flatrs = Catalyst.flatten(rs) error_if_constraints(SDESystem, flatrs) - if any(isnoisescalingparameter, get_ps(flatrs)) - any(is_noise_scaling_parameter.(parameters(rs))) && error("You have declared some paraemters as noise scaling parameters, and also given a \"noise_scaling\" argument to SDEProblem. Please remove the \"noise_scaling\", as this way of scaling CLE noise is being depricated.") - @warn "Passing noise scaling input into SDEProblem will be deprecated. New standard is to declare one (or several) paraemter as noise scaling parameters when the ReactionSystem is created. Please read https://docs.sciml.ai/Catalyst/stable/catalyst_applications/advanced_simulations/#Scaling-the-noise-magnitude-in-the-chemical-Langevin-equations." - end + isnothing(noise_scaling) && (noise_scaling = get_noise_scaling(rs)) # Required until passing nosie into SDEProblem can be depricated. if noise_scaling isa AbstractArray (length(noise_scaling) != numreactions(flatrs)) && error("The number of elements in 'noise_scaling' must be equal " * @@ -1650,6 +1647,12 @@ function DiffEqBase.SDEProblem(rs::ReactionSystem, u0, tspan, include_zero_odes = true, checks = false, check_length = false, remove_conserved = false, kwargs...) + + if !isnothing(noise_scaling) + !isnothing(get_noise_scaling(rs)) && error("You have declared some parameters as noise scaling parameters, and also given a \"noise_scaling\" argument to SDEProblem. Please remove the \"noise_scaling\", as this way of scaling CLE noise is being depricated.") + @warn "Passing noise scaling input into SDEProblem will be deprecated. New standard is to declare one (or several) paraemter as noise scaling parameters when the ReactionSystem is created. Please read https://docs.sciml.ai/Catalyst/stable/catalyst_applications/advanced_simulations/#Scaling-the-noise-magnitude-in-the-chemical-Langevin-equations." + end + u0map = symmap_to_varmap(rs, u0) pmap = symmap_to_varmap(rs, p) sde_sys = convert(SDESystem, rs; noise_scaling, name, combinatoric_ratelaws, diff --git a/test/model_simulation/simulate_SDEs.jl b/test/model_simulation/simulate_SDEs.jl index 6480382371..8424d0d9fc 100644 --- a/test/model_simulation/simulate_SDEs.jl +++ b/test/model_simulation/simulate_SDEs.jl @@ -119,6 +119,19 @@ let end end +### Checks Simulations Don't Error ### + +#Tries to create a large number of problem, ensuring there are no errors (cannot solve as solution likely to go into negatives). +let + for reaction_network in reaction_networks_all + for factor in [1e-2, 1e-1, 1e0, 1e1] + u0 = factor * rand(rng, length(states(reaction_network))) + p = factor * rand(rng, length(parameters(reaction_network))) + prob = SDEProblem(reaction_network, u0, (0.0, 1.0), p) + end + end +end + ### Noise Scaling ### # Tests with multiple noise scaling parameters directly in the macro. From 68465d89da07a980c946864e184db413a19164a2 Mon Sep 17 00:00:00 2001 From: Torkel Date: Mon, 25 Sep 2023 20:51:11 -0400 Subject: [PATCH 15/41] update --- src/reaction_network.jl | 2 +- src/reactionsystem.jl | 17 +++++++++-------- 2 files changed, 10 insertions(+), 9 deletions(-) diff --git a/src/reaction_network.jl b/src/reaction_network.jl index 9aa2c10482..86bf42565a 100644 --- a/src/reaction_network.jl +++ b/src/reaction_network.jl @@ -410,6 +410,7 @@ function make_reaction_system(ex::Expr; name = :(gensym(:ReactionSystem))) pexprs = get_pexpr(parameters_extracted, options) ns_ps, ns_pssym = scalarize_macro(haskey(options, :noise_scaling_parameters), noise_scaling_pexpr, "ns_ps") ps, pssym = scalarize_macro(!isempty(parameters), pexprs, "ps") + rs_p_syms = haskey(options, :noise_scaling_parameters) ? :(union($pssym, $ns_pssym)) : pssym vars, varssym = scalarize_macro(!isempty(variables), vexprs, "vars") sps, spssym = scalarize_macro(!isempty(species), sexprs, "specs") comps, compssym = scalarize_macro(!isempty(compound_species), compound_expr, "comps") @@ -417,7 +418,6 @@ function make_reaction_system(ex::Expr; name = :(gensym(:ReactionSystem))) for reaction in reactions push!(rxexprs.args, get_rxexprs(reaction)) end - # Returns the rephrased expression. quote $ps diff --git a/src/reactionsystem.jl b/src/reactionsystem.jl index bfd52502a3..aec89a9ef3 100644 --- a/src/reactionsystem.jl +++ b/src/reactionsystem.jl @@ -1500,7 +1500,7 @@ Notes: differential equations. """ function Base.convert(::Type{<:SDESystem}, rs::ReactionSystem; - noise_scaling = get_noise_scaling(rs), name = nameof(rs), + noise_scaling = nothing, name = nameof(rs), combinatoric_ratelaws = get_combinatoric_ratelaws(rs), include_zero_odes = true, checks = false, remove_conserved = false, default_u0 = Dict(), default_p = Dict(), defaults = _merge(Dict(default_u0), Dict(default_p)), @@ -1510,7 +1510,13 @@ function Base.convert(::Type{<:SDESystem}, rs::ReactionSystem; flatrs = Catalyst.flatten(rs) error_if_constraints(SDESystem, flatrs) - isnothing(noise_scaling) && (noise_scaling = get_noise_scaling(rs)) # Required until passing nosie into SDEProblem can be depricated. + # Required until passing nosie into SDEProblem can be depricated. When properly deprecated, change the kwarg to noise_scaling = get_noise_scaling(rs). + if !isnothing(noise_scaling) + !isnothing(get_noise_scaling(rs)) && error("You have declared some parameters as noise scaling parameters, and also given a \"noise_scaling\" argument to SDEProblem. Please remove the \"noise_scaling\", as this way of scaling CLE noise is being depricated.") + @warn "Passing noise scaling input into SDEProblem will be deprecated. New standard is to declare one (or several) paraemter as noise scaling parameters when the ReactionSystem is created. Please read https://docs.sciml.ai/Catalyst/stable/catalyst_applications/advanced_simulations/#Scaling-the-noise-magnitude-in-the-chemical-Langevin-equations." + end + isnothing(noise_scaling) && (noise_scaling = get_noise_scaling(flatrs)) + if noise_scaling isa AbstractArray (length(noise_scaling) != numreactions(flatrs)) && error("The number of elements in 'noise_scaling' must be equal " * @@ -1551,7 +1557,7 @@ function get_noise_scaling(rs::ReactionSystem) if isempty(ns_params) return nothing elseif length(ns_params) == 1 - return Num(ns_params[1]) + return ns_params[1] elseif length(ns_params) == length(reactions(rs)) return ns_params else @@ -1648,11 +1654,6 @@ function DiffEqBase.SDEProblem(rs::ReactionSystem, u0, tspan, check_length = false, remove_conserved = false, kwargs...) - if !isnothing(noise_scaling) - !isnothing(get_noise_scaling(rs)) && error("You have declared some parameters as noise scaling parameters, and also given a \"noise_scaling\" argument to SDEProblem. Please remove the \"noise_scaling\", as this way of scaling CLE noise is being depricated.") - @warn "Passing noise scaling input into SDEProblem will be deprecated. New standard is to declare one (or several) paraemter as noise scaling parameters when the ReactionSystem is created. Please read https://docs.sciml.ai/Catalyst/stable/catalyst_applications/advanced_simulations/#Scaling-the-noise-magnitude-in-the-chemical-Langevin-equations." - end - u0map = symmap_to_varmap(rs, u0) pmap = symmap_to_varmap(rs, p) sde_sys = convert(SDESystem, rs; noise_scaling, name, combinatoric_ratelaws, From 18bc7ea3ae03420647f3bc7421e3c8cac2fb1fef Mon Sep 17 00:00:00 2001 From: Torkel Date: Mon, 25 Sep 2023 21:25:58 -0400 Subject: [PATCH 16/41] test update. --- test/model_simulation/simulate_SDEs.jl | 20 +++++++++++++------- 1 file changed, 13 insertions(+), 7 deletions(-) diff --git a/test/model_simulation/simulate_SDEs.jl b/test/model_simulation/simulate_SDEs.jl index 8424d0d9fc..72a760207a 100644 --- a/test/model_simulation/simulate_SDEs.jl +++ b/test/model_simulation/simulate_SDEs.jl @@ -158,15 +158,21 @@ let end # Tests using default values for noise scaling. -let - noise_scaling_network = @reaction_network begin +# Tests when reaction system is created programmatically. +# Tests @noise_scaling_parameters macro. + let + @variables t + @species X1(t) X2(t) @noise_scaling_parameters η=0.0 - (k1, k2), X1 ↔ X2 + @parameters k1 k2 + r1 = Reaction(k1,[X1],[X2],[1],[1]) + r2 = Reaction(k2,[X2],[X1],[1],[1]) + @named noise_scaling_network = ReactionSystem([r1, r2], t, [X1, X2], [k1, k2, η]) + + u0 = [:X1 => 1100.0, :X2 => 3900.0] + p = [:k1 => 2.0, :k2 => 0.5, :η=>0.0] + @test SDEProblem(noise_scaling_network, u0, (0.0, 1000.0), p)[:η] == 0.0 end - u0 = [:X1 => 1100.0, :X2 => 3900.0] - p = [:k1 => 2.0, :k2 => 0.5, :η=>0.0] - @test SDEProblem(noise_scaling_network, u0, (0.0, 1000.0), p)[:η] == 0.0 -end # Complicated test with many combinations of options. let From 023f213d536075066f916593973f6bd6ab6ab4e8 Mon Sep 17 00:00:00 2001 From: Torkel Date: Tue, 26 Sep 2023 11:13:54 -0400 Subject: [PATCH 17/41] update --- docs/src/api/catalyst_api.md | 1 + .../advanced_simulations.md | 68 +++++++++++++------ src/Catalyst.jl | 2 +- src/networkapi.jl | 9 +++ test/model_simulation/simulate_SDEs.jl | 33 +++++---- 5 files changed, 79 insertions(+), 34 deletions(-) diff --git a/docs/src/api/catalyst_api.md b/docs/src/api/catalyst_api.md index 5085a60d3f..36a3b0fd48 100644 --- a/docs/src/api/catalyst_api.md +++ b/docs/src/api/catalyst_api.md @@ -155,6 +155,7 @@ accessor functions. ```@docs species nonspecies +noise_scaling_parameters reactionparams reactions numspecies diff --git a/docs/src/catalyst_applications/advanced_simulations.md b/docs/src/catalyst_applications/advanced_simulations.md index 0c1a2699d7..bf5c6c4dba 100644 --- a/docs/src/catalyst_applications/advanced_simulations.md +++ b/docs/src/catalyst_applications/advanced_simulations.md @@ -306,7 +306,7 @@ tspan = (0.0, 10.0) p_1 = [:k1 => 1.0, :k2 => 1.0] sprob_1 = SDEProblem(rn_1, u0, tspan, p_1) -sol_1 = solve(sprob_1) +sol_1 = solve(sprob_1, ImplicitEM()) plot(sol_1; idxs = :X1, ylimit = (0.0, 20.0)) ``` Here we can see that the `X` concentration fluctuates around a steady state of $X≈10.0$. @@ -319,46 +319,76 @@ Setting $η<1.0$ will reduce noise and $η>1.0$ will increase noise. The syntax for setting a noise scaling parameter `η` is ```@example ex3 rn_2 = @reaction_network begin - @parameters η + @noise_scaling_parameters η (k1,k2), X1 <--> X2 end +``` +The `@noise_scaling_parameter` option is equivalent to the `@parameters` option, but also designates subsequent parameter as a *noise scaling parameter*. *η* becomes a parameter of the system, and we can now modualte its value to scale simualtion noise: +```@example ex3 u0 = [:X1 => 10.0, :X2 => 10.0] tspan = (0.0, 10.0) p_2 = [:k1 => 1.0, :k2 => 1.0, :η => 0.1] -sprob_2 = SDEProblem(rn_2, u0, tspan, p_2; noise_scaling = (@parameters η)[1]) -``` -Here, we first need to add `η` as a parameter to the system using the -`@parameters η` option. Next, we pass the `noise_scaling = (@parameters η)[1]` -argument to the `SDEProblem`. We can now simulate our system and confirm that -noise is reduced: -```@example ex3 -sol_2 = solve(sprob_2) +sprob_2 = SDEProblem(rn_2, u0, tspan, p_2) +sol_2 = solve(sprob_2, ImplicitEM()) plot(sol_2; idxs = :X1, ylimit = (0.0, 20.0)) ``` -Finally, it is possible to set individual noise scaling parameters for each -reaction of the system. Our model has two reactions (`X1 --> X2` and `X2 --> -X1`) so we will use two noise scaling parameters, `η1` and `η2`. We use the -following syntax: +It is worth noting that in the CLE, nosie is tied to *reactions* (and not species, which is a common missperception). If only a single noise scaling parameter is given, it will scale the noise for all reaction. However, it is also possible to set several nosie scaling parameters, with each scaling the noise of a single reaction. Our model has two reactions (`X1 --> X2` and `X2 --> X1`) so we will use two noise scaling parameters (`η1` and `η2`): ```@example ex3 rn_3 = @reaction_network begin - @parameters η1 η2 + @noise_scaling_parameters η1 η2 (k1,k2), X1 <--> X2 end u0 = [:X1 => 10.0, :X2 => 10.0] tspan = (0.0, 10.0) p_3 = [:k1 => 1.0, :k2 => 1.0, :η1 => 0.1, :η2 => 1.0] -sprob_3 = SDEProblem(rn_3, u0, tspan, p_3; noise_scaling = @parameters η1 η2) +sprob_3 = SDEProblem(rn_3, u0, tspan, p_3) ``` -plotting the results, we see that we have less fluctuation than for the first -simulation, but more as compared to the second one (which is as expected): +Both the noise scaling parameters and the reaction are ordered (these orders can be seen by calling `reactions(rn_3)` and `noise_scaling_parameters(rn_3)`, respectively). The i'th noise scaling parameter scales the noise of the i'th reaction. Plotting the results, we see that we have less fluctuation than for the first simulation, but more as compared to the second one (which is as expected): ```@example ex3 -sol_3 = solve(sprob_3) +sol_3 = solve(sprob_3, ImplicitEM()) plot(sol_3; idxs = :X1, ylimit = (0.0, 20.0)) ``` +For systems with many reactions, the `η[1:n]` (where `n` is equal to the number of reactions) notation can be useful (this however, requires `@unpack`'ing the system's parameters): +```@example ex3 +rn_4 = @reaction_network begin + @noise_scaling_parameters η[1:6] + (p, d), 0 <--> X + (p, d), 0 <--> Y + (p, d), 0 <--> Z +end +@unpack p, d, η = rn_4 + +u0_4 = [:X => 10.0, :Y => 10.0, :Z => 10.0] +tspan = (0.0, 10.0) +p_4 = [p => 10.0, d => 1., η[1]=>0.1, η[2]=>0.1, η[3]=>1., η[4]=>1., η[5]=>1., η[6]=>1.] + +sprob_4 = SDEProblem(rn_4, u0_4, tspan, p_4) +sol_4 = solve(sprob_4, ImplicitEM()) +plot(sol_4; ylimit = (0.0, 20.0)) +``` +Here, we have reduced the noise of the reactions involving `X` only. + +Finally, it is possible to use the `@noise_scaling_parameter` option as a normal macro when creating reaction systems programmatically: +```@example ex3 +@variables t +@species X1(t) X2(t) +@noise_scaling_parameters η +@parameters k1 k2 +r1 = Reaction(k1,[X1],[X2],[1],[1]) +r2 = Reaction(k2,[X2],[X1],[1],[1]) +@named rn_5 = ReactionSystem([r1, r2], t, [X1, X2], [k1, k2, η]) +nothing # hide +``` +which is equivalent to `rn_2`. In this example, calling `@noise_scaling_parameters η` is equivalent to calling `parameters η` with the `noise_scaling_parameter` metadata: +```@example ex3 +@parameters η [noise_scaling_parameter=true] +nothing # hide +``` + ## Useful plotting options Catalyst, just like DifferentialEquations, uses the Plots package for all plotting. For a detailed description of differential equation plotting, see diff --git a/src/Catalyst.jl b/src/Catalyst.jl index 3f0573fd6a..d17b927f70 100644 --- a/src/Catalyst.jl +++ b/src/Catalyst.jl @@ -85,7 +85,7 @@ export mm, mmr, hill, hillr, hillar # functions to query network properties include("networkapi.jl") -export species, nonspecies, reactionparams, reactions, speciesmap, paramsmap +export species, nonspecies, noise_scaling_parameters, reactionparams, reactions, speciesmap, paramsmap export numspecies, numreactions, numreactionparams, setdefaults!, symmap_to_varmap export make_empty_network, addspecies!, addparam!, addreaction!, reactionparamsmap export dependants, dependents, substoichmat, prodstoichmat, netstoichmat diff --git a/src/networkapi.jl b/src/networkapi.jl index 1926fd34e2..2281763ea5 100644 --- a/src/networkapi.jl +++ b/src/networkapi.jl @@ -45,6 +45,15 @@ function nonspecies(network) unknowns(network)[(numspecies(network) + 1):end] end +""" + noise_scaling_parameters(network) + +Given a [`ReactionSystem`](@ref), return a vector of all noise scaling parameters defined in the system. +""" +function noise_scaling_parameters(network) + filter(is_noise_scaling_parameter, parameters(network)) +end + """ reactionparams(network) diff --git a/test/model_simulation/simulate_SDEs.jl b/test/model_simulation/simulate_SDEs.jl index 72a760207a..37c3400fa8 100644 --- a/test/model_simulation/simulate_SDEs.jl +++ b/test/model_simulation/simulate_SDEs.jl @@ -121,7 +121,7 @@ end ### Checks Simulations Don't Error ### -#Tries to create a large number of problem, ensuring there are no errors (cannot solve as solution likely to go into negatives). +# Tries to create a large number of problem, ensuring there are no errors (cannot solve as solution likely to go into negatives). let for reaction_network in reaction_networks_all for factor in [1e-2, 1e-1, 1e0, 1e1] @@ -160,21 +160,23 @@ end # Tests using default values for noise scaling. # Tests when reaction system is created programmatically. # Tests @noise_scaling_parameters macro. - let - @variables t - @species X1(t) X2(t) - @noise_scaling_parameters η=0.0 - @parameters k1 k2 - r1 = Reaction(k1,[X1],[X2],[1],[1]) - r2 = Reaction(k2,[X2],[X1],[1],[1]) - @named noise_scaling_network = ReactionSystem([r1, r2], t, [X1, X2], [k1, k2, η]) - - u0 = [:X1 => 1100.0, :X2 => 3900.0] - p = [:k1 => 2.0, :k2 => 0.5, :η=>0.0] - @test SDEProblem(noise_scaling_network, u0, (0.0, 1000.0), p)[:η] == 0.0 - end +let + η_stored = :η + @variables t + @species X1(t) X2(t) + ηs = @noise_scaling_parameters $(η_stored)=0.0 + @parameters k1 k2 + r1 = Reaction(k1,[X1],[X2],[1],[1]) + r2 = Reaction(k2,[X2],[X1],[1],[1]) + @named noise_scaling_network = ReactionSystem([r1, r2], t, [X1, X2], [k1, k2, ηs[1]]) + + u0 = [:X1 => 1100.0, :X2 => 3900.0] + p = [:k1 => 2.0, :k2 => 0.5, :η=>0.0] + @test SDEProblem(noise_scaling_network, u0, (0.0, 1000.0), p)[:η] == 0.0 +end # Complicated test with many combinations of options. +# Tests the noise_scaling_parameters getter. let noise_scaling_network = @reaction_network begin @parameters k1 par1 [description="Parameter par1"] par2 η1 [noise_scaling_parameter=true] @@ -190,6 +192,9 @@ let sprob = SDEProblem(noise_scaling_network, u0, (0.0, 1000.0), p) @test sprob[:η1] == sprob[:η2] == sprob[:η3] == sprob[:η4] == 0.0 + + @unpack η1, η2, η3, η4 = noise_scaling_network + isequal([η1, η2, η3, η4], noise_scaling_parameters(noise_scaling_network)) end From b8834f0a4aab6bbaa6eed82f4c622f725a19ae7b Mon Sep 17 00:00:00 2001 From: Torkel Date: Wed, 27 Sep 2023 12:14:08 -0400 Subject: [PATCH 18/41] Update docs/src/catalyst_applications/advanced_simulations.md Co-authored-by: Sam Isaacson --- docs/src/catalyst_applications/advanced_simulations.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/src/catalyst_applications/advanced_simulations.md b/docs/src/catalyst_applications/advanced_simulations.md index bf5c6c4dba..9f1bd1be57 100644 --- a/docs/src/catalyst_applications/advanced_simulations.md +++ b/docs/src/catalyst_applications/advanced_simulations.md @@ -323,7 +323,7 @@ rn_2 = @reaction_network begin (k1,k2), X1 <--> X2 end ``` -The `@noise_scaling_parameter` option is equivalent to the `@parameters` option, but also designates subsequent parameter as a *noise scaling parameter*. *η* becomes a parameter of the system, and we can now modualte its value to scale simualtion noise: +The `@noise_scaling_parameter` option creates one or more new parameters with additional metadata that allows Catalyst to know they represent a noise scaling. *η* becomes a parameter of the system, and we can now modulate its value to scale simulation noise: ```@example ex3 u0 = [:X1 => 10.0, :X2 => 10.0] tspan = (0.0, 10.0) From 6ee1ae02ddeea3f9dae0422fa79c6d3a04501f0e Mon Sep 17 00:00:00 2001 From: Torkel Date: Wed, 27 Sep 2023 12:14:17 -0400 Subject: [PATCH 19/41] Update docs/src/catalyst_applications/advanced_simulations.md Co-authored-by: Sam Isaacson --- docs/src/catalyst_applications/advanced_simulations.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/src/catalyst_applications/advanced_simulations.md b/docs/src/catalyst_applications/advanced_simulations.md index 9f1bd1be57..1cec092dca 100644 --- a/docs/src/catalyst_applications/advanced_simulations.md +++ b/docs/src/catalyst_applications/advanced_simulations.md @@ -334,7 +334,7 @@ sol_2 = solve(sprob_2, ImplicitEM()) plot(sol_2; idxs = :X1, ylimit = (0.0, 20.0)) ``` -It is worth noting that in the CLE, nosie is tied to *reactions* (and not species, which is a common missperception). If only a single noise scaling parameter is given, it will scale the noise for all reaction. However, it is also possible to set several nosie scaling parameters, with each scaling the noise of a single reaction. Our model has two reactions (`X1 --> X2` and `X2 --> X1`) so we will use two noise scaling parameters (`η1` and `η2`): +It is worth noting that in the CLE, nosie is tied to *reactions* (and not species, which is a common missperception). If only a single noise scaling parameter is given, it will scale the noise for all reactions. However, it is also possible to set several noise scaling parameters, with each scaling the noise of a single reaction. Our model has two reactions (`X1 --> X2` and `X2 --> X1`) so we will use two noise scaling parameters (`η1` and `η2`): ```@example ex3 rn_3 = @reaction_network begin @noise_scaling_parameters η1 η2 From 786fd3cefe47e24b1a0494ba274422e702ae3938 Mon Sep 17 00:00:00 2001 From: Torkel Date: Fri, 2 Feb 2024 17:12:47 -0500 Subject: [PATCH 20/41] save --- src/reaction_network.jl | 4 ---- 1 file changed, 4 deletions(-) diff --git a/src/reaction_network.jl b/src/reaction_network.jl index 86bf42565a..de428d5778 100644 --- a/src/reaction_network.jl +++ b/src/reaction_network.jl @@ -722,10 +722,6 @@ function extract_metadata(metadata_line::Expr) return metadata end -# macro, similar to @parameters, but all paraemters become noise scaling parameters. -macro noise_scaling_parameters(ex...) - vars = Symbolics._parse_vars(:parameters, Real, ex, toparam) - # Reads the observables options. Outputs an expression ofr creating the obervable variables, and a vector of observable equations. function read_observed_options(options, species_n_vars_declared, ivs_sorted) if haskey(options, :observables) From 61a84b763cb28eb172d2f3a8c6a7ed32021a10e4 Mon Sep 17 00:00:00 2001 From: Torkel Date: Sat, 3 Feb 2024 17:19:27 -0500 Subject: [PATCH 21/41] Update to use the new approach. --- docs/src/api/catalyst_api.md | 1 - .../advanced_simulations.md | 74 +++++---------- src/Catalyst.jl | 4 +- src/networkapi.jl | 9 -- src/reaction_network.jl | 61 ++++++++---- src/reactionsystem.jl | 95 ++++++------------- test/model_simulation/simulate_SDEs.jl | 73 ++++++++++---- .../reactionsystem.jl | 43 --------- 8 files changed, 151 insertions(+), 209 deletions(-) diff --git a/docs/src/api/catalyst_api.md b/docs/src/api/catalyst_api.md index 36a3b0fd48..5085a60d3f 100644 --- a/docs/src/api/catalyst_api.md +++ b/docs/src/api/catalyst_api.md @@ -155,7 +155,6 @@ accessor functions. ```@docs species nonspecies -noise_scaling_parameters reactionparams reactions numspecies diff --git a/docs/src/catalyst_applications/advanced_simulations.md b/docs/src/catalyst_applications/advanced_simulations.md index 1cec092dca..7fbf729d10 100644 --- a/docs/src/catalyst_applications/advanced_simulations.md +++ b/docs/src/catalyst_applications/advanced_simulations.md @@ -292,8 +292,9 @@ plot(sol) ## Scaling the noise magnitude in the chemical Langevin equations When using the CLE to generate SDEs from a CRN, it can sometimes be desirable to -scale the magnitude of the noise terms. This can be done by introducing a *noise -scaling parameter*. First, we simulate a simple two-state CRN model using the +scale the magnitude of the noise terms. Here, each reaction of the system generates a separate noise term in the CLE. If you require identical scaling for all reaction, the `@default_noise_scaling` option can be used. Else, you can supply a `noise_scaling` metadata for each individual reaction, describing how to scale the noise for that reaction. + +We begin with considering the first approach. First, we simulate a simple two-state CRN model using the CLE: ```@example ex3 using Catalyst, StochasticDiffEq, Plots @@ -309,85 +310,58 @@ sprob_1 = SDEProblem(rn_1, u0, tspan, p_1) sol_1 = solve(sprob_1, ImplicitEM()) plot(sol_1; idxs = :X1, ylimit = (0.0, 20.0)) ``` -Here we can see that the `X` concentration fluctuates around a steady state of $X≈10.0$. - -Next, we wish to introduce a noise scaling parameter ,`η`. This will scale the -noise magnitude so that for $η=0.0$ the system lacks noise (and its SDE -simulations are identical to its ODE simulations) and for $η=1.0$ noise is not -scaled (and SDE simulations are identical to as if no noise scaling was used). -Setting $η<1.0$ will reduce noise and $η>1.0$ will increase noise. The syntax -for setting a noise scaling parameter `η` is +Here we can see that the $X$ concentration fluctuates around a steady state of $X≈10.0$. + +Next, we wish increase the amount of noise with by a factor 2. To do so, we use the `@default_noise_scaling` option, to which we provide the desired scaling ```@example ex3 rn_2 = @reaction_network begin - @noise_scaling_parameters η + @default_noise_scaling 2 (k1,k2), X1 <--> X2 end ``` -The `@noise_scaling_parameter` option creates one or more new parameters with additional metadata that allows Catalyst to know they represent a noise scaling. *η* becomes a parameter of the system, and we can now modulate its value to scale simulation noise: +If we re-simualte the system we see that the amount of noise have increased: ```@example ex3 -u0 = [:X1 => 10.0, :X2 => 10.0] -tspan = (0.0, 10.0) -p_2 = [:k1 => 1.0, :k2 => 1.0, :η => 0.1] - -sprob_2 = SDEProblem(rn_2, u0, tspan, p_2) -sol_2 = solve(sprob_2, ImplicitEM()) -plot(sol_2; idxs = :X1, ylimit = (0.0, 20.0)) +sprob_1 = SDEProblem(rn_2, u0, tspan, p_1) +sol_1 = solve(sprob_1, ImplicitEM()) +plot(sol_1; idxs = :X1, ylimit = (0.0, 20.0)) ``` -It is worth noting that in the CLE, nosie is tied to *reactions* (and not species, which is a common missperception). If only a single noise scaling parameter is given, it will scale the noise for all reactions. However, it is also possible to set several noise scaling parameters, with each scaling the noise of a single reaction. Our model has two reactions (`X1 --> X2` and `X2 --> X1`) so we will use two noise scaling parameters (`η1` and `η2`): +It is possible to scale the amount of noise using any expression. A common use of this is to set a parameter which determines the amount of noise. Here we create a parameter $η$, and uses its value to scale the noise. ```@example ex3 +using Catalyst, StochasticDiffEq, Plots + rn_3 = @reaction_network begin - @noise_scaling_parameters η1 η2 + @parameters η + @default_noise_scaling η (k1,k2), X1 <--> X2 end u0 = [:X1 => 10.0, :X2 => 10.0] tspan = (0.0, 10.0) -p_3 = [:k1 => 1.0, :k2 => 1.0, :η1 => 0.1, :η2 => 1.0] +p_3 = [:k1 => 1.0, :k2 => 1.0, :η => 0.2] sprob_3 = SDEProblem(rn_3, u0, tspan, p_3) -``` -Both the noise scaling parameters and the reaction are ordered (these orders can be seen by calling `reactions(rn_3)` and `noise_scaling_parameters(rn_3)`, respectively). The i'th noise scaling parameter scales the noise of the i'th reaction. Plotting the results, we see that we have less fluctuation than for the first simulation, but more as compared to the second one (which is as expected): -```@example ex3 sol_3 = solve(sprob_3, ImplicitEM()) plot(sol_3; idxs = :X1, ylimit = (0.0, 20.0)) ``` +Here we saw how, by setting a small $η$ value, the amount of noise was reduced. + +It is possible to use a different noise scaling expression for each reaction. Here, each reaction's noise scaling expression is provided using the `noise_scaling` metadata. In the following example, we use this to turn the noise of for both reactions involving the species $Y$. -For systems with many reactions, the `η[1:n]` (where `n` is equal to the number of reactions) notation can be useful (this however, requires `@unpack`'ing the system's parameters): ```@example ex3 rn_4 = @reaction_network begin - @noise_scaling_parameters η[1:6] (p, d), 0 <--> X - (p, d), 0 <--> Y - (p, d), 0 <--> Z + (p, d), 0 <--> Y, ([noise_scaling=0.0, noise_scaling=0.0]) end -@unpack p, d, η = rn_4 -u0_4 = [:X => 10.0, :Y => 10.0, :Z => 10.0] +u0_4 = [:X => 10.0, :Y => 10.0] tspan = (0.0, 10.0) -p_4 = [p => 10.0, d => 1., η[1]=>0.1, η[2]=>0.1, η[3]=>1., η[4]=>1., η[5]=>1., η[6]=>1.] +p_4 = [p => 10.0, d => 1.] sprob_4 = SDEProblem(rn_4, u0_4, tspan, p_4) sol_4 = solve(sprob_4, ImplicitEM()) plot(sol_4; ylimit = (0.0, 20.0)) ``` -Here, we have reduced the noise of the reactions involving `X` only. - -Finally, it is possible to use the `@noise_scaling_parameter` option as a normal macro when creating reaction systems programmatically: -```@example ex3 -@variables t -@species X1(t) X2(t) -@noise_scaling_parameters η -@parameters k1 k2 -r1 = Reaction(k1,[X1],[X2],[1],[1]) -r2 = Reaction(k2,[X2],[X1],[1],[1]) -@named rn_5 = ReactionSystem([r1, r2], t, [X1, X2], [k1, k2, η]) -nothing # hide -``` -which is equivalent to `rn_2`. In this example, calling `@noise_scaling_parameters η` is equivalent to calling `parameters η` with the `noise_scaling_parameter` metadata: -```@example ex3 -@parameters η [noise_scaling_parameter=true] -nothing # hide -``` +Here, we not that there is n fluctuation in the value of $Y$. If the `@default_noise_scaling` option is used, its value is used for all reactions for which the `noise_scaling` metadata is unused. If `@default_noise_scaling` is not used, teh default noise scaling value is `1.0`. ## Useful plotting options Catalyst, just like DifferentialEquations, uses the Plots package for all diff --git a/src/Catalyst.jl b/src/Catalyst.jl index d17b927f70..b0f053bef7 100644 --- a/src/Catalyst.jl +++ b/src/Catalyst.jl @@ -77,7 +77,7 @@ export ODEProblem, const ExprValues = Union{Expr, Symbol, Float64, Int, Bool} include("expression_utils.jl") include("reaction_network.jl") -export @reaction_network, @add_reactions, @reaction, @species, @noise_scaling_parameters +export @reaction_network, @add_reactions, @reaction, @species # registers CRN specific functions using Symbolics.jl include("registered_functions.jl") @@ -85,7 +85,7 @@ export mm, mmr, hill, hillr, hillar # functions to query network properties include("networkapi.jl") -export species, nonspecies, noise_scaling_parameters, reactionparams, reactions, speciesmap, paramsmap +export species, nonspecies, eactionparams, reactions, speciesmap, paramsmap export numspecies, numreactions, numreactionparams, setdefaults!, symmap_to_varmap export make_empty_network, addspecies!, addparam!, addreaction!, reactionparamsmap export dependants, dependents, substoichmat, prodstoichmat, netstoichmat diff --git a/src/networkapi.jl b/src/networkapi.jl index 2281763ea5..1926fd34e2 100644 --- a/src/networkapi.jl +++ b/src/networkapi.jl @@ -45,15 +45,6 @@ function nonspecies(network) unknowns(network)[(numspecies(network) + 1):end] end -""" - noise_scaling_parameters(network) - -Given a [`ReactionSystem`](@ref), return a vector of all noise scaling parameters defined in the system. -""" -function noise_scaling_parameters(network) - filter(is_noise_scaling_parameter, parameters(network)) -end - """ reactionparams(network) diff --git a/src/reaction_network.jl b/src/reaction_network.jl index de428d5778..fa3d697868 100644 --- a/src/reaction_network.jl +++ b/src/reaction_network.jl @@ -123,7 +123,7 @@ const forbidden_variables_error = let end # Declares the keys used for various options. -const option_keys = (:species, :parameters, :variables, :ivs, :compounds, :observables) +const option_keys = (:species, :parameters, :variables, :ivs, :compounds, :observables, :default_noise_scaling) ### The @species macro, basically a copy of the @variables macro. ### macro species(ex...) @@ -361,17 +361,16 @@ function make_reaction_system(ex::Expr; name = :(gensym(:ReactionSystem))) option_lines)) # Reads options. + default_reaction_metadata = :([]) compound_expr, compound_species = read_compound_options(options) + check_default_noise_scaling!(default_reaction_metadata, options) # Parses reactions, species, and parameters. - reactions = get_reactions(reaction_lines) + reactions = get_reactions(reaction_lines; default_reaction_metadata) species_declared = [extract_syms(options, :species); compound_species] parameters_declared = extract_syms(options, :parameters) variables = extract_syms(options, :variables) - # Handles noise scaling parameters. - noise_scaling_pexpr = handle_noise_scaling_parameters!(parameters_declared, options) - # handle independent variables if haskey(options, :ivs) ivs = Tuple(extract_syms(options, :ivs)) @@ -394,7 +393,7 @@ function make_reaction_system(ex::Expr; name = :(gensym(:ReactionSystem))) declared_syms) species = vcat(species_declared, species_extracted) parameters = vcat(parameters_declared, parameters_extracted) - + # Checks for input errors. (sum(length.([reaction_lines, option_lines])) != length(ex.args)) && error("@reaction_network input contain $(length(ex.args) - sum(length.([reaction_lines,option_lines]))) malformed lines.") @@ -408,9 +407,7 @@ function make_reaction_system(ex::Expr; name = :(gensym(:ReactionSystem))) sexprs = get_sexpr(species_extracted, options; iv_symbols = ivs) vexprs = haskey(options, :variables) ? options[:variables] : :() pexprs = get_pexpr(parameters_extracted, options) - ns_ps, ns_pssym = scalarize_macro(haskey(options, :noise_scaling_parameters), noise_scaling_pexpr, "ns_ps") ps, pssym = scalarize_macro(!isempty(parameters), pexprs, "ps") - rs_p_syms = haskey(options, :noise_scaling_parameters) ? :(union($pssym, $ns_pssym)) : pssym vars, varssym = scalarize_macro(!isempty(variables), vexprs, "vars") sps, spssym = scalarize_macro(!isempty(species), sexprs, "specs") comps, compssym = scalarize_macro(!isempty(compound_species), compound_expr, "comps") @@ -418,10 +415,9 @@ function make_reaction_system(ex::Expr; name = :(gensym(:ReactionSystem))) for reaction in reactions push!(rxexprs.args, get_rxexprs(reaction)) end - # Returns the rephrased expression. + quote $ps - $ns_ps $ivexpr $vars $sps @@ -597,7 +593,8 @@ function get_reaction(line) end # Generates a vector containing a number of reaction structures, each containing the information about one reaction. -function get_reactions(exprs::Vector{Expr}, reactions = Vector{ReactionStruct}(undef, 0)) +function get_reactions(exprs::Vector{Expr}, reactions = Vector{ReactionStruct}(undef, 0); + default_reaction_metadata = []) for line in exprs # Reads core reaction information. arrow, rate, reaction, metadata = read_reaction_line(line) @@ -607,12 +604,16 @@ function get_reactions(exprs::Vector{Expr}, reactions = Vector{ReactionStruct}(u if typeof(rate) != Expr || rate.head != :tuple error("Error: Must provide a tuple of reaction rates when declaring a bi-directional reaction.") end - push_reactions!(reactions, reaction.args[2], reaction.args[3], rate.args[1], metadata.args[1], arrow) - push_reactions!(reactions, reaction.args[3], reaction.args[2], rate.args[2], metadata.args[2], arrow) + push_reactions!(reactions, reaction.args[2], reaction.args[3], rate.args[1], metadata.args[1], arrow; + default_reaction_metadata) + push_reactions!(reactions, reaction.args[3], reaction.args[2], rate.args[2], metadata.args[2], arrow; + default_reaction_metadata) elseif in(arrow, fwd_arrows) - push_reactions!(reactions, reaction.args[2], reaction.args[3], rate, metadata, arrow) + push_reactions!(reactions, reaction.args[2], reaction.args[3], rate, metadata, arrow; + default_reaction_metadata) elseif in(arrow, bwd_arrows) - push_reactions!(reactions, reaction.args[3], reaction.args[2], rate, metadata, arrow) + push_reactions!(reactions, reaction.args[3], reaction.args[2], rate, metadata, arrow; + default_reaction_metadata) else throw("Malformed reaction, invalid arrow type used in: $(MacroTools.striplines(line))") end @@ -645,7 +646,7 @@ end # Used to create multiple reactions from, for instance, `k, (X,Y) --> 0`. # Handles metadata, e.g. `1.0, Z --> 0, [noisescaling=η]`. function push_reactions!(reactions::Vector{ReactionStruct}, sub_line::ExprValues, prod_line::ExprValues, - rate::ExprValues, metadata::ExprValues, arrow::Symbol) + rate::ExprValues, metadata::ExprValues, arrow::Symbol; default_reaction_metadata::Expr) # The rates, substrates, products, and metadata may be in a tupple form (e.g. `k, (X,Y) --> 0`). # This finds the length of these tuples (or 1 if not in tuple forms). Errors if lengs inconsistent. lengs = (tup_leng(sub_line), tup_leng(prod_line), tup_leng(rate), tup_leng(metadata)) @@ -657,7 +658,8 @@ function push_reactions!(reactions::Vector{ReactionStruct}, sub_line::ExprValues for i in 1:maximum(lengs) # If the `only_use_rate` metadata was not provided, this has to be infered from the arrow used. metadata_i = get_tup_arg(metadata, i) - if all(arg.args[1] != :only_use_rate for arg in metadata_i.args) + merge_metadata!(metadata_i, default_reaction_metadata) + if all([arg.args[1] != :only_use_rate for arg in metadata_i.args]) push!(metadata_i.args, :(only_use_rate = $(in(arrow, pure_rate_arrows)))) end @@ -671,6 +673,16 @@ function push_reactions!(reactions::Vector{ReactionStruct}, sub_line::ExprValues end end +# Merge the metadata encoded by a reaction with the default metadata (for all reactions) given by the system. +# If a metadata value is present in both vectors, uses the non-default value. +function merge_metadata!(metadata, default_reaction_metadata) + for entry in default_reaction_metadata.args + if !in(entry.args[1], arg.args[1] for arg in metadata.args) + push!(metadata.args, entry) + end + end +end + function processmult(op, mult, stoich) if (mult isa Number) && (stoich isa Number) op(mult, stoich) @@ -722,6 +734,9 @@ function extract_metadata(metadata_line::Expr) return metadata end +### DSL Options Handling ### +# Most options handled in previous sections, when code re-organised, these should ideally be moved to the same place. + # Reads the observables options. Outputs an expression ofr creating the obervable variables, and a vector of observable equations. function read_observed_options(options, species_n_vars_declared, ivs_sorted) if haskey(options, :observables) @@ -801,6 +816,17 @@ function make_observed_eqs(observables_expr) end end +# Checks if the `@default_noise_scaling` option is used. If so, read its input and adds it as a +# default metadata value to the `default_reaction_metadata` vector. +function check_default_noise_scaling!(default_reaction_metadata, options) + if haskey(options, :default_noise_scaling) + if (length(options[:default_noise_scaling].args) != 3) # Becasue of how expressions are, 3 is used. + error("@default_noise_scaling should only have a single input, this appears not to be the case: \"$(options[:default_noise_scaling])\"") + end + push!(default_reaction_metadata.args, :(noise_scaling=$(options[:default_noise_scaling].args[3]))) + end +end + ### Functionality for expanding function call to custom and specific functions ### #Recursively traverses an expression and replaces special function call like "hill(...)" with the actual corresponding expression. @@ -814,7 +840,6 @@ function recursive_expand_functions!(expr::ExprValues) expr end - # ### Old functions (for deleting). # function get_rx_species_deletethis(rxs, ps) diff --git a/src/reactionsystem.jl b/src/reactionsystem.jl index aec89a9ef3..9d16a25501 100644 --- a/src/reactionsystem.jl +++ b/src/reactionsystem.jl @@ -54,17 +54,6 @@ end # variables drop_dynamics(s) = isconstant(s) || isbc(s) || (!isspecies(s)) -### Other Metadata Fields ### - -# Denotes that a parameter controls the scaling of noise in the CLE. -struct NoiseScalingParameter end -Symbolics.option_to_metadata_type(::Val{:noise_scaling_parameter}) = NoiseScalingParameter - -is_noise_scaling_parameter(s::Num) = is_noise_scaling_parameter(MT.value(s)) -function is_noise_scaling_parameter(s) - MT.getmetadata(s, NoiseScalingParameter, false) -end - """ $(TYPEDEF) @@ -906,7 +895,7 @@ Arguments: Example: ```julia -reaction = @reaction k, 0 --> X, [noise_scaling=0.0] +reaction = @reaction k, 0 --> X, [description="Production reaction"] get_metadata_dict(reaction) ``` """ @@ -925,8 +914,8 @@ Arguments: Example: ```julia -reaction = @reaction k, 0 --> X, [noise_scaling=0.0] -has_metadata(reaction, :noise_scaling) +reaction = @reaction k, 0 --> X, [description="Production reaction"] +has_metadata(reaction, :description) ``` """ function has_metadata(reaction::Reaction, md_key::Symbol) @@ -944,8 +933,8 @@ Arguments: Example: ```julia -reaction = @reaction k, 0 --> X, [noise_scaling=0.0] -get_metadata(reaction, :noise_scaling) +reaction = @reaction k, 0 --> X, [description="Production reaction"] +get_metadata(reaction, :description) ``` """ function get_metadata(reaction::Reaction, md_key::Symbol) @@ -956,6 +945,17 @@ function get_metadata(reaction::Reaction, md_key::Symbol) return metadata[findfirst(isequal(md_key, entry[1]) for entry in get_metadata_dict(reaction))][2] end +""" + get_noise_scaling(reaction::Reaction) + +Retrives the nosie scaling factor for the reaction's noise term (as generated by the reaction in an SDE). +If a specific noise scaling expression have been provided in the reactions "noise_scaling" metadata, +returns that. If not, returns `1.0`. +""" +function get_noise_scaling(reaction::Reaction) + has_metadata(reaction, :noise_scaling) ? get_metadata(reaction, :noise_scaling) : 1.0 +end + ######################## Conversion to ODEs/SDEs/jump, etc ############################## """ @@ -1057,7 +1057,7 @@ function assemble_drift(rs, ispcs; combinatoric_ratelaws = true, as_odes = true, end # this doesn't work with constraint equations currently -function assemble_diffusion(rs, sts, ispcs, noise_scaling; combinatoric_ratelaws = true, +function assemble_diffusion(rs, sts, ispcs; combinatoric_ratelaws = true, remove_conserved = false) # as BC species should ultimately get an equation, we include them in the noise matrix num_bcsts = count(isbc, get_unknowns(rs)) @@ -1075,7 +1075,7 @@ function assemble_diffusion(rs, sts, ispcs, noise_scaling; combinatoric_ratelaws for (j, rx) in enumerate(get_rxs(rs)) rlsqrt = sqrt(abs(oderatelaw(rx; combinatoric_ratelaw = combinatoric_ratelaws))) - (noise_scaling !== nothing) && (rlsqrt *= noise_scaling[j]) + rlsqrt *= get_noise_scaling(rx) remove_conserved && (rlsqrt = substitute(rlsqrt, depspec_submap)) for (spec, stoich) in rx.netstoich @@ -1485,14 +1485,6 @@ Notes: `combinatoric_ratelaws=false` for a ratelaw of `k*S^2`, i.e. the scaling factor is ignored. Defaults to the value given when the `ReactionSystem` was constructed (which itself defaults to true). -- `noise_scaling=nothing::Union{Vector{Num},Num,Nothing}` allows for linear scaling of the - noise in the chemical Langevin equations. If `nothing` is given, the default value as in - Gillespie 2000 is used. Alternatively, a `Num` can be given, this is added as a parameter - to the system (at the end of the parameter array). All noise terms are linearly scaled - with this value. The parameter may be one already declared in the `ReactionSystem`. - Finally, a `Vector{Num}` can be provided (the length must be equal to the number of - reactions). Here the noise for each reaction is scaled by the corresponding parameter in - the input vector. This input may contain repeat parameters. - `remove_conserved=false`, if set to `true` will calculate conservation laws of the underlying set of reactions (ignoring constraint equations), and then apply them to reduce the number of equations. @@ -1500,42 +1492,23 @@ Notes: differential equations. """ function Base.convert(::Type{<:SDESystem}, rs::ReactionSystem; - noise_scaling = nothing, name = nameof(rs), - combinatoric_ratelaws = get_combinatoric_ratelaws(rs), + name = nameof(rs), combinatoric_ratelaws = get_combinatoric_ratelaws(rs), include_zero_odes = true, checks = false, remove_conserved = false, default_u0 = Dict(), default_p = Dict(), defaults = _merge(Dict(default_u0), Dict(default_p)), + noise_scaling = nothing, # To be removed (used for deprication message only). kwargs...) + isnothing(noise_scaling) || error("Supplying a noise_scaling agument to SDESystems have been depricated. Please see the \"Scaling the noise magnitude in the chemical Langevin equations\" section of the \"Advanced Simulation Options\" page of the documentation for the new approach.") spatial_convert_err(rs::ReactionSystem, SDESystem) flatrs = Catalyst.flatten(rs) error_if_constraints(SDESystem, flatrs) - # Required until passing nosie into SDEProblem can be depricated. When properly deprecated, change the kwarg to noise_scaling = get_noise_scaling(rs). - if !isnothing(noise_scaling) - !isnothing(get_noise_scaling(rs)) && error("You have declared some parameters as noise scaling parameters, and also given a \"noise_scaling\" argument to SDEProblem. Please remove the \"noise_scaling\", as this way of scaling CLE noise is being depricated.") - @warn "Passing noise scaling input into SDEProblem will be deprecated. New standard is to declare one (or several) paraemter as noise scaling parameters when the ReactionSystem is created. Please read https://docs.sciml.ai/Catalyst/stable/catalyst_applications/advanced_simulations/#Scaling-the-noise-magnitude-in-the-chemical-Langevin-equations." - end - isnothing(noise_scaling) && (noise_scaling = get_noise_scaling(flatrs)) - - if noise_scaling isa AbstractArray - (length(noise_scaling) != numreactions(flatrs)) && - error("The number of elements in 'noise_scaling' must be equal " * - "to the number of reactions in the flattened reaction system.") - if !(noise_scaling isa Symbolics.Arr) - noise_scaling = value.(noise_scaling) - end - elseif !isnothing(noise_scaling) - noise_scaling = fill(value(noise_scaling), numreactions(flatrs)) - end - remove_conserved && conservationlaws(flatrs) ists, ispcs = get_indep_sts(flatrs, remove_conserved) eqs = assemble_drift(flatrs, ispcs; combinatoric_ratelaws, include_zero_odes, remove_conserved) - noiseeqs = assemble_diffusion(flatrs, ists, ispcs, noise_scaling; combinatoric_ratelaws, - remove_conserved) + noiseeqs = assemble_diffusion(flatrs, ists, ispcs; combinatoric_ratelaws, remove_conserved) eqs, sts, ps, obs, defs = addconstraints!(eqs, flatrs, ists, ispcs; remove_conserved) - ps = (noise_scaling === nothing) ? ps : vcat(ps, toparam(noise_scaling)) if any(isbc, get_unknowns(flatrs)) @info "Boundary condition species detected. As constraint equations are not currently supported when converting to SDESystems, the resulting system will be undetermined. Consider using constant species instead." @@ -1551,20 +1524,6 @@ function Base.convert(::Type{<:SDESystem}, rs::ReactionSystem; kwargs...) end -# Extracts any noise scaling parameters from a reaction system. -function get_noise_scaling(rs::ReactionSystem) - ns_params = filter(p -> is_noise_scaling_parameter(p), parameters(rs)) - if isempty(ns_params) - return nothing - elseif length(ns_params) == 1 - return ns_params[1] - elseif length(ns_params) == length(reactions(rs)) - return ns_params - else - error("The system have $(length(ns_params)) noise scaling parameters. This number should be equal to 0, 1, or the number of reactions ($(length(reactions(rs)))).") - end -end - """ ```julia Base.convert(::Type{<:JumpSystem},rs::ReactionSystem; combinatoric_ratelaws=true) @@ -1587,7 +1546,9 @@ function Base.convert(::Type{<:JumpSystem}, rs::ReactionSystem; name = nameof(rs combinatoric_ratelaws = get_combinatoric_ratelaws(rs), remove_conserved = nothing, checks = false, default_u0 = Dict(), default_p = Dict(), defaults = _merge(Dict(default_u0), Dict(default_p)), + noise_scaling = nothing, # To be removed (used for deprication message only). kwargs...) + isnothing(noise_scaling) || error("Supplying a noise_scaling agument to SDESystems have been depricated. Please see the \"Scaling the noise magnitude in the chemical Langevin equations\" section of the \"Advanced Simulation Options\" page of the documentation for the new approach.") spatial_convert_err(rs::ReactionSystem, JumpSystem) (remove_conserved !== nothing) && @@ -1648,15 +1609,13 @@ end # SDEProblem from AbstractReactionNetwork function DiffEqBase.SDEProblem(rs::ReactionSystem, u0, tspan, p = DiffEqBase.NullParameters(), args...; - noise_scaling = nothing, name = nameof(rs), - combinatoric_ratelaws = get_combinatoric_ratelaws(rs), - include_zero_odes = true, checks = false, - check_length = false, + name = nameof(rs), combinatoric_ratelaws = get_combinatoric_ratelaws(rs), + include_zero_odes = true, checks = false, check_length = false, remove_conserved = false, kwargs...) u0map = symmap_to_varmap(rs, u0) pmap = symmap_to_varmap(rs, p) - sde_sys = convert(SDESystem, rs; noise_scaling, name, combinatoric_ratelaws, + sde_sys = convert(SDESystem, rs; name, combinatoric_ratelaws, include_zero_odes, checks, remove_conserved) p_matrix = zeros(length(get_unknowns(sde_sys)), numreactions(rs)) return SDEProblem(sde_sys, u0map, tspan, pmap, args...; check_length, diff --git a/test/model_simulation/simulate_SDEs.jl b/test/model_simulation/simulate_SDEs.jl index 37c3400fa8..d4d25f27d1 100644 --- a/test/model_simulation/simulate_SDEs.jl +++ b/test/model_simulation/simulate_SDEs.jl @@ -137,8 +137,8 @@ end # Tests with multiple noise scaling parameters directly in the macro. let noise_scaling_network_1 = @reaction_network begin - @noise_scaling_parameters η1 η2 - (k1, k2), X1 ↔ X2 + @parameters η1 η2 + (k1, k2), X1 ↔ X2, ([noise_scaling=η1],[noise_scaling=η2]) end u0 = [:X1 => 1000.0, :X2 => 3000.0] sol_1_1 = solve(SDEProblem(noise_scaling_network_1, u0, (0.0, 1000.0), [:k1 => 2.0, :k2 => 0.66, :η1 => 2.0, :η2 => 2.0]), ImplicitEM(); saveat=1.0) @@ -148,7 +148,7 @@ let noise_scaling_network_2 = @reaction_network begin @noise_scaling_parameters η[1:2] - (k1, k2), X1 ↔ X2 + (k1, k2), X1 ↔ X2, ([noise_scaling=η[1]],[noise_scaling=η[2]]) end @unpack k1, k2, η = noise_scaling_network_2 sol_2_1 = solve(SDEProblem(noise_scaling_network_2, u0, (0.0, 1000.0), [k1 => 2.0, k2 => 0.66, η[1] => 2.0, η[2] => 2.0]), ImplicitEM(); saveat=1.0) @@ -164,37 +164,74 @@ let η_stored = :η @variables t @species X1(t) X2(t) - ηs = @noise_scaling_parameters $(η_stored)=0.0 - @parameters k1 k2 - r1 = Reaction(k1,[X1],[X2],[1],[1]) - r2 = Reaction(k2,[X2],[X1],[1],[1]) - @named noise_scaling_network = ReactionSystem([r1, r2], t, [X1, X2], [k1, k2, ηs[1]]) + p_syms = @parameters $(η_stored) k1 k2 + + r1 = Reaction(k1,[X1],[X2],[1],[1]; metadata = [:noise_scaling => η_stored]) + r2 = Reaction(k2,[X2],[X1],[1],[1]; metadata = [:noise_scaling => η_stored]) + @named noise_scaling_network = ReactionSystem([r1, r2], t, [X1, X2], [k1, k2, p_syms[1]]) u0 = [:X1 => 1100.0, :X2 => 3900.0] p = [:k1 => 2.0, :k2 => 0.5, :η=>0.0] - @test SDEProblem(noise_scaling_network, u0, (0.0, 1000.0), p)[:η] == 0.0 + @test_boken SDEProblem(noise_scaling_network, u0, (0.0, 1000.0), p).ps[:η] == 0.0 # Broken due to SII/MTK stuff. end # Complicated test with many combinations of options. # Tests the noise_scaling_parameters getter. let noise_scaling_network = @reaction_network begin - @parameters k1 par1 [description="Parameter par1"] par2 η1 [noise_scaling_parameter=true] - @noise_scaling_parameters η2=0.0 [description="Parameter η2"] η3=1.0 η4 - (p, d), 0 ↔ X1 - (k1, k2), X1 ↔ X2 + @parameters k1 par1 [description="Parameter par1"] par2 η1 η2=0.0 [description="Parameter η2"] η3=1.0 η4 + (p, d), 0 ↔ X1, ([noise_scaling=η1],[noise_scaling=η2]) + (k1, k2), X1 ↔ X2, ([noise_scaling=η3],[noise_scaling=η4]) end u0 = [:X1 => 500.0, :X2 => 500.0] p = [:p => 20.0, :d => 0.1, :η1 => 0.0, :η3 => 0.0, :η4 => 0.0, :k1 => 2.0, :k2 => 2.0, :par1 => 1000.0, :par2 => 1000.0] @test getdescription(parameters(noise_scaling_network)[2]) == "Parameter par1" - @test getdescription(parameters(noise_scaling_network)[8]) == "Parameter η2" - + @test getdescription(parameters(noise_scaling_network)[5]) == "Parameter η2" + sprob = SDEProblem(noise_scaling_network, u0, (0.0, 1000.0), p) - @test sprob[:η1] == sprob[:η2] == sprob[:η3] == sprob[:η4] == 0.0 + @test_broken sprob.ps[:η1] == sprob.ps[:η2] == sprob.ps[:η3] == sprob.ps[:η4] == 0.0 # Broken due to SII/MTK stuff. +end + +# Tests default_noise_scaling_option. +# Tests using sys. indexing for setting noise scaling parameter values. +# tests for default value used in noise scaling parameter. +let + noise_scaling_network = @reaction_network begin + @parameters η1 η2=0.1 + @default_noise_scaling η1 + (p,d), 0 <--> X1, ([noise_scaling=η2],[noise_scaling=η2]) + (p,d), 0 <--> X2, ([description="Y"],[description="Y"]) + (p,d), 0 <--> X3 + end + noise_scaling_network = complete(noise_scaling_network) + + u0 = [:X1 => 1000.0, :X2 => 1000.0, :X3 => 1000.0] + ps = [noise_scaling_network.p => 1000.0, noise_scaling_network.d => 1.0, noise_scaling_network.η1 => 1.0] + sol = solve(SDEProblem(noise_scaling_network, u0, (0.0, 1000.0), ps), ImplicitEM(); saveat=1.0) + @test var(sol[:X1]) < var(sol[:X2]) + @test var(sol[:X1]) < var(sol[:X3]) +end + +# Tests using complicated noise scaling expressions +let + noise_scaling_network = @reaction_network begin + @parameters η1 η2 η3 η4 + @species N1(t) N2(t)=0.5 + @variables N3(t) + @default_noise_scaling η1 + N1 + 5.0 + p, 0 --> X1 + p, 0 --> X2, [noise_scaling=0.33η2^1.2 + N2] + p, 0 --> X3, [noise_scaling=N3*η3] + p, 0 --> X4, [noise_scaling=exp(-η4) - 0.008] + p, 0 --> X5, [noise_scaling=0.0] + d, (X1, X2, X3, X4, X5) --> 0, ([], [noise_scaling=0.33η2^1.2 + N2], [noise_scaling=N3*η3], [noise_scaling=exp(-η4) - 0.008], [noise_scaling=0.0]) + end - @unpack η1, η2, η3, η4 = noise_scaling_network - isequal([η1, η2, η3, η4], noise_scaling_parameters(noise_scaling_network)) + u0 = [:X1 => 1000.0, :X2 => 1000.0, :X3 => 1000.0, :X4 => 1000.0, :X5 => 1000.0, :N1 => 3.0, :N3 => 0.33] + ps = [:p => 1000.0, :d => 1.0, :η1 => 1.0, :η2 => 1.4, :η3 => 0.33, :η4 => 4.0] + sol = solve(SDEProblem(noise_scaling_network, u0, (0.0, 1000.0), ps), ImplicitEM(); saveat=1.0, adaptive=false, dt=0.1) + @test var(sol[:X1]) > var(sol[:X2]) > var(sol[:X3]) > var(sol[:X4]) > var(sol[:X5]) end diff --git a/test/reactionsystem_structure/reactionsystem.jl b/test/reactionsystem_structure/reactionsystem.jl index 10f1129033..a977143ce7 100644 --- a/test/reactionsystem_structure/reactionsystem.jl +++ b/test/reactionsystem_structure/reactionsystem.jl @@ -162,49 +162,6 @@ let @test norm(du - dunl) < 100 * eps() end -# Tests the noise_scaling argument. -let - p = rand(rng, length(k) + 1) - u = rand(rng, length(k)) - t = 0.0 - G = p[21] * sdenoise(u, p, t) - @variables η - sdesys_noise_scaling = convert(SDESystem, rs; noise_scaling = η) - sf = SDEFunction{false}(sdesys_noise_scaling, unknowns(rs), - parameters(sdesys_noise_scaling)) - G2 = sf.g(u, p, t) - @test norm(G - G2) < 100 * eps() -end - -# Tests the noise_scaling vector argument. -let - p = rand(rng, length(k) + 3) - u = rand(rng, length(k)) - t = 0.0 - G = vcat(fill(p[21], 8), fill(p[22], 3), fill(p[23], 9))' .* sdenoise(u, p, t) - @variables η[1:3] - sdesys_noise_scaling = convert(SDESystem, rs; - noise_scaling = vcat(fill(η[1], 8), fill(η[2], 3), - fill(η[3], 9))) - sf = SDEFunction{false}(sdesys_noise_scaling, unknowns(rs), - parameters(sdesys_noise_scaling)) - G2 = sf.g(u, p, t) - @test norm(G - G2) < 100 * eps() -end - -# Tests using previous parameter for noise scaling -let - p = rand(rng, length(k)) - u = rand(rng, length(k)) - t = 0.0 - G = [p p p p]' .* sdenoise(u, p, t) - sdesys_noise_scaling = convert(SDESystem, rs; noise_scaling = k) - sf = SDEFunction{false}(sdesys_noise_scaling, unknowns(rs), - parameters(sdesys_noise_scaling)) - G2 = sf.g(u, p, t) - @test norm(G - G2) < 100 * eps() -end - # Test with JumpSystem. let p = rand(rng, length(k)) From ebb21e42a35fc9080587a6fc8e1eee2ca72f662d Mon Sep 17 00:00:00 2001 From: Torkel Date: Sat, 3 Feb 2024 17:38:50 -0500 Subject: [PATCH 22/41] up --- src/reaction_network.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/reaction_network.jl b/src/reaction_network.jl index fa3d697868..048bab3d05 100644 --- a/src/reaction_network.jl +++ b/src/reaction_network.jl @@ -594,7 +594,7 @@ end # Generates a vector containing a number of reaction structures, each containing the information about one reaction. function get_reactions(exprs::Vector{Expr}, reactions = Vector{ReactionStruct}(undef, 0); - default_reaction_metadata = []) + default_reaction_metadata = :([])) for line in exprs # Reads core reaction information. arrow, rate, reaction, metadata = read_reaction_line(line) From 2322fb8aecc7e05412247be28b4f0c6d191ac4fa Mon Sep 17 00:00:00 2001 From: Torkel Date: Sun, 4 Feb 2024 16:03:51 -0500 Subject: [PATCH 23/41] up --- src/Catalyst.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Catalyst.jl b/src/Catalyst.jl index b0f053bef7..6bd26b35c9 100644 --- a/src/Catalyst.jl +++ b/src/Catalyst.jl @@ -85,7 +85,7 @@ export mm, mmr, hill, hillr, hillar # functions to query network properties include("networkapi.jl") -export species, nonspecies, eactionparams, reactions, speciesmap, paramsmap +export species, nonspecies, reactionparams, reactions, speciesmap, paramsmap export numspecies, numreactions, numreactionparams, setdefaults!, symmap_to_varmap export make_empty_network, addspecies!, addparam!, addreaction!, reactionparamsmap export dependants, dependents, substoichmat, prodstoichmat, netstoichmat From 35afb4aedb7c9bfbf412606f7123aefe70a6db56 Mon Sep 17 00:00:00 2001 From: Torkel Date: Sun, 4 Feb 2024 16:30:09 -0500 Subject: [PATCH 24/41] up --- test/model_simulation/simulate_SDEs.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/model_simulation/simulate_SDEs.jl b/test/model_simulation/simulate_SDEs.jl index d4d25f27d1..6a15ff48bb 100644 --- a/test/model_simulation/simulate_SDEs.jl +++ b/test/model_simulation/simulate_SDEs.jl @@ -147,7 +147,7 @@ let @test var(sol_1_1[1,:]) > var(sol_1_2[1,:]) > var(sol_1_3[1,:]) noise_scaling_network_2 = @reaction_network begin - @noise_scaling_parameters η[1:2] + @parameters η[1:2] (k1, k2), X1 ↔ X2, ([noise_scaling=η[1]],[noise_scaling=η[2]]) end @unpack k1, k2, η = noise_scaling_network_2 From 950e09866a909e6af80f92dac7a05d24613e23af Mon Sep 17 00:00:00 2001 From: Torkel Date: Sun, 4 Feb 2024 16:46:21 -0500 Subject: [PATCH 25/41] test tweak --- test/model_simulation/simulate_SDEs.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/model_simulation/simulate_SDEs.jl b/test/model_simulation/simulate_SDEs.jl index 6a15ff48bb..ec529bebb3 100644 --- a/test/model_simulation/simulate_SDEs.jl +++ b/test/model_simulation/simulate_SDEs.jl @@ -172,7 +172,7 @@ let u0 = [:X1 => 1100.0, :X2 => 3900.0] p = [:k1 => 2.0, :k2 => 0.5, :η=>0.0] - @test_boken SDEProblem(noise_scaling_network, u0, (0.0, 1000.0), p).ps[:η] == 0.0 # Broken due to SII/MTK stuff. + @test_broken SDEProblem(noise_scaling_network, u0, (0.0, 1000.0), p).ps[:η] == 0.0 # Broken due to SII/MTK stuff. end # Complicated test with many combinations of options. From 14b7e734e65ea77669dbcc3c6613b71b15d3b991 Mon Sep 17 00:00:00 2001 From: Torkel Date: Wed, 7 Feb 2024 10:46:35 -0500 Subject: [PATCH 26/41] broken test is no longer broken due to SII/MTK stuff --- test/model_simulation/simulate_SDEs.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/model_simulation/simulate_SDEs.jl b/test/model_simulation/simulate_SDEs.jl index ec529bebb3..a6018fbd4d 100644 --- a/test/model_simulation/simulate_SDEs.jl +++ b/test/model_simulation/simulate_SDEs.jl @@ -190,7 +190,7 @@ let @test getdescription(parameters(noise_scaling_network)[5]) == "Parameter η2" sprob = SDEProblem(noise_scaling_network, u0, (0.0, 1000.0), p) - @test_broken sprob.ps[:η1] == sprob.ps[:η2] == sprob.ps[:η3] == sprob.ps[:η4] == 0.0 # Broken due to SII/MTK stuff. + @test sprob.ps[:η1] == sprob.ps[:η2] == sprob.ps[:η3] == sprob.ps[:η4] == 0.0 end # Tests default_noise_scaling_option. From 99507e695fbba1787017b35121a6390a0d254740 Mon Sep 17 00:00:00 2001 From: Torkel Date: Mon, 18 Mar 2024 15:14:54 -0400 Subject: [PATCH 27/41] up --- HISTORY.md | 1 + src/reactionsystem.jl | 14 +------------- 2 files changed, 2 insertions(+), 13 deletions(-) diff --git a/HISTORY.md b/HISTORY.md index 4577d231a8..fdbbfea5a4 100644 --- a/HISTORY.md +++ b/HISTORY.md @@ -24,6 +24,7 @@ rn = @reaction_network begin end getnoisescaling(rn) ``` +- `SDEProblem` no longer takes the `noise_scaling` argument (she above for new approach to handle noise scaling). - Changed fields of internal `Reaction` structure. `ReactionSystems`s saved using `serialize` on previous Catalyst versions cannot be loaded using this (or later) versions. - Simulation of spatial ODEs now supported. For full details, please see https://github.com/SciML/Catalyst.jl/pull/644 and upcoming documentation. Note that these methods are currently considered alpha, with the interface and approach changing even in non-breaking Catalyst releases. - LatticeReactionSystem structure represents a spatial reaction network: diff --git a/src/reactionsystem.jl b/src/reactionsystem.jl index 9d16a25501..49f6982ef5 100644 --- a/src/reactionsystem.jl +++ b/src/reactionsystem.jl @@ -945,16 +945,6 @@ function get_metadata(reaction::Reaction, md_key::Symbol) return metadata[findfirst(isequal(md_key, entry[1]) for entry in get_metadata_dict(reaction))][2] end -""" - get_noise_scaling(reaction::Reaction) - -Retrives the nosie scaling factor for the reaction's noise term (as generated by the reaction in an SDE). -If a specific noise scaling expression have been provided in the reactions "noise_scaling" metadata, -returns that. If not, returns `1.0`. -""" -function get_noise_scaling(reaction::Reaction) - has_metadata(reaction, :noise_scaling) ? get_metadata(reaction, :noise_scaling) : 1.0 -end ######################## Conversion to ODEs/SDEs/jump, etc ############################## @@ -1075,7 +1065,7 @@ function assemble_diffusion(rs, sts, ispcs; combinatoric_ratelaws = true, for (j, rx) in enumerate(get_rxs(rs)) rlsqrt = sqrt(abs(oderatelaw(rx; combinatoric_ratelaw = combinatoric_ratelaws))) - rlsqrt *= get_noise_scaling(rx) + rlsqrt *= getnoisescaling(rx) remove_conserved && (rlsqrt = substitute(rlsqrt, depspec_submap)) for (spec, stoich) in rx.netstoich @@ -1497,7 +1487,6 @@ function Base.convert(::Type{<:SDESystem}, rs::ReactionSystem; default_u0 = Dict(), default_p = Dict(), defaults = _merge(Dict(default_u0), Dict(default_p)), noise_scaling = nothing, # To be removed (used for deprication message only). kwargs...) - isnothing(noise_scaling) || error("Supplying a noise_scaling agument to SDESystems have been depricated. Please see the \"Scaling the noise magnitude in the chemical Langevin equations\" section of the \"Advanced Simulation Options\" page of the documentation for the new approach.") spatial_convert_err(rs::ReactionSystem, SDESystem) flatrs = Catalyst.flatten(rs) @@ -1548,7 +1537,6 @@ function Base.convert(::Type{<:JumpSystem}, rs::ReactionSystem; name = nameof(rs default_u0 = Dict(), default_p = Dict(), defaults = _merge(Dict(default_u0), Dict(default_p)), noise_scaling = nothing, # To be removed (used for deprication message only). kwargs...) - isnothing(noise_scaling) || error("Supplying a noise_scaling agument to SDESystems have been depricated. Please see the \"Scaling the noise magnitude in the chemical Langevin equations\" section of the \"Advanced Simulation Options\" page of the documentation for the new approach.") spatial_convert_err(rs::ReactionSystem, JumpSystem) (remove_conserved !== nothing) && From 7dbb02cf811e5dee7481bc5f8a9e05fe0fab97e4 Mon Sep 17 00:00:00 2001 From: Torkel Date: Mon, 18 Mar 2024 17:46:08 -0400 Subject: [PATCH 28/41] up --- HISTORY.md | 2 +- src/reaction_network.jl | 2 +- src/reactionsystem.jl | 68 ++++++++++++------- test/dsl/dsl_basics.jl | 28 ++++---- test/reactionsystem_structure/reactions.jl | 44 ++++++------ .../reactionsystem.jl | 2 +- 6 files changed, 82 insertions(+), 64 deletions(-) diff --git a/HISTORY.md b/HISTORY.md index fdbbfea5a4..5921f0193c 100644 --- a/HISTORY.md +++ b/HISTORY.md @@ -22,7 +22,7 @@ rn = @reaction_network begin @parameters η k, 2X --> X2, [noise_scaling=η] end -getnoisescaling(rn) +get_noise_scaling(rn) ``` - `SDEProblem` no longer takes the `noise_scaling` argument (she above for new approach to handle noise scaling). - Changed fields of internal `Reaction` structure. `ReactionSystems`s saved using `serialize` on previous Catalyst versions cannot be loaded using this (or later) versions. diff --git a/src/reaction_network.jl b/src/reaction_network.jl index 048bab3d05..a97e47eb01 100644 --- a/src/reaction_network.jl +++ b/src/reaction_network.jl @@ -644,7 +644,7 @@ end # Takes a reaction line and creates reaction(s) from it and pushes those to the reaction array. # Used to create multiple reactions from, for instance, `k, (X,Y) --> 0`. -# Handles metadata, e.g. `1.0, Z --> 0, [noisescaling=η]`. +# Handles metadata, e.g. `1.0, Z --> 0, [noise_scaling=η]`. function push_reactions!(reactions::Vector{ReactionStruct}, sub_line::ExprValues, prod_line::ExprValues, rate::ExprValues, metadata::ExprValues, arrow::Symbol; default_reaction_metadata::Expr) # The rates, substrates, products, and metadata may be in a tupple form (e.g. `k, (X,Y) --> 0`). diff --git a/src/reactionsystem.jl b/src/reactionsystem.jl index 49f6982ef5..1f395602a6 100644 --- a/src/reactionsystem.jl +++ b/src/reactionsystem.jl @@ -864,10 +864,27 @@ end ######################## Other accessors ############################## """ - getnoisescaling(reaction::Reaction) +has_noise_scaling(reaction::Reaction) -Returns the noise scaling associated with a specific reaction. If the `:noise_scaling` metadata has -set, returns that. Else, returns `1.0`. +Checks whether a specific reaction has the metadata field `noise_scaing`. If so, returns `true`, else +returns `false`. + +Arguments: +- `reaction`: The reaction for which we wish to check. + +Example: +```julia +reaction = @reaction k, 0 --> X, [noise_scaling=0.0] +has_noise_scaling(reaction) +""" +function has_noise_scaling(reaction::Reaction) + return hasmetadata(reaction, :noise_scaling) +end + +""" +get_noise_scaling(reaction::Reaction) + +Returns the noise_scaling metadata from a specific reaction. Arguments: - `reaction`: The reaction for which we wish to retrive all metadata. @@ -875,18 +892,21 @@ Arguments: Example: ```julia reaction = @reaction k, 0 --> X, [noise_scaling=0.0] -getnoisescaling(reaction) +get_noise_scaling(reaction) """ -function getnoisescaling(reaction::Reaction) - has_metadata(reaction, :noise_scaling) && (return get_metadata(reaction, :noise_scaling)) - return 1.0 +function get_noise_scaling(reaction::Reaction) + if has_noise_scaling(reaction) + return getmetadata(reaction, :noise_scaling) + else + error("Attempts to access noise_scaling metadata field for a reaction which does not have a value assigned for this metadata.") + end end -# These are currently considered internal, but can be used by public accessor functions like getnoisescaling. +# These are currently considered internal, but can be used by public accessor functions like get_noise_scaling. """ - get_metadata_dict(reaction::Reaction) + getmetadata_dict(reaction::Reaction) Retrives the `ImmutableDict` containing all of the metadata associated with a specific reaction. @@ -896,15 +916,15 @@ Arguments: Example: ```julia reaction = @reaction k, 0 --> X, [description="Production reaction"] -get_metadata_dict(reaction) +getmetadata_dict(reaction) ``` """ -function get_metadata_dict(reaction::Reaction) +function getmetadata_dict(reaction::Reaction) return reaction.metadata end """ - has_metadata(reaction::Reaction, md_key::Symbol) + hasmetadata(reaction::Reaction, md_key::Symbol) Checks if a `Reaction` have a certain metadata field. If it does, returns `true` (else returns `false`). @@ -915,15 +935,15 @@ Arguments: Example: ```julia reaction = @reaction k, 0 --> X, [description="Production reaction"] -has_metadata(reaction, :description) +hasmetadata(reaction, :description) ``` """ -function has_metadata(reaction::Reaction, md_key::Symbol) - return any(isequal(md_key, entry[1]) for entry in get_metadata_dict(reaction)) +function hasmetadata(reaction::Reaction, md_key::Symbol) + return any(isequal(md_key, entry[1]) for entry in getmetadata_dict(reaction)) end """ - get_metadata(reaction::Reaction, md_key::Symbol) +getmetadata(reaction::Reaction, md_key::Symbol) Retrives a certain metadata value from a `Reaction`. If the metadata does not exists, throws an error. @@ -934,15 +954,15 @@ Arguments: Example: ```julia reaction = @reaction k, 0 --> X, [description="Production reaction"] -get_metadata(reaction, :description) +getmetadata(reaction, :description) ``` """ -function get_metadata(reaction::Reaction, md_key::Symbol) - if !has_metadata(reaction, md_key) - error("The reaction does not have a metadata field $md_key. It does have the following metadata fields: $(keys(get_metadata_dict(reaction))).") +function getmetadata(reaction::Reaction, md_key::Symbol) + if !hasmetadata(reaction, md_key) + error("The reaction does not have a metadata field $md_key. It does have the following metadata fields: $(keys(getmetadata_dict(reaction))).") end - metadata = get_metadata_dict(reaction) - return metadata[findfirst(isequal(md_key, entry[1]) for entry in get_metadata_dict(reaction))][2] + metadata = getmetadata_dict(reaction) + return metadata[findfirst(isequal(md_key, entry[1]) for entry in getmetadata_dict(reaction))][2] end @@ -1065,7 +1085,7 @@ function assemble_diffusion(rs, sts, ispcs; combinatoric_ratelaws = true, for (j, rx) in enumerate(get_rxs(rs)) rlsqrt = sqrt(abs(oderatelaw(rx; combinatoric_ratelaw = combinatoric_ratelaws))) - rlsqrt *= getnoisescaling(rx) + has_noise_scaling(rx) && (rlsqrt *= get_noise_scaling(rx)) remove_conserved && (rlsqrt = substitute(rlsqrt, depspec_submap)) for (spec, stoich) in rx.netstoich @@ -1485,7 +1505,6 @@ function Base.convert(::Type{<:SDESystem}, rs::ReactionSystem; name = nameof(rs), combinatoric_ratelaws = get_combinatoric_ratelaws(rs), include_zero_odes = true, checks = false, remove_conserved = false, default_u0 = Dict(), default_p = Dict(), defaults = _merge(Dict(default_u0), Dict(default_p)), - noise_scaling = nothing, # To be removed (used for deprication message only). kwargs...) spatial_convert_err(rs::ReactionSystem, SDESystem) @@ -1535,7 +1554,6 @@ function Base.convert(::Type{<:JumpSystem}, rs::ReactionSystem; name = nameof(rs combinatoric_ratelaws = get_combinatoric_ratelaws(rs), remove_conserved = nothing, checks = false, default_u0 = Dict(), default_p = Dict(), defaults = _merge(Dict(default_u0), Dict(default_p)), - noise_scaling = nothing, # To be removed (used for deprication message only). kwargs...) spatial_convert_err(rs::ReactionSystem, JumpSystem) diff --git a/test/dsl/dsl_basics.jl b/test/dsl/dsl_basics.jl index 053bfcddf9..d970242062 100644 --- a/test/dsl/dsl_basics.jl +++ b/test/dsl/dsl_basics.jl @@ -169,20 +169,20 @@ let # Checks DSL reactions are correct. rxs = reactions(rs) @test isequal([r1, r2, r3], rxs) - @test isequal(Catalyst.get_metadata_dict(r1), Catalyst.get_metadata_dict(rxs[1])) - @test isequal(Catalyst.get_metadata_dict(r2), Catalyst.get_metadata_dict(rxs[2])) - @test isequal(Catalyst.get_metadata_dict(r3), Catalyst.get_metadata_dict(rxs[3])) + @test isequal(Catalyst.getmetadata_dict(r1), Catalyst.getmetadata_dict(rxs[1])) + @test isequal(Catalyst.getmetadata_dict(r2), Catalyst.getmetadata_dict(rxs[2])) + @test isequal(Catalyst.getmetadata_dict(r3), Catalyst.getmetadata_dict(rxs[3])) # Checks that accessor functions works on the DSL. - @test Catalyst.has_metadata(rxs[1], :noise_scaling) - @test !Catalyst.has_metadata(rxs[1], :md_1) - @test !Catalyst.has_metadata(rxs[2], :noise_scaling) - @test Catalyst.has_metadata(rxs[2], :md_1) - @test !Catalyst.has_metadata(rxs[3], :noise_scaling) - @test !Catalyst.has_metadata(rxs[3], :md_1) + @test Catalyst.hasmetadata(rxs[1], :noise_scaling) + @test !Catalyst.hasmetadata(rxs[1], :md_1) + @test !Catalyst.hasmetadata(rxs[2], :noise_scaling) + @test Catalyst.hasmetadata(rxs[2], :md_1) + @test !Catalyst.hasmetadata(rxs[3], :noise_scaling) + @test !Catalyst.hasmetadata(rxs[3], :md_1) - @test isequal(Catalyst.get_metadata(rxs[1], :noise_scaling), η) - @test isequal(Catalyst.get_metadata(rxs[2], :md_1), 1.0) + @test isequal(Catalyst.getmetadata(rxs[1], :noise_scaling), η) + @test isequal(Catalyst.getmetadata(rxs[2], :md_1), 1.0) # Test that metadata works for @reaction macro. rx1 = @reaction k, 2X --> X2, [noise_scaling=$η] @@ -190,9 +190,9 @@ let rx3 = @reaction k, 2X --> X2 @test isequal([rx1, rx2, rx3], rxs) - @test isequal(Catalyst.get_metadata_dict(rx1), Catalyst.get_metadata_dict(rxs[1])) - @test isequal(Catalyst.get_metadata_dict(rx2), Catalyst.get_metadata_dict(rxs[2])) - @test isequal(Catalyst.get_metadata_dict(rx3), Catalyst.get_metadata_dict(rxs[3])) + @test isequal(Catalyst.getmetadata_dict(rx1), Catalyst.getmetadata_dict(rxs[1])) + @test isequal(Catalyst.getmetadata_dict(rx2), Catalyst.getmetadata_dict(rxs[2])) + @test isequal(Catalyst.getmetadata_dict(rx3), Catalyst.getmetadata_dict(rxs[3])) end # Checks that repeated metadata throws errors. diff --git a/test/reactionsystem_structure/reactions.jl b/test/reactionsystem_structure/reactions.jl index d51623156a..0746e3dfbe 100644 --- a/test/reactionsystem_structure/reactions.jl +++ b/test/reactionsystem_structure/reactions.jl @@ -16,10 +16,10 @@ let metadata = [:noise_scaling => 0.0] r = Reaction(k, [X], [X2], [2], [1]; metadata=metadata) - @test Catalyst.get_metadata_dict(r) == [:noise_scaling => 0.0] - @test Catalyst.has_metadata(r, :noise_scaling) - @test !Catalyst.has_metadata(r, :nonexisting_metadata) - @test Catalyst.get_metadata(r, :noise_scaling) == 0.0 + @test Catalyst.getmetadata_dict(r) == [:noise_scaling => 0.0] + @test Catalyst.hasmetadata(r, :noise_scaling) + @test !Catalyst.hasmetadata(r, :nonexisting_metadata) + @test Catalyst.getmetadata(r, :noise_scaling) == 0.0 metadata_repeated = [:noise_scaling => 0.0, :noise_scaling => 1.0, :metadata_entry => "unused"] @test_throws Exception Reaction(k, [X], [X2], [2], [1]; metadata=metadata_repeated) @@ -36,8 +36,8 @@ let r2 = Reaction(k, [X], [X2], [2], [1]; metadata=metadata) @test isequal(r1, r2) - @test Catalyst.get_metadata_dict(r1) == Pair{Symbol,Any}[] - @test !Catalyst.has_metadata(r1, :md) + @test Catalyst.getmetadata_dict(r1) == Pair{Symbol,Any}[] + @test !Catalyst.hasmetadata(r1, :md) end # Tests creation. @@ -57,21 +57,21 @@ let push!(metadata, :md_6 => (0.1, 2.0)) r = Reaction(k, [X], [X2], [2], [1]; metadata=metadata) - @test Catalyst.get_metadata_dict(r) isa Vector{Pair{Symbol,Any}} - @test Catalyst.has_metadata(r, :md_1) - @test Catalyst.has_metadata(r, :md_2) - @test Catalyst.has_metadata(r, :md_3) - @test Catalyst.has_metadata(r, :md_4) - @test Catalyst.has_metadata(r, :md_5) - @test Catalyst.has_metadata(r, :md_6) - @test !Catalyst.has_metadata(r, :md_8) + @test Catalyst.getmetadata_dict(r) isa Vector{Pair{Symbol,Any}} + @test Catalyst.hasmetadata(r, :md_1) + @test Catalyst.hasmetadata(r, :md_2) + @test Catalyst.hasmetadata(r, :md_3) + @test Catalyst.hasmetadata(r, :md_4) + @test Catalyst.hasmetadata(r, :md_5) + @test Catalyst.hasmetadata(r, :md_6) + @test !Catalyst.hasmetadata(r, :md_8) - @test isequal(Catalyst.get_metadata(r, :md_1), 1.0) - @test isequal(Catalyst.get_metadata(r, :md_2), false) - @test isequal(Catalyst.get_metadata(r, :md_3), "Hello world") - @test isequal(Catalyst.get_metadata(r, :md_4), :sym) - @test isequal(Catalyst.get_metadata(r, :md_5), X + X2^k -1) - @test isequal(Catalyst.get_metadata(r, :md_6), (0.1, 2.0)) + @test isequal(Catalyst.getmetadata(r, :md_1), 1.0) + @test isequal(Catalyst.getmetadata(r, :md_2), false) + @test isequal(Catalyst.getmetadata(r, :md_3), "Hello world") + @test isequal(Catalyst.getmetadata(r, :md_4), :sym) + @test isequal(Catalyst.getmetadata(r, :md_5), X + X2^k -1) + @test isequal(Catalyst.getmetadata(r, :md_6), (0.1, 2.0)) end # Noise scaling metadata. @@ -85,6 +85,6 @@ let r1 = Reaction(k, [X], [X2], [2], [1]) r2 = Reaction(k, [X], [X2], [2], [1]; metadata=metadata) - @test isequal(Catalyst.getnoisescaling(r1), 1.0) - @test isequal(Catalyst.getnoisescaling(r2), η) + @test isequal(Catalyst.get_noise_scaling(r1), 1.0) + @test isequal(Catalyst.get_noise_scaling(r2), η) end \ No newline at end of file diff --git a/test/reactionsystem_structure/reactionsystem.jl b/test/reactionsystem_structure/reactionsystem.jl index a977143ce7..6a51b9b5d8 100644 --- a/test/reactionsystem_structure/reactionsystem.jl +++ b/test/reactionsystem_structure/reactionsystem.jl @@ -707,6 +707,6 @@ end # Tests metadata. let - @test isnothing(ModelingToolkit.get_metadata(rs)) + @test isnothing(ModelingToolkit.getmetadata(rs)) end From 53f19e1df2c70d25462f8044003c6618f010d6fb Mon Sep 17 00:00:00 2001 From: Torkel Date: Mon, 18 Mar 2024 18:01:27 -0400 Subject: [PATCH 29/41] up --- test/reactionsystem_structure/reactionsystem.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/reactionsystem_structure/reactionsystem.jl b/test/reactionsystem_structure/reactionsystem.jl index 6a51b9b5d8..a977143ce7 100644 --- a/test/reactionsystem_structure/reactionsystem.jl +++ b/test/reactionsystem_structure/reactionsystem.jl @@ -707,6 +707,6 @@ end # Tests metadata. let - @test isnothing(ModelingToolkit.getmetadata(rs)) + @test isnothing(ModelingToolkit.get_metadata(rs)) end From 1327b5f0535ceb38b8544ae27e38d780b3bec303 Mon Sep 17 00:00:00 2001 From: Torkel Loman Date: Fri, 22 Mar 2024 08:44:24 -0400 Subject: [PATCH 30/41] Update HISTORY.md Co-authored-by: Sam Isaacson --- HISTORY.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/HISTORY.md b/HISTORY.md index 5921f0193c..4f0cfc55b8 100644 --- a/HISTORY.md +++ b/HISTORY.md @@ -24,7 +24,7 @@ rn = @reaction_network begin end get_noise_scaling(rn) ``` -- `SDEProblem` no longer takes the `noise_scaling` argument (she above for new approach to handle noise scaling). +- `SDEProblem` no longer takes the `noise_scaling` argument (see above for new approach to handle noise scaling). - Changed fields of internal `Reaction` structure. `ReactionSystems`s saved using `serialize` on previous Catalyst versions cannot be loaded using this (or later) versions. - Simulation of spatial ODEs now supported. For full details, please see https://github.com/SciML/Catalyst.jl/pull/644 and upcoming documentation. Note that these methods are currently considered alpha, with the interface and approach changing even in non-breaking Catalyst releases. - LatticeReactionSystem structure represents a spatial reaction network: From e144ea409ce02f4df5f76ebfb6b81c7e0ba4de3e Mon Sep 17 00:00:00 2001 From: Torkel Loman Date: Fri, 22 Mar 2024 08:44:35 -0400 Subject: [PATCH 31/41] Update docs/src/catalyst_applications/advanced_simulations.md Co-authored-by: Sam Isaacson --- docs/src/catalyst_applications/advanced_simulations.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/src/catalyst_applications/advanced_simulations.md b/docs/src/catalyst_applications/advanced_simulations.md index 7fbf729d10..b0676bb320 100644 --- a/docs/src/catalyst_applications/advanced_simulations.md +++ b/docs/src/catalyst_applications/advanced_simulations.md @@ -292,7 +292,7 @@ plot(sol) ## Scaling the noise magnitude in the chemical Langevin equations When using the CLE to generate SDEs from a CRN, it can sometimes be desirable to -scale the magnitude of the noise terms. Here, each reaction of the system generates a separate noise term in the CLE. If you require identical scaling for all reaction, the `@default_noise_scaling` option can be used. Else, you can supply a `noise_scaling` metadata for each individual reaction, describing how to scale the noise for that reaction. +scale the magnitude of the noise terms. Here, each reaction of the system generates a separate noise term in the CLE. If you require identical scaling for all reactions, the `@default_noise_scaling` option can be used. Else, you can supply a `noise_scaling` metadata for each individual reaction, describing how to scale the noise for that reaction. We begin with considering the first approach. First, we simulate a simple two-state CRN model using the CLE: From dcc092d9a60d329d5d86069532fe14d38f15be08 Mon Sep 17 00:00:00 2001 From: Torkel Loman Date: Fri, 22 Mar 2024 08:44:49 -0400 Subject: [PATCH 32/41] Update docs/src/catalyst_applications/advanced_simulations.md Co-authored-by: Sam Isaacson --- docs/src/catalyst_applications/advanced_simulations.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/src/catalyst_applications/advanced_simulations.md b/docs/src/catalyst_applications/advanced_simulations.md index b0676bb320..ee2b52b4e2 100644 --- a/docs/src/catalyst_applications/advanced_simulations.md +++ b/docs/src/catalyst_applications/advanced_simulations.md @@ -312,7 +312,7 @@ plot(sol_1; idxs = :X1, ylimit = (0.0, 20.0)) ``` Here we can see that the $X$ concentration fluctuates around a steady state of $X≈10.0$. -Next, we wish increase the amount of noise with by a factor 2. To do so, we use the `@default_noise_scaling` option, to which we provide the desired scaling +Next, we wish increase the amount of noise by a factor 2. To do so, we use the `@default_noise_scaling` option, to which we provide the desired scaling ```@example ex3 rn_2 = @reaction_network begin @default_noise_scaling 2 From 8f48831ea86ba6a1eed59cdc64a299bcda1df93e Mon Sep 17 00:00:00 2001 From: Torkel Date: Fri, 22 Mar 2024 09:52:20 -0400 Subject: [PATCH 33/41] doc update --- docs/src/catalyst_applications/advanced_simulations.md | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/docs/src/catalyst_applications/advanced_simulations.md b/docs/src/catalyst_applications/advanced_simulations.md index ee2b52b4e2..2df209f660 100644 --- a/docs/src/catalyst_applications/advanced_simulations.md +++ b/docs/src/catalyst_applications/advanced_simulations.md @@ -319,7 +319,7 @@ rn_2 = @reaction_network begin (k1,k2), X1 <--> X2 end ``` -If we re-simualte the system we see that the amount of noise have increased: +If we re-simulate the system we see that the amount of noise have increased: ```@example ex3 sprob_1 = SDEProblem(rn_2, u0, tspan, p_1) sol_1 = solve(sprob_1, ImplicitEM()) @@ -345,7 +345,7 @@ plot(sol_3; idxs = :X1, ylimit = (0.0, 20.0)) ``` Here we saw how, by setting a small $η$ value, the amount of noise was reduced. -It is possible to use a different noise scaling expression for each reaction. Here, each reaction's noise scaling expression is provided using the `noise_scaling` metadata. In the following example, we use this to turn the noise of for both reactions involving the species $Y$. +It is possible to use a different noise scaling expression for each reaction. Here, each reaction's noise scaling expression is provided using the `noise_scaling` metadata. In the following example, we use this to tune the noise of for both reactions involving the species $Y$. ```@example ex3 rn_4 = @reaction_network begin @@ -361,7 +361,7 @@ sprob_4 = SDEProblem(rn_4, u0_4, tspan, p_4) sol_4 = solve(sprob_4, ImplicitEM()) plot(sol_4; ylimit = (0.0, 20.0)) ``` -Here, we not that there is n fluctuation in the value of $Y$. If the `@default_noise_scaling` option is used, its value is used for all reactions for which the `noise_scaling` metadata is unused. If `@default_noise_scaling` is not used, teh default noise scaling value is `1.0`. +Here, we not that there is n fluctuation in the value of $Y$. If the `@default_noise_scaling` option is used, its value is used for all reactions for which the `noise_scaling` metadata is unused. If `@default_noise_scaling` is not used, the default noise scaling value is `1.0` (i.e. no scaling). ## Useful plotting options Catalyst, just like DifferentialEquations, uses the Plots package for all From 3e37fc4cfdcc7fc85e451892769873f893453c0a Mon Sep 17 00:00:00 2001 From: Torkel Date: Wed, 27 Mar 2024 18:16:26 -0400 Subject: [PATCH 34/41] updates --- src/expression_utils.jl | 10 ++++++++++ src/reaction_network.jl | 40 ++++++++++++++-------------------------- src/reactionsystem.jl | 33 +++++++++++++++++++++++++++++++++ 3 files changed, 57 insertions(+), 26 deletions(-) diff --git a/src/expression_utils.jl b/src/expression_utils.jl index 28ed668704..0ca7bbdeb9 100644 --- a/src/expression_utils.jl +++ b/src/expression_utils.jl @@ -93,6 +93,16 @@ function is_escaped_expr(expr) return (expr isa Expr) && (expr.head == :escape) && (length(expr.args) == 1) end + +### Generic Expression Handling ### + +# Convert an expression with equal signs (e.g. :(a=1.0, b=2.0)) to one with pairs (e.g. :(a=>1.0, b=>2.0)) +function expr_equal_vector_to_pairs(expr_vec) + pair_vector = :([]) + foreach(arg -> push!(pair_vector.args, arg.args[1] => arg.args[2]), expr_vec.args) + return pair_vector +end + ### Old Stuff ### #This will be called whenever a function stored in funcdict is called. diff --git a/src/reaction_network.jl b/src/reaction_network.jl index a97e47eb01..a2cc9e8790 100644 --- a/src/reaction_network.jl +++ b/src/reaction_network.jl @@ -364,9 +364,10 @@ function make_reaction_system(ex::Expr; name = :(gensym(:ReactionSystem))) default_reaction_metadata = :([]) compound_expr, compound_species = read_compound_options(options) check_default_noise_scaling!(default_reaction_metadata, options) + default_reaction_metadata = expr_equal_vector_to_pairs(default_reaction_metadata) # Parses reactions, species, and parameters. - reactions = get_reactions(reaction_lines; default_reaction_metadata) + reactions = get_reactions(reaction_lines) species_declared = [extract_syms(options, :species); compound_species] parameters_declared = extract_syms(options, :parameters) variables = extract_syms(options, :variables) @@ -424,9 +425,12 @@ function make_reaction_system(ex::Expr; name = :(gensym(:ReactionSystem))) $observed_vars $comps - Catalyst.make_ReactionSystem_internal($rxexprs, $tiv, setdiff(union($spssym, $varssym, $compssym), $obs_syms), - $pssym; name = $name, spatial_ivs = $sivs, - observed = $observed_eqs) + Catalyst.remake_ReactionSystem_internal( + Catalyst.make_ReactionSystem_internal( + $rxexprs, $tiv, setdiff(union($spssym, $varssym, $compssym), $obs_syms), + $pssym; name = $name, spatial_ivs = $sivs, observed = $observed_eqs); + default_reaction_metadata = $default_reaction_metadata + ) end end @@ -593,8 +597,7 @@ function get_reaction(line) end # Generates a vector containing a number of reaction structures, each containing the information about one reaction. -function get_reactions(exprs::Vector{Expr}, reactions = Vector{ReactionStruct}(undef, 0); - default_reaction_metadata = :([])) +function get_reactions(exprs::Vector{Expr}, reactions = Vector{ReactionStruct}(undef, 0)) for line in exprs # Reads core reaction information. arrow, rate, reaction, metadata = read_reaction_line(line) @@ -604,16 +607,12 @@ function get_reactions(exprs::Vector{Expr}, reactions = Vector{ReactionStruct}(u if typeof(rate) != Expr || rate.head != :tuple error("Error: Must provide a tuple of reaction rates when declaring a bi-directional reaction.") end - push_reactions!(reactions, reaction.args[2], reaction.args[3], rate.args[1], metadata.args[1], arrow; - default_reaction_metadata) - push_reactions!(reactions, reaction.args[3], reaction.args[2], rate.args[2], metadata.args[2], arrow; - default_reaction_metadata) + push_reactions!(reactions, reaction.args[2], reaction.args[3], rate.args[1], metadata.args[1], arrow) + push_reactions!(reactions, reaction.args[3], reaction.args[2], rate.args[2], metadata.args[2], arrow) elseif in(arrow, fwd_arrows) - push_reactions!(reactions, reaction.args[2], reaction.args[3], rate, metadata, arrow; - default_reaction_metadata) + push_reactions!(reactions, reaction.args[2], reaction.args[3], rate, metadata, arrow) elseif in(arrow, bwd_arrows) - push_reactions!(reactions, reaction.args[3], reaction.args[2], rate, metadata, arrow; - default_reaction_metadata) + push_reactions!(reactions, reaction.args[3], reaction.args[2], rate, metadata, arrow) else throw("Malformed reaction, invalid arrow type used in: $(MacroTools.striplines(line))") end @@ -646,7 +645,7 @@ end # Used to create multiple reactions from, for instance, `k, (X,Y) --> 0`. # Handles metadata, e.g. `1.0, Z --> 0, [noise_scaling=η]`. function push_reactions!(reactions::Vector{ReactionStruct}, sub_line::ExprValues, prod_line::ExprValues, - rate::ExprValues, metadata::ExprValues, arrow::Symbol; default_reaction_metadata::Expr) + rate::ExprValues, metadata::ExprValues, arrow::Symbol) # The rates, substrates, products, and metadata may be in a tupple form (e.g. `k, (X,Y) --> 0`). # This finds the length of these tuples (or 1 if not in tuple forms). Errors if lengs inconsistent. lengs = (tup_leng(sub_line), tup_leng(prod_line), tup_leng(rate), tup_leng(metadata)) @@ -658,7 +657,6 @@ function push_reactions!(reactions::Vector{ReactionStruct}, sub_line::ExprValues for i in 1:maximum(lengs) # If the `only_use_rate` metadata was not provided, this has to be infered from the arrow used. metadata_i = get_tup_arg(metadata, i) - merge_metadata!(metadata_i, default_reaction_metadata) if all([arg.args[1] != :only_use_rate for arg in metadata_i.args]) push!(metadata_i.args, :(only_use_rate = $(in(arrow, pure_rate_arrows)))) end @@ -673,16 +671,6 @@ function push_reactions!(reactions::Vector{ReactionStruct}, sub_line::ExprValues end end -# Merge the metadata encoded by a reaction with the default metadata (for all reactions) given by the system. -# If a metadata value is present in both vectors, uses the non-default value. -function merge_metadata!(metadata, default_reaction_metadata) - for entry in default_reaction_metadata.args - if !in(entry.args[1], arg.args[1] for arg in metadata.args) - push!(metadata.args, entry) - end - end -end - function processmult(op, mult, stoich) if (mult isa Number) && (stoich isa Number) op(mult, stoich) diff --git a/src/reactionsystem.jl b/src/reactionsystem.jl index 1f395602a6..da4158153f 100644 --- a/src/reactionsystem.jl +++ b/src/reactionsystem.jl @@ -756,6 +756,39 @@ function ReactionSystem(iv; kwargs...) ReactionSystem(Reaction[], iv, [], []; kwargs...) end + +""" + remake_ReactionSystem_internal(rs::ReactionSystem; + default_reaction_metadata::Vector{Pair{Symbol, T}} = Vector{Pair{Symbol, Any}}()) where {T} + +Takes a reaction systems and reameks it, returning a modified reaction system. Modifications depend +on which arguments are provided. + +Arguments: +- `rs::ReactionSystem`: The `ReactionSystem` which you wish to make a remade version of. +- `default_reaction_metadata::Vector{Pair{Symbol, T}}`: A vector with default `Reaction` metadata values. + Each metadata in each `Reaction` of the updated `ReactionSystem` will have teh value desiganted in + `default_reaction_metadata` (however, `Reaction`s that already have that metadata designated will not + have their value updated). +""" +#function remake_ReactionSystem_internal(rs::ReactionSystem; default_reaction_metadata::Vector{Pair{Symbol, T}} = Vector{Pair{Symbol, Any}}()) where {T} +function remake_ReactionSystem_internal(rs; default_reaction_metadata = default_reaction_metadata) + # Updates the metadata for all reactions (reactions are ignored). + updated_reactions = [set_default_metadata(rx, default_reaction_metadata) for rx in reactions(rs)] + rs = @set rs.rxs = updated_reactions + return rs +end + +# For a `Reaction`, adds missing default metadata values. Equations are passed back unmodified. +function set_default_metadata(rx, default_metadata) + missing_metadata = filter(md -> !in(md[1], entry[1] for entry in rx.metadata), default_metadata) + updated_metadata = vcat(rx.metadata, missing_metadata) + updated_metadata = convert(Vector{Pair{Symbol, Any}}, updated_metadata) + return @set rx.metadata = updated_metadata +end + + + """ isspatial(rn::ReactionSystem) From 5b98f807ae236c48356ca6e3500f8d866ebddb95 Mon Sep 17 00:00:00 2001 From: Torkel Date: Wed, 27 Mar 2024 18:35:57 -0400 Subject: [PATCH 35/41] up --- src/Catalyst.jl | 1 + src/expression_utils.jl | 4 ++- src/reaction_network.jl | 1 - src/reactionsystem.jl | 34 +++++++++++++++++++------- test/model_simulation/simulate_SDEs.jl | 16 ++++++++++++ 5 files changed, 45 insertions(+), 11 deletions(-) diff --git a/src/Catalyst.jl b/src/Catalyst.jl index 6bd26b35c9..6a64ef7f0b 100644 --- a/src/Catalyst.jl +++ b/src/Catalyst.jl @@ -91,6 +91,7 @@ export make_empty_network, addspecies!, addparam!, addreaction!, reactionparamsm export dependants, dependents, substoichmat, prodstoichmat, netstoichmat export conservationlaws, conservedquantities, conservedequations, conservationlaw_constants export isequivalent +export remake_noise_scaling # depreciated functions to remove in future releases export params, numparams diff --git a/src/expression_utils.jl b/src/expression_utils.jl index 0ca7bbdeb9..5effe48ed8 100644 --- a/src/expression_utils.jl +++ b/src/expression_utils.jl @@ -96,7 +96,9 @@ end ### Generic Expression Handling ### -# Convert an expression with equal signs (e.g. :(a=1.0, b=2.0)) to one with pairs (e.g. :(a=>1.0, b=>2.0)) +# Convert an expression that is a vector with symbols that have values assigned using `=` to vector +# where the assignment instead uses pairs (e.g. :([a=1.0, b=2.0])). Used to e.g. convert default +# reaction metadata to a form that can be evaluated as actual code. function expr_equal_vector_to_pairs(expr_vec) pair_vector = :([]) foreach(arg -> push!(pair_vector.args, arg.args[1] => arg.args[2]), expr_vec.args) diff --git a/src/reaction_network.jl b/src/reaction_network.jl index a2cc9e8790..68061e8596 100644 --- a/src/reaction_network.jl +++ b/src/reaction_network.jl @@ -643,7 +643,6 @@ end # Takes a reaction line and creates reaction(s) from it and pushes those to the reaction array. # Used to create multiple reactions from, for instance, `k, (X,Y) --> 0`. -# Handles metadata, e.g. `1.0, Z --> 0, [noise_scaling=η]`. function push_reactions!(reactions::Vector{ReactionStruct}, sub_line::ExprValues, prod_line::ExprValues, rate::ExprValues, metadata::ExprValues, arrow::Symbol) # The rates, substrates, products, and metadata may be in a tupple form (e.g. `k, (X,Y) --> 0`). diff --git a/src/reactionsystem.jl b/src/reactionsystem.jl index da4158153f..c3bd5f738c 100644 --- a/src/reactionsystem.jl +++ b/src/reactionsystem.jl @@ -761,32 +761,48 @@ end remake_ReactionSystem_internal(rs::ReactionSystem; default_reaction_metadata::Vector{Pair{Symbol, T}} = Vector{Pair{Symbol, Any}}()) where {T} -Takes a reaction systems and reameks it, returning a modified reaction system. Modifications depend -on which arguments are provided. +Takes a `ReactionSystem` and remakes it, returning a modified `ReactionSystem`. Modifications depend +on which additional arguments are provided. The input `ReactionSystem` is not mutated. Arguments: -- `rs::ReactionSystem`: The `ReactionSystem` which you wish to make a remade version of. +- `rs::ReactionSystem`: The `ReactionSystem` which you wish to remake. - `default_reaction_metadata::Vector{Pair{Symbol, T}}`: A vector with default `Reaction` metadata values. - Each metadata in each `Reaction` of the updated `ReactionSystem` will have teh value desiganted in + Each metadata in each `Reaction` of the updated `ReactionSystem` will have the value desiganted in `default_reaction_metadata` (however, `Reaction`s that already have that metadata designated will not have their value updated). """ -#function remake_ReactionSystem_internal(rs::ReactionSystem; default_reaction_metadata::Vector{Pair{Symbol, T}} = Vector{Pair{Symbol, Any}}()) where {T} -function remake_ReactionSystem_internal(rs; default_reaction_metadata = default_reaction_metadata) - # Updates the metadata for all reactions (reactions are ignored). +function remake_ReactionSystem_internal(rs::ReactionSystem; default_reaction_metadata = []) + # Updates the metadata for all reactions (equation are ignored). + updated_equations = [set_default_metadata(rx, default_reaction_metadata) for rx in reactions(rs)] updated_reactions = [set_default_metadata(rx, default_reaction_metadata) for rx in reactions(rs)] - rs = @set rs.rxs = updated_reactions + @set! rs.eqs = updated_equations + @set! rs.rxs = updated_reactions return rs end # For a `Reaction`, adds missing default metadata values. Equations are passed back unmodified. -function set_default_metadata(rx, default_metadata) +function set_default_metadata(rx::Reaction, default_metadata) missing_metadata = filter(md -> !in(md[1], entry[1] for entry in rx.metadata), default_metadata) updated_metadata = vcat(rx.metadata, missing_metadata) updated_metadata = convert(Vector{Pair{Symbol, Any}}, updated_metadata) return @set rx.metadata = updated_metadata end +set_default_metadata(eq::Equation, default_metadata) = eq +""" + remake_noise_scaling(rs::ReactionSystem, noise_scaling) + +Creates an updated `ReactionSystem`. This is the old `ReactionSystem`, but each `Reaction` that does +not have a `Noise_scaling` metadata have its noise_scaling metadata updated. The input `ReactionSystem` +is not mutated. + +Arguments: +- `rs::ReactionSystem`: The `ReactionSystem` which you wish to remake. +- `noise_scaling`: The updated noise scaling terms +""" +function remake_noise_scaling(rs::ReactionSystem, noise_scaling) + return remake_ReactionSystem_internal(rs, default_reaction_metadata = [:noise_scaling => noise_scaling]) +end """ diff --git a/test/model_simulation/simulate_SDEs.jl b/test/model_simulation/simulate_SDEs.jl index a6018fbd4d..5101e5aa0c 100644 --- a/test/model_simulation/simulate_SDEs.jl +++ b/test/model_simulation/simulate_SDEs.jl @@ -234,6 +234,22 @@ let @test var(sol[:X1]) > var(sol[:X2]) > var(sol[:X3]) > var(sol[:X4]) > var(sol[:X5]) end +# Tests the `remake_noise_scaling` function. +let + # Creates noise scaling networks. + noise_scaling_network1 = @reaction_network begin + p, 0 --> X, [noise_scaling=2.0] + d, X --> 0 + end + noise_scaling_network2 = remake_noise_scaling(noise_scaling_network1, 0.5) + + # Checks that the two networks' reactions have the correct metadata. + @test reactions(noise_scaling_network1)[1].metadata == [:noise_scaling => 2.0] + @test reactions(noise_scaling_network1)[2].metadata == [] + @test reactions(noise_scaling_network2)[1].metadata == [:noise_scaling => 2.0] + @test reactions(noise_scaling_network2)[2].metadata == [:noise_scaling => 0.5] +end + ### Checks Simulations Don't Error ### From 69d7250fb587e63e889b602b14411291b702b6fe Mon Sep 17 00:00:00 2001 From: Torkel Date: Wed, 27 Mar 2024 18:38:58 -0400 Subject: [PATCH 36/41] up --- src/expression_utils.jl | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/src/expression_utils.jl b/src/expression_utils.jl index 5effe48ed8..d5fafe64a2 100644 --- a/src/expression_utils.jl +++ b/src/expression_utils.jl @@ -96,9 +96,10 @@ end ### Generic Expression Handling ### -# Convert an expression that is a vector with symbols that have values assigned using `=` to vector -# where the assignment instead uses pairs (e.g. :([a=1.0, b=2.0])). Used to e.g. convert default -# reaction metadata to a form that can be evaluated as actual code. +# Convert an expression that is a vector with symbols that have values assigned using `=` +# (e.g. :([a=1.0, b=2.0])) to a vector where the assignment instead uses symbols and pairs +# (e.g. :([a=>1.0, b=>2.0])). Used to e.g. convert default reaction metadata to a form that can be +# evaluated as actual code. function expr_equal_vector_to_pairs(expr_vec) pair_vector = :([]) foreach(arg -> push!(pair_vector.args, arg.args[1] => arg.args[2]), expr_vec.args) From 2c716e81722e1f0acccdb371ecdaa469e19ed876 Mon Sep 17 00:00:00 2001 From: Torkel Loman Date: Thu, 28 Mar 2024 15:31:26 -0400 Subject: [PATCH 37/41] Update src/reactionsystem.jl Co-authored-by: Sam Isaacson --- src/reactionsystem.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/reactionsystem.jl b/src/reactionsystem.jl index c3bd5f738c..f4cd148a5c 100644 --- a/src/reactionsystem.jl +++ b/src/reactionsystem.jl @@ -773,10 +773,10 @@ Arguments: """ function remake_ReactionSystem_internal(rs::ReactionSystem; default_reaction_metadata = []) # Updates the metadata for all reactions (equation are ignored). - updated_equations = [set_default_metadata(rx, default_reaction_metadata) for rx in reactions(rs)] - updated_reactions = [set_default_metadata(rx, default_reaction_metadata) for rx in reactions(rs)] + eqtransform(eq) = eq isa Reaction ? set_default_metadata(eq, default_reaction_metadata) : eq + updated_equations = map(eqtransform, get_eqs(rs)) @set! rs.eqs = updated_equations - @set! rs.rxs = updated_reactions + @set! rs.rxs = Reaction[rx for rx in updated_equations if rx isa Reaction] return rs end From b41bb7bc81b00ba719a1f64d7aaf526fc9f53efc Mon Sep 17 00:00:00 2001 From: Torkel Date: Thu, 28 Mar 2024 16:01:21 -0400 Subject: [PATCH 38/41] remake --- src/Catalyst.jl | 2 +- src/reactionsystem.jl | 28 ++++++++++++++++++++------ test/model_simulation/simulate_SDEs.jl | 22 ++++++++++++++++++++ 3 files changed, 45 insertions(+), 7 deletions(-) diff --git a/src/Catalyst.jl b/src/Catalyst.jl index 6a64ef7f0b..55efaef221 100644 --- a/src/Catalyst.jl +++ b/src/Catalyst.jl @@ -91,7 +91,7 @@ export make_empty_network, addspecies!, addparam!, addreaction!, reactionparamsm export dependants, dependents, substoichmat, prodstoichmat, netstoichmat export conservationlaws, conservedquantities, conservedequations, conservationlaw_constants export isequivalent -export remake_noise_scaling +export remake_noise_scaling, get_noise_scaling, has_noise_scaling # depreciated functions to remove in future releases export params, numparams diff --git a/src/reactionsystem.jl b/src/reactionsystem.jl index f4cd148a5c..b84a7938cb 100644 --- a/src/reactionsystem.jl +++ b/src/reactionsystem.jl @@ -762,7 +762,8 @@ end default_reaction_metadata::Vector{Pair{Symbol, T}} = Vector{Pair{Symbol, Any}}()) where {T} Takes a `ReactionSystem` and remakes it, returning a modified `ReactionSystem`. Modifications depend -on which additional arguments are provided. The input `ReactionSystem` is not mutated. +on which additional arguments are provided. The input `ReactionSystem` is not mutated. Updating +default reaction metadata is currently the only supported feature. Arguments: - `rs::ReactionSystem`: The `ReactionSystem` which you wish to remake. @@ -772,11 +773,26 @@ Arguments: have their value updated). """ function remake_ReactionSystem_internal(rs::ReactionSystem; default_reaction_metadata = []) - # Updates the metadata for all reactions (equation are ignored). + rs = set_default_metadata(rs; default_reaction_metadata) + return rs +end + +# For a `ReactionSystem`, updates all `Reaction`'s default metadata. +function set_default_metadata(rs::ReactionSystem; default_reaction_metadata = []) + # Updates reaction metadata for for reactions in this specific system. eqtransform(eq) = eq isa Reaction ? set_default_metadata(eq, default_reaction_metadata) : eq updated_equations = map(eqtransform, get_eqs(rs)) @set! rs.eqs = updated_equations @set! rs.rxs = Reaction[rx for rx in updated_equations if rx isa Reaction] + + # Updates reaction metadata for all its subsystems. + new_sub_systems = similar(get_systems(rs)) + for (i, sub_system) in enumerate(get_systems(rs)) + new_sub_systems[i] = set_default_metadata(sub_system; default_reaction_metadata) + end + @set! rs.systems = new_sub_systems + + # Returns the updated system. return rs end @@ -790,17 +806,17 @@ end set_default_metadata(eq::Equation, default_metadata) = eq """ - remake_noise_scaling(rs::ReactionSystem, noise_scaling) +set_default_noise_scaling(rs::ReactionSystem, noise_scaling) Creates an updated `ReactionSystem`. This is the old `ReactionSystem`, but each `Reaction` that does -not have a `Noise_scaling` metadata have its noise_scaling metadata updated. The input `ReactionSystem` -is not mutated. +not have a `noise_scaling` metadata have its noise_scaling metadata updated. The input `ReactionSystem` +is not mutated. Any subsystems of `rs` have their `noise_scaling` metadata updated as well. Arguments: - `rs::ReactionSystem`: The `ReactionSystem` which you wish to remake. - `noise_scaling`: The updated noise scaling terms """ -function remake_noise_scaling(rs::ReactionSystem, noise_scaling) +function set_default_noise_scaling(rs::ReactionSystem, noise_scaling) return remake_ReactionSystem_internal(rs, default_reaction_metadata = [:noise_scaling => noise_scaling]) end diff --git a/test/model_simulation/simulate_SDEs.jl b/test/model_simulation/simulate_SDEs.jl index 5101e5aa0c..29a417f811 100644 --- a/test/model_simulation/simulate_SDEs.jl +++ b/test/model_simulation/simulate_SDEs.jl @@ -250,6 +250,28 @@ let @test reactions(noise_scaling_network2)[2].metadata == [:noise_scaling => 0.5] end +# Tests the `remake_noise_scaling` function on a hierarchical model. +let + # Creates hierarchical model. + rn1 = @reaction_network begin + p, 0 --> X, [noise_scaling=2.0] + d, X --> 0 + end + rn2 = @reaction_network begin + k1, X1 --> X2, [noise_scaling=5.0] + k2, X2 --> X1 + end + rn = compose(rn1, [rn2]) + + # Checks that systems have the correct noise scaling terms. + rn = Catalyst.remake_ReactionSystem_internal(rn, default_reaction_metadata = [:noise_scaling => 0.5]) + rn1_noise_scaling = [get_noise_scaling(rx) for rx in rn.rxs] + rn2_noise_scaling = [get_noise_scaling(rx) for rx in Catalyst.get_systems(rn)[1].rxs] + rn_noise_scaling = [get_noise_scaling(rx) for rx in reactions(rn)] + @test issetequal(rn1_noise_scaling, [2.0, 0.5]) + @test issetequal(rn2_noise_scaling, [5.0, 0.5]) + @test issetequal(rn_noise_scaling, [2.0, 0.5, 5.0, 0.5]) +end ### Checks Simulations Don't Error ### From 1a6d537d854db76fcd7700333d313d36b58297ab Mon Sep 17 00:00:00 2001 From: Torkel Date: Thu, 28 Mar 2024 16:05:39 -0400 Subject: [PATCH 39/41] up --- src/Catalyst.jl | 2 +- test/model_simulation/simulate_SDEs.jl | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/src/Catalyst.jl b/src/Catalyst.jl index 55efaef221..91b403da36 100644 --- a/src/Catalyst.jl +++ b/src/Catalyst.jl @@ -91,7 +91,7 @@ export make_empty_network, addspecies!, addparam!, addreaction!, reactionparamsm export dependants, dependents, substoichmat, prodstoichmat, netstoichmat export conservationlaws, conservedquantities, conservedequations, conservationlaw_constants export isequivalent -export remake_noise_scaling, get_noise_scaling, has_noise_scaling +export set_default_noise_scaling, get_noise_scaling, has_noise_scaling # depreciated functions to remove in future releases export params, numparams diff --git a/test/model_simulation/simulate_SDEs.jl b/test/model_simulation/simulate_SDEs.jl index 29a417f811..e2a812dd7a 100644 --- a/test/model_simulation/simulate_SDEs.jl +++ b/test/model_simulation/simulate_SDEs.jl @@ -241,7 +241,7 @@ let p, 0 --> X, [noise_scaling=2.0] d, X --> 0 end - noise_scaling_network2 = remake_noise_scaling(noise_scaling_network1, 0.5) + noise_scaling_network2 = set_default_noise_scaling(noise_scaling_network1, 0.5) # Checks that the two networks' reactions have the correct metadata. @test reactions(noise_scaling_network1)[1].metadata == [:noise_scaling => 2.0] @@ -264,7 +264,7 @@ let rn = compose(rn1, [rn2]) # Checks that systems have the correct noise scaling terms. - rn = Catalyst.remake_ReactionSystem_internal(rn, default_reaction_metadata = [:noise_scaling => 0.5]) + rn = set_default_noise_scaling(rn, 0.5) rn1_noise_scaling = [get_noise_scaling(rx) for rx in rn.rxs] rn2_noise_scaling = [get_noise_scaling(rx) for rx in Catalyst.get_systems(rn)[1].rxs] rn_noise_scaling = [get_noise_scaling(rx) for rx in reactions(rn)] From de84a61f123c813acfbe896ab95558dde4a5e37e Mon Sep 17 00:00:00 2001 From: Torkel Date: Thu, 28 Mar 2024 16:16:50 -0400 Subject: [PATCH 40/41] use `get_rxs` in test --- test/model_simulation/simulate_SDEs.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/model_simulation/simulate_SDEs.jl b/test/model_simulation/simulate_SDEs.jl index e2a812dd7a..704e2de378 100644 --- a/test/model_simulation/simulate_SDEs.jl +++ b/test/model_simulation/simulate_SDEs.jl @@ -265,8 +265,8 @@ let # Checks that systems have the correct noise scaling terms. rn = set_default_noise_scaling(rn, 0.5) - rn1_noise_scaling = [get_noise_scaling(rx) for rx in rn.rxs] - rn2_noise_scaling = [get_noise_scaling(rx) for rx in Catalyst.get_systems(rn)[1].rxs] + rn1_noise_scaling = [get_noise_scaling(rx) for rx in get_rxs(rn)] + rn2_noise_scaling = [get_noise_scaling(rx) for rx in get_rxs(Catalyst.get_systems(rn)[1])] rn_noise_scaling = [get_noise_scaling(rx) for rx in reactions(rn)] @test issetequal(rn1_noise_scaling, [2.0, 0.5]) @test issetequal(rn2_noise_scaling, [5.0, 0.5]) From b4560e8862858ae6df6737bd6a7427e6e6469cb3 Mon Sep 17 00:00:00 2001 From: Torkel Loman Date: Fri, 29 Mar 2024 01:40:49 +0100 Subject: [PATCH 41/41] Update src/reaction_network.jl Co-authored-by: Sam Isaacson --- src/reaction_network.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/reaction_network.jl b/src/reaction_network.jl index 68061e8596..d2d5b29594 100644 --- a/src/reaction_network.jl +++ b/src/reaction_network.jl @@ -656,7 +656,7 @@ function push_reactions!(reactions::Vector{ReactionStruct}, sub_line::ExprValues for i in 1:maximum(lengs) # If the `only_use_rate` metadata was not provided, this has to be infered from the arrow used. metadata_i = get_tup_arg(metadata, i) - if all([arg.args[1] != :only_use_rate for arg in metadata_i.args]) + if all(arg -> arg.args[1] != :only_use_rate, metadata_i.args) push!(metadata_i.args, :(only_use_rate = $(in(arrow, pure_rate_arrows)))) end