Skip to content

Commit

Permalink
improved noise scaling
Browse files Browse the repository at this point in the history
  • Loading branch information
TorkelE committed Sep 17, 2023
1 parent af70af1 commit aa2fe23
Show file tree
Hide file tree
Showing 3 changed files with 69 additions and 9 deletions.
26 changes: 21 additions & 5 deletions src/networkapi.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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!!!
Expand Down Expand Up @@ -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 ###############################
Expand Down
39 changes: 37 additions & 2 deletions src/reaction_network.jl
Original file line number Diff line number Diff line change
Expand Up @@ -123,7 +123,7 @@ const forbidden_variables_error = let
end

# Declares the keys used for various options.
const option_keys = (:species, :parameters, :variables, :ivs)
const option_keys = (:species, :parameters, :variables, :ivs, :noise_scaling_parameters)

### The @species macro, basically a copy of the @variables macro. ###
macro species(ex...)
Expand Down Expand Up @@ -364,6 +364,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))
Expand All @@ -382,7 +385,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.")
Expand All @@ -404,6 +407,9 @@ 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)

# Returns the rephrased expression.
quote
$ps
Expand Down Expand Up @@ -531,8 +537,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)
Expand Down Expand Up @@ -676,6 +696,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)
Expand Down
13 changes: 11 additions & 2 deletions src/reactionsystem.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down Expand Up @@ -1427,7 +1427,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

"""
Expand Down

0 comments on commit aa2fe23

Please sign in to comment.