From 5e805b91674b64d93be4ca9a614e81787f5b955b Mon Sep 17 00:00:00 2001 From: Torkel Date: Mon, 18 Sep 2023 18:45:39 -0400 Subject: [PATCH] update --- src/networkapi.jl | 11 ++--- src/reaction_network.jl | 8 ++-- src/reactionsystem.jl | 7 ++- test/model_simulation/simulate_SDEs.jl | 44 ++++--------------- .../model_simulation/u0_n_parameter_inputs.jl | 20 --------- 5 files changed, 22 insertions(+), 68 deletions(-) diff --git a/src/networkapi.jl b/src/networkapi.jl index e0a3341a37..3ec3359669 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 e85a3f4e30..7e2da25989 100644 --- a/src/reaction_network.jl +++ b/src/reaction_network.jl @@ -365,7 +365,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) @@ -407,8 +407,8 @@ function make_reaction_system(ex::Expr; name = :(gensym(:ReactionSystem))) push!(rxexprs.args, get_rxexprs(reaction)) end - # Adds potentital nosie scaling parameters. - #append!(ps.args, nosie_scaling_p_args) + # Adds potentital noise scaling parameters. + #append!(ps.args, noise_scaling_p_args) # Returns the rephrased expression. quote @@ -540,7 +540,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 0a5f935c79..b49cb41f28 100644 --- a/src/reactionsystem.jl +++ b/src/reactionsystem.jl @@ -1387,6 +1387,7 @@ function Base.convert(::Type{<:SDESystem}, rs::ReactionSystem; include_zero_odes = true, checks = false, remove_conserved = false, 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) @@ -1520,11 +1521,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 8fbe2e7b37..a9716402a1 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 0068fd8616..b5d556a23d 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