Skip to content

Commit

Permalink
update
Browse files Browse the repository at this point in the history
  • Loading branch information
@TorkelE
TorkelE committed Sep 18, 2023
1 parent 16b9d35 commit 5e805b9
Show file tree
Hide file tree
Showing 5 changed files with 22 additions and 68 deletions.
11 changes: 3 additions & 8 deletions src/networkapi.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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 ###############################
Expand Down
8 changes: 4 additions & 4 deletions src/reaction_network.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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]
Expand Down
7 changes: 6 additions & 1 deletion src/reactionsystem.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down Expand Up @@ -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,
Expand Down
44 changes: 9 additions & 35 deletions test/model_simulation/simulate_SDEs.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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

Expand All @@ -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.
Expand All @@ -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 ###

Expand Down
20 changes: 0 additions & 20 deletions test/model_simulation/u0_n_parameter_inputs.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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

0 comments on commit 5e805b9

Please sign in to comment.