Skip to content
Snippets Groups Projects
Commit e6944e0b authored by Bentriou Mahmoud's avatar Bentriou Mahmoud
Browse files

Implementation of the possibility of naming a model

parent 70407040
No related branches found
No related tags found
No related merge requests found
...@@ -35,7 +35,8 @@ function get_str_propensity(propensity::Symbol, dict_species::Dict, dict_params: ...@@ -35,7 +35,8 @@ function get_str_propensity(propensity::Symbol, dict_species::Dict, dict_params:
return str_propensity return str_propensity
end end
macro biochemical_network(expr_name,expr_network) macro biochemical_network(expr_network,expr_name...)
model_name = isempty(expr_name) ? "Unnamed macro generated" : expr_name[1]
transitions = String[] transitions = String[]
dict_species = Dict{String,Int}() dict_species = Dict{String,Int}()
dict_params = Dict{String,Int}() dict_params = Dict{String,Int}()
...@@ -103,7 +104,7 @@ macro biochemical_network(expr_name,expr_network) ...@@ -103,7 +104,7 @@ macro biochemical_network(expr_name,expr_network)
# Let's write some lines that creates the function f! (step of a simulation) for this biochemical network # Let's write some lines that creates the function f! (step of a simulation) for this biochemical network
nbr_rand = rand(1:1000) nbr_rand = rand(1:1000)
nbr_reactions = length(list_expr_reactions) nbr_reactions = length(list_expr_reactions)
basename_func = "$(replace(expr_name, ' '=>'_'))_$(nbr_rand)" basename_func = "$(replace(model_name, ' '=>'_'))_$(nbr_rand)"
expr_model_f! = "function $(basename_func)_f!(xnplus1::Vector{Int}, l_t::Vector{Float64}, l_tr::Vector{Union{Nothing,String}}, xn::Vector{Int}, tn::Float64, p::Vector{Float64})\n\t" expr_model_f! = "function $(basename_func)_f!(xnplus1::Vector{Int}, l_t::Vector{Float64}, l_tr::Vector{Union{Nothing,String}}, xn::Vector{Int}, tn::Float64, p::Vector{Float64})\n\t"
# Computation of nu and propensity functions in f! # Computation of nu and propensity functions in f!
str_l_a = "l_a = (" str_l_a = "l_a = ("
...@@ -178,6 +179,7 @@ macro biochemical_network(expr_name,expr_network) ...@@ -178,6 +179,7 @@ macro biochemical_network(expr_name,expr_network)
model_f! = eval(Meta.parse(expr_model_f!)) model_f! = eval(Meta.parse(expr_model_f!))
model_isabsorbing = eval(Meta.parse(expr_model_isabsorbing)) model_isabsorbing = eval(Meta.parse(expr_model_isabsorbing))
return :(ContinuousTimeModel($dim_state, $dim_params, $dict_species, $dict_params, $transitions, return :(ContinuousTimeModel($dim_state, $dim_params, $dict_species, $dict_params, $transitions,
$(zeros(dim_params)), $(zeros(Int, dim_state)), 0.0, $model_f!, $model_isabsorbing; g = $list_species)) $(zeros(dim_params)), $(zeros(Int, dim_state)), 0.0, $model_f!, $model_isabsorbing;
g = $list_species, name=$model_name))
end end
...@@ -10,6 +10,7 @@ const Location = String ...@@ -10,6 +10,7 @@ const Location = String
const VariableAutomaton = String const VariableAutomaton = String
mutable struct ContinuousTimeModel <: Model mutable struct ContinuousTimeModel <: Model
name::String
dim_state::Int # state space dim dim_state::Int # state space dim
dim_params::Int # parameter space dim dim_params::Int # parameter space dim
map_var_idx::Dict{String,Int} # maps variable str to index in the state space map_var_idx::Dict{String,Int} # maps variable str to index in the state space
...@@ -85,7 +86,8 @@ end ...@@ -85,7 +86,8 @@ end
function ContinuousTimeModel(dim_state::Int, dim_params::Int, map_var_idx::Dict, map_param_idx::Dict, transitions::Vector{String}, function ContinuousTimeModel(dim_state::Int, dim_params::Int, map_var_idx::Dict, map_param_idx::Dict, transitions::Vector{String},
p::Vector{Float64}, x0::Vector{Int}, t0::Float64, p::Vector{Float64}, x0::Vector{Int}, t0::Float64,
f!::Function, isabsorbing::Function; f!::Function, isabsorbing::Function;
g::Vector{String} = keys(map_var_idx), time_bound::Float64 = Inf, buffer_size::Int = 10, estim_min_states::Int = 50) g::Vector{String} = keys(map_var_idx), time_bound::Float64 = Inf,
buffer_size::Int = 10, estim_min_states::Int = 50, name::String = "Unnamed")
dim_obs_state = length(g) dim_obs_state = length(g)
_map_obs_var_idx = Dict() _map_obs_var_idx = Dict()
_g_idx = Vector{Int}(undef, dim_obs_state) _g_idx = Vector{Int}(undef, dim_obs_state)
...@@ -100,7 +102,8 @@ function ContinuousTimeModel(dim_state::Int, dim_params::Int, map_var_idx::Dict, ...@@ -100,7 +102,8 @@ function ContinuousTimeModel(dim_state::Int, dim_params::Int, map_var_idx::Dict,
if length(methods(isabsorbing)) >= 2 if length(methods(isabsorbing)) >= 2
@warn "You have possibly redefined a function Model.isabsorbing used in a previously instantiated model." @warn "You have possibly redefined a function Model.isabsorbing used in a previously instantiated model."
end end
new_model = ContinuousTimeModel(dim_state, dim_params, map_var_idx, _map_obs_var_idx, map_param_idx, transitions, p, x0, t0, f!, g, _g_idx, isabsorbing, time_bound, buffer_size, estim_min_states) new_model = ContinuousTimeModel(name, dim_state, dim_params, map_var_idx, _map_obs_var_idx, map_param_idx, transitions,
p, x0, t0, f!, g, _g_idx, isabsorbing, time_bound, buffer_size, estim_min_states)
@assert check_consistency(new_model) @assert check_consistency(new_model)
return new_model return new_model
end end
......
...@@ -364,6 +364,23 @@ function distribute_prob_accept_lha(sm::SynchronizedModel, nbr_sim::Int) ...@@ -364,6 +364,23 @@ function distribute_prob_accept_lha(sm::SynchronizedModel, nbr_sim::Int)
return sum_val / nbr_sim return sum_val / nbr_sim
end end
function Base.show(io::IO, m::ContinuousTimeModel)
print(io, "$(m.name) model\n")
print(io, "- variables :\n")
for (var, idx) in m.map_var_idx
print(io, "* $var (idx = $idx in state space)\n")
end
print(io, "- parameters :\n")
for (param, idx) in m.map_param_idx
print(io, "* $param (idx = $idx in parameter space)\n")
end
print(io, "- transitions : $(m.transitions)\n")
print(io, "- observed variables :\n")
for i in eachindex(m.g)
print(io, "* $(m.g[i]) (idx = $i in observed state space)\n")
end
end
isbounded(m::Model) = get_proba_model(m).time_bound < Inf isbounded(m::Model) = get_proba_model(m).time_bound < Inf
function check_consistency(m::ContinuousTimeModel) function check_consistency(m::ContinuousTimeModel)
@assert m.dim_state == length(m.map_var_idx) @assert m.dim_state == length(m.map_var_idx)
......
...@@ -52,7 +52,7 @@ isabsorbing_ER(p::Vector{Float64},xn::Vector{Int}) = ...@@ -52,7 +52,7 @@ isabsorbing_ER(p::Vector{Float64},xn::Vector{Int}) =
(p[1]*xn[1]*xn[2] + (p[2]+p[3])*xn[3]) === 0.0 (p[1]*xn[1]*xn[2] + (p[2]+p[3])*xn[3]) === 0.0
g_ER = ["P"] g_ER = ["P"]
ER = ContinuousTimeModel(d,k,dict_var_ER,dict_p_ER,l_tr_ER,p_ER,x0_ER,t0_ER,ER_f!,isabsorbing_ER; g=g_ER) ER = ContinuousTimeModel(d,k,dict_var_ER,dict_p_ER,l_tr_ER,p_ER,x0_ER,t0_ER,ER_f!,isabsorbing_ER; g=g_ER,name="ER pkg")
function create_ER(new_p::Vector{Float64}) function create_ER(new_p::Vector{Float64})
ER_new = deepcopy(ER) ER_new = deepcopy(ER)
......
...@@ -50,7 +50,7 @@ end ...@@ -50,7 +50,7 @@ end
isabsorbing_SIR(p::Vector{Float64}, xn::Vector{Int}) = (p[1]*xn[1]*xn[2] + p[2]*xn[2]) === 0.0 isabsorbing_SIR(p::Vector{Float64}, xn::Vector{Int}) = (p[1]*xn[1]*xn[2] + p[2]*xn[2]) === 0.0
g_SIR = ["I"] g_SIR = ["I"]
SIR = ContinuousTimeModel(d,k,dict_var,dict_p,l_tr_SIR,p_SIR,x0_SIR,t0_SIR,SIR_f!,isabsorbing_SIR; g=g_SIR) SIR = ContinuousTimeModel(d,k,dict_var,dict_p,l_tr_SIR,p_SIR,x0_SIR,t0_SIR,SIR_f!,isabsorbing_SIR; g=g_SIR, name="SIR pkg")
function create_SIR(new_p::Vector{Float64}) function create_SIR(new_p::Vector{Float64})
SIR_new = deepcopy(SIR) SIR_new = deepcopy(SIR)
......
...@@ -37,7 +37,7 @@ g_SIR_tauleap = ["I"] ...@@ -37,7 +37,7 @@ g_SIR_tauleap = ["I"]
SIR_tauleap = ContinuousTimeModel(d,k,dict_var_SIR_tauleap,dict_p_SIR_tauleap,l_tr_SIR_tauleap, SIR_tauleap = ContinuousTimeModel(d,k,dict_var_SIR_tauleap,dict_p_SIR_tauleap,l_tr_SIR_tauleap,
p_SIR_tauleap,x0_SIR_tauleap,t0_SIR_tauleap, p_SIR_tauleap,x0_SIR_tauleap,t0_SIR_tauleap,
SIR_tauleap_f!,isabsorbing_SIR_tauleap; g=g_SIR_tauleap) SIR_tauleap_f!,isabsorbing_SIR_tauleap; g=g_SIR_tauleap, name="SIR tauleap pkg")
function create_SIR_tauleap(new_p::Vector{Float64}) function create_SIR_tauleap(new_p::Vector{Float64})
SIR_tauleap_new = deepcopy(SIR_tauleap) SIR_tauleap_new = deepcopy(SIR_tauleap)
@assert length(SIR_tauleap_new.p) == length(new_p) @assert length(SIR_tauleap_new.p) == length(new_p)
......
...@@ -24,7 +24,8 @@ g_poisson = ["N"] ...@@ -24,7 +24,8 @@ g_poisson = ["N"]
poisson = ContinuousTimeModel(d,k,dict_var_poisson,dict_p_poisson,l_tr_poisson, poisson = ContinuousTimeModel(d,k,dict_var_poisson,dict_p_poisson,l_tr_poisson,
p_poisson,x0_poisson,t0_poisson, p_poisson,x0_poisson,t0_poisson,
poisson_f!,isabsorbing_poisson; g=g_poisson, time_bound=1.0) poisson_f!,isabsorbing_poisson;
g=g_poisson, time_bound=1.0, name="Poisson process pkg")
function create_poisson(new_p::Vector{Float64}) function create_poisson(new_p::Vector{Float64})
poisson_new = deepcopy(poisson) poisson_new = deepcopy(poisson)
@assert length(poisson_new.p) == length(new_p) @assert length(poisson_new.p) == length(new_p)
......
using MarkovProcesses using MarkovProcesses
model_SIR = @biochemical_network "SIR" begin model_SIR = @biochemical_network begin
R1: (S+I => 2I, ki*S*I) R1: (S+I => 2I, ki*S*I)
R2: (I => R, kr*I) R2: (I => R, kr*I)
end end "SIR"
set_x0!(model_SIR, [95,5,0]) set_x0!(model_SIR, [95,5,0])
set_param!(model_SIR, [0.012, 0.05]) set_param!(model_SIR, [0.012, 0.05])
model_ER = @biochemical_network "ER" begin model_unnamed_SIR = @biochemical_network begin
R1: (S+I => 2I, ki*S*I)
R2: (I => R, kr*I)
end
set_x0!(model_unnamed_SIR, [95,5,0])
set_param!(model_unnamed_SIR, [0.012, 0.05])
model_ER = @biochemical_network begin
R1: (E+S => ES, k1*E*S) R1: (E+S => ES, k1*E*S)
R2: (ES => E+S, k2*ES) R2: (ES => E+S, k2*ES)
R3: (ES => E+P, k3*ES) R3: (ES => E+P, k3*ES)
end end "ER"
set_x0!(model_ER, [100,100,0,0]) set_x0!(model_ER, [100,100,0,0])
set_param!(model_ER, [1.0,1.0,1.0]) set_param!(model_ER, [1.0,1.0,1.0])
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment