Skip to content

Commit

Permalink
update
Browse files Browse the repository at this point in the history
  • Loading branch information
TorkelE committed Sep 19, 2023
1 parent c5824ee commit 916b1cb
Show file tree
Hide file tree
Showing 5 changed files with 63 additions and 32 deletions.
2 changes: 1 addition & 1 deletion src/Catalyst.jl
Original file line number Diff line number Diff line change
Expand Up @@ -65,7 +65,7 @@ export ODEProblem,
const ExprValues = Union{Expr, Symbol, Float64, Int}
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")
Expand Down
1 change: 1 addition & 0 deletions src/networkapi.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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!!!
Expand Down
73 changes: 52 additions & 21 deletions src/reaction_network.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down Expand Up @@ -364,6 +365,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))
Expand Down Expand Up @@ -396,6 +400,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")
Expand All @@ -404,17 +409,20 @@ 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

Catalyst.make_ReactionSystem_internal($rxexprs, $tiv, union($spssym, $varssym),
$pssym; name = $name,
union($pssym, $ns_pssym); name = $name,
spatial_ivs = $sivs)
end
end
Expand Down Expand Up @@ -533,26 +541,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)
Expand Down Expand Up @@ -683,6 +673,47 @@ function recursive_find_reactants!(ex::ExprValues, mult::ExprValues,
reactants
end

### Handle Options ###

# macro, similar to @parameters, but all paraemters become noise scaling parameters.
macro noise_scaling_parameters(ex...)
vars = Symbolics._parse_parameters(:parameters, Real, ex)

# vector of symbols that get defined
lastarg = vars.args[end]

# start adding metadata statements where the vector of symbols was previously declared
idx = length(vars.args)
resize!(vars.args, idx + length(lastarg.args) + 1)
for sym in lastarg.args
vars.args[idx] = :($sym = ModelingToolkit.wrap(setmetadata(ModelingToolkit.value($sym),
Catalyst.NoiseScalingParameter,
true)))
idx += 1
end

# check nothing was declared isconstantspecies
ex = quote end
vars.args[idx] = ex
idx += 1

# put back the vector of the new species symbols
vars.args[idx] = lastarg
println(vars)
esc(vars)
end

# Creates an expression listing any noise scaling parameters (and also return a list of them).
function handle_noise_scaling_parameters!(parameters_declared, options)
haskey(options, :noise_scaling_parameters) || return :()

noise_scaling_parameters = extract_syms(options, :noise_scaling_parameters)
isempty(intersect(parameters_declared, noise_scaling_parameters)) || error("Parameters: $(intersect(parameters_declared, noise_scaling_parameters_declared)) were declared using both @parameters and @noise_scaling_parameters, please only use one (the later for nosie scaling parameters, else the latter).")
append!(parameters_declared, noise_scaling_parameters)

return options[:noise_scaling_parameters]
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.
Expand Down
17 changes: 8 additions & 9 deletions src/reactionsystem.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down Expand Up @@ -1387,10 +1387,13 @@ 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)
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)) &&
Expand Down Expand Up @@ -1428,7 +1431,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
Expand Down Expand Up @@ -1526,10 +1529,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,
Expand Down
2 changes: 1 addition & 1 deletion test/model_simulation/simulate_SDEs.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down

0 comments on commit 916b1cb

Please sign in to comment.