From 5d886fc4a6d0297d1d67a27599c699fbb34c9ba4 Mon Sep 17 00:00:00 2001
From: Mahmoud Bentriou <mahmoud.bentriou@centralesupelec.fr>
Date: Sat, 20 Feb 2021 11:31:02 +0100
Subject: [PATCH] The package becomes more meta to reach higher performance.
This commit groups the change operated to the creation of models and simulate
function of a ContinuousTimeModel. The general idea is to create a concrete
type and a simulate function per model creation by metaprogramming. - Now,
ContinuousTimeModel is an abstract type. Each creation of a model defines a
concrete type T <: ContinuousTimeModel by meta programming. - f! and
isabsorbing ContinuousTimeModel fields are Symbols. -
simulate(::ContinuousTimeModel) is run by multiple dispatch, according to the
type of the model.
Can't run the whole tests for now but unit/simulate_available_models.jl
runs properly (i've updated the list of models in this commit), and I've
manually checked in the repl that simulations run correctly (distributed
/ plots).
---
core/MarkovProcesses.jl | 3 +-
core/common.jl | 94 ++--
core/lha.jl | 352 +++++++-------
core/model.jl | 591 ++++++++++++------------
core/network_model.jl | 104 +++--
core/trajectory.jl | 3 +-
core/utils.jl | 6 +-
models/ER.jl | 9 +-
models/SIR.jl | 10 +-
models/SIR_tauleap.jl | 17 +-
models/doping_3way_oscillator.jl | 2 -
models/intracellular_viral_infection.jl | 2 -
models/poisson.jl | 18 +-
models/repressilator.jl | 6 +-
models/square_wave_oscillator.jl | 15 +-
tests/unit/simulate_available_models.jl | 8 +
16 files changed, 653 insertions(+), 587 deletions(-)
diff --git a/core/MarkovProcesses.jl b/core/MarkovProcesses.jl
index a1746b8..bd76447 100644
--- a/core/MarkovProcesses.jl
+++ b/core/MarkovProcesses.jl
@@ -4,9 +4,10 @@ import Base: +, -, *
import Base: copy, getfield, getindex, getproperty, lastindex
import Base: setindex!, setproperty!, fill!, copyto!
-import StaticArrays: SVector
import Distributed: @everywhere, @distributed
import Distributions: Distribution, Product, Uniform, Normal
+import Dates
+import StaticArrays: SVector
export Distribution, Product, Uniform, Normal
diff --git a/core/common.jl b/core/common.jl
index 1623ef4..e7909d5 100644
--- a/core/common.jl
+++ b/core/common.jl
@@ -3,6 +3,7 @@ import Distributions: Distribution, Univariate, Continuous, UnivariateDistributi
MultivariateDistribution, product_distribution
abstract type Model end
+abstract type ContinuousTimeModel <: Model end
abstract type AbstractTrajectory end
const VariableModel = Symbol
@@ -11,24 +12,27 @@ const Transition = Union{Symbol,Nothing}
const Location = Symbol
const VariableAutomaton = Symbol
-mutable struct ContinuousTimeModel <: Model
- name::String
- dim_state::Int # state space dim
- dim_params::Int # parameter space dim
- map_var_idx::Dict{VariableModel,Int} # maps variable str to index in the state space
- _map_obs_var_idx::Dict{VariableModel,Int} # maps variable str to index in the observed state space
- map_param_idx::Dict{ParameterModel,Int} # maps parameter str to index in the parameter space
- transitions::Vector{Transition}
- p::Vector{Float64}
- x0::Vector{Int}
- t0::Float64
- f!::Function
- g::Vector{VariableModel} # of dimension dim_obs_state
- _g_idx::Vector{Int} # of dimension dim_obs_state
- isabsorbing::Function
- time_bound::Float64
- buffer_size::Int
- estim_min_states::Int
+function generate_code_model_type_def(model_name::Symbol)
+ return quote
+ mutable struct $(model_name) <: ContinuousTimeModel
+ dim_state::Int # state space dim
+ dim_params::Int # parameter space dim
+ map_var_idx::Dict{VariableModel,Int} # maps variable str to index in the state space
+ _map_obs_var_idx::Dict{VariableModel,Int} # maps variable str to index in the observed state space
+ map_param_idx::Dict{ParameterModel,Int} # maps parameter str to index in the parameter space
+ transitions::Vector{Transition}
+ p::Vector{Float64}
+ x0::Vector{Int}
+ t0::Float64
+ f!::Symbol
+ g::Vector{VariableModel} # of dimension dim_obs_state
+ _g_idx::Vector{Int} # of dimension dim_obs_state
+ isabsorbing::Symbol
+ time_bound::Float64
+ buffer_size::Int
+ estim_min_states::Int
+ end
+ end
end
struct Trajectory <: AbstractTrajectory
@@ -38,11 +42,14 @@ struct Trajectory <: AbstractTrajectory
transitions::Vector{Transition}
end
+#=
struct Edge
transitions::Union{Nothing,Vector{Symbol}}
check_constraints::Function
update_state!::Function
end
+=#
+abstract type Edge end
struct LHA
name::String
@@ -86,31 +93,34 @@ struct ParametricModel
end
# Constructors
-function ContinuousTimeModel(dim_state::Int, dim_params::Int, map_var_idx::Dict{VariableModel,Int},
- map_param_idx::Dict{ParameterModel,Int}, transitions::Vector{<:Transition},
- p::Vector{Float64}, x0::Vector{Int}, t0::Float64,
- f!::Function, isabsorbing::Function;
- g::Vector{VariableModel} = [var for var in 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)
- transitions = convert(Vector{Transition}, transitions)
- _map_obs_var_idx = Dict()
- _g_idx = Vector{Int}(undef, dim_obs_state)
- for i = 1:dim_obs_state
- _g_idx[i] = map_var_idx[g[i]] # = ( (g[i] = i-th obs var)::VariableModel => idx in state space )
- _map_obs_var_idx[g[i]] = i
- end
-
- if length(methods(f!)) >= 2
- @warn "You have possibly redefined a function Model.f! used in a previously instantiated model."
- end
- if length(methods(isabsorbing)) >= 2
- @warn "You have possibly redefined a function Model.isabsorbing used in a previously instantiated model."
+function generate_code_model_type_constructor(model_name::Symbol)
+ return quote
+ function $(model_name)(dim_state::Int, dim_params::Int, map_var_idx::Dict{VariableModel,Int},
+ map_param_idx::Dict{ParameterModel,Int}, transitions::Vector{<:Transition},
+ p::Vector{Float64}, x0::Vector{Int}, t0::Float64,
+ f!::Symbol, isabsorbing::Symbol;
+ g::Vector{VariableModel} = [var for var in keys(map_var_idx)], time_bound::Float64 = Inf,
+ buffer_size::Int = 10, estim_min_states::Int = 50)
+ dim_obs_state = length(g)
+ transitions = convert(Vector{Transition}, transitions)
+ _map_obs_var_idx = Dict()
+ _g_idx = Vector{Int}(undef, dim_obs_state)
+ for i = 1:dim_obs_state
+ _g_idx[i] = map_var_idx[g[i]] # = ( (g[i] = i-th obs var)::VariableModel => idx in state space )
+ _map_obs_var_idx[g[i]] = i
+ end
+ if length(methods(getfield(Main, f!))) >= 2
+ @warn "You have possibly redefined a function Model.f! used in a previously instantiated model."
+ end
+ if length(methods(getfield(Main, isabsorbing))) >= 2
+ @warn "You have possibly redefined a function Model.isabsorbing used in a previously instantiated model."
+ end
+ new_model = $(model_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)
+ return new_model
+ end
end
- 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)
- return new_model
end
LHA(A::LHA, map_var::Dict{VariableModel,Int}) = LHA(A.name, A.transitions, A.locations, A.Λ,
diff --git a/core/lha.jl b/core/lha.jl
index 0a5297c..54f86ca 100644
--- a/core/lha.jl
+++ b/core/lha.jl
@@ -58,12 +58,12 @@ function Base.show(io::IO, S::StateLHA)
end
function Base.show(io::IO, E::Edge)
- print(io, "(Edge: ")
- print(io, (E.transitions == nothing) ? "Asynchronous #, " : ("Synchronized with " * join(E.transitions,',') * ", "))
- print(io, Symbol(E.check_constraints))
+ print(io, "($(typeof(E)): ")
+ print(io, (E.transitions == nothing) ? "Asynchronous #)" : ("Synchronized with " * join(E.transitions,',') * ")"))
+ #=print(io, Symbol(E.check_constraints))
print(io, ", ")
print(io, Symbol(E.update_state!))
- print(io, ")")
+ print(io, ")")=#
end
function Base.copyto!(Sdest::StateLHA, Ssrc::StateLHA)
@@ -101,154 +101,172 @@ function _push_edge!(edge_candidates::Vector{Edge}, edge::Edge, nbr_candidates::
end
end
-function _find_edge_candidates!(edge_candidates::Vector{Edge},
- edges_from_current_loc::Dict{Location,Vector{Edge}},
- Λ::Dict{Location,Function},
- S_time::Float64, S_values::Vector{Float64},
- x::Vector{Int}, p::Vector{Float64},
- only_asynchronous::Bool)
- nbr_candidates = 0
- for target_loc in keys(edges_from_current_loc)
- if !Λ[target_loc](x) continue end
- for edge in edges_from_current_loc[target_loc]
- if getfield(edge, :check_constraints)(S_time, S_values, x, p)
- if getfield(edge, :transitions) == nothing
- _push_edge!(edge_candidates, edge, nbr_candidates)
- nbr_candidates += 1
- return nbr_candidates
- else
- if !only_asynchronous
- _push_edge!(edge_candidates, edge, nbr_candidates)
- nbr_candidates += 1
+function generate_next_state_code()
+ return quote
+ function _find_edge_candidates!(edge_candidates::Vector{Edge},
+ edges_from_current_loc::Dict{Location,Vector{Edge}},
+ Λ::Dict{Location,Function},
+ S_time::Float64, S_values::Vector{Float64},
+ x::Vector{Int}, p::Vector{Float64},
+ only_asynchronous::Bool)
+ nbr_candidates = 0
+ for target_loc in keys(edges_from_current_loc)
+ if !Λ[target_loc](x) continue end
+ for edge in edges_from_current_loc[target_loc]
+ if check_constraints(edge, S_time, S_values, x, p)
+ if getfield(edge, :transitions) == nothing
+ MarkovProcesses._push_edge!(edge_candidates, edge, nbr_candidates)
+ nbr_candidates += 1
+ return nbr_candidates
+ else
+ if !only_asynchronous
+ MarkovProcesses._push_edge!(edge_candidates, edge, nbr_candidates)
+ nbr_candidates += 1
+ end
+ end
end
end
end
+ return nbr_candidates
end
- end
- return nbr_candidates
-end
-function _get_edge_index(edge_candidates::Vector{Edge}, nbr_candidates::Int,
- detected_event::Bool, tr_nplus1::Transition)
- ind_edge = 0
- bool_event = detected_event
- for i = 1:nbr_candidates
- edge = edge_candidates[i]
- # Asynchronous edge detection: we fire it
- if getfield(edge, :transitions) == nothing
- return (i, detected_event)
- end
- # Synchronous detection
- if !detected_event && tr_nplus1 != nothing
- if (getfield(edge, :transitions)[1] == :ALL) ||
- (tr_nplus1 in getfield(edge, :transitions))
- ind_edge = i
- bool_event = true
+ function _get_edge_index(edge_candidates::Vector{Edge}, nbr_candidates::Int,
+ detected_event::Bool, tr_nplus1::Transition)
+ ind_edge = 0
+ bool_event = detected_event
+ for i = 1:nbr_candidates
+ edge = edge_candidates[i]
+ # Asynchronous edge detection: we fire it
+ if getfield(edge, :transitions) == nothing
+ return (i, detected_event)
+ end
+ # Synchronous detection
+ if !detected_event && tr_nplus1 != nothing
+ if (getfield(edge, :transitions)[1] == :ALL) ||
+ (tr_nplus1 in getfield(edge, :transitions))
+ ind_edge = i
+ bool_event = true
+ end
+ end
end
+ return (ind_edge, bool_event)
end
- end
- return (ind_edge, bool_event)
-end
-function next_state!(Snplus1::StateLHA, A::LHA,
- xnplus1::Vector{Int}, tnplus1::Float64, tr_nplus1::Transition,
- Sn::StateLHA, xn::Vector{Int}, p::Vector{Float64},
- edge_candidates::Vector{Edge}; verbose::Bool = false)
- # En fait d'apres observation de Cosmos, après qu'on ait lu la transition on devrait stop.
- detected_event::Bool = false
- turns = 0
- current_values = getfield(Snplus1, :values)
- current_time = getfield(Snplus1, :time)
- current_loc = getfield(Snplus1, :loc)
- Λ = getfield(A, :Λ)
- flow = getfield(A, :flow)
- map_edges = getfield(A, :map_edges)
+ function next_state!(Snplus1::StateLHA, A::LHA,
+ xnplus1::Vector{Int}, tnplus1::Float64, tr_nplus1::Transition,
+ Sn::StateLHA, xn::Vector{Int}, p::Vector{Float64},
+ edge_candidates::Vector{Edge}; verbose::Bool = false)
+ # En fait d'apres observation de Cosmos, après qu'on ait lu la transition on devrait stop.
+ detected_event::Bool = false
+ turns = 0
+ current_values = getfield(Snplus1, :values)
+ current_time = getfield(Snplus1, :time)
+ current_loc = getfield(Snplus1, :loc)
+ Λ = getfield(A, :Λ)
+ flow = getfield(A, :flow)
+ map_edges = getfield(A, :map_edges)
- if verbose
- println("##### Begin next_state!")
- @show xnplus1, tnplus1, tr_nplus1
- @show Sn
- end
- # In terms of values not reference, Snplus1 == Sn
- # First, we check the asynchronous transitions
- while true
- turns += 1
- #edge_candidates = empty!(edge_candidates)
- edges_from_current_loc = map_edges[current_loc]
- # Save all edges that satisfies transition predicate (asynchronous ones)
- nbr_candidates = _find_edge_candidates!(edge_candidates, edges_from_current_loc, Λ,
- current_time, current_values, xn, p, true)
- # Search the one we must chose, here the event is nothing because
- # we're not processing yet the next event
- ind_edge, detected_event = _get_edge_index(edge_candidates, nbr_candidates, detected_event, nothing)
- # Update the state with the chosen one (if it exists)
- # Should be xn here
- #first_round = false
- if ind_edge > 0
- current_loc = getfield(edge_candidates[ind_edge], :update_state!)(current_time, current_values, xn, p)
- else
- if verbose println("No edge fired") end
- break
- end
- if verbose
- @show turns
- @show edge_candidates
- @show ind_edge, detected_event, nbr_candidates
- @show current_loc
- @show current_time
- @show current_values
- if turns == 500
- @warn "We've reached 500 turns"
+ if verbose
+ println("##### Begin next_state!")
+ @show xnplus1, tnplus1, tr_nplus1
+ @show Sn
end
- end
- # For debug
- #=
- if turns > 100
- println("Number of turns in next_state! is suspicious")
- @show first_round, detected_event
- @show length(edge_candidates)
- @show tnplus1, tr_nplus1, xnplus1
- @show edge_candidates
- for edge in edge_candidates
+ # In terms of values not reference, Snplus1 == Sn
+ # First, we check the asynchronous transitions
+ while true
+ turns += 1
+ #edge_candidates = empty!(edge_candidates)
+ edges_from_current_loc = map_edges[current_loc]
+ # Save all edges that satisfies transition predicate (asynchronous ones)
+ nbr_candidates = _find_edge_candidates!(edge_candidates, edges_from_current_loc, Λ,
+ current_time, current_values, xn, p, true)
+ # Search the one we must chose, here the event is nothing because
+ # we're not processing yet the next event
+ ind_edge, detected_event = _get_edge_index(edge_candidates, nbr_candidates, detected_event, nothing)
+ # Update the state with the chosen one (if it exists)
+ # Should be xn here
+ #first_round = false
+ if ind_edge > 0
+ firing_edge = edge_candidates[ind_edge]
+ current_loc = getfield(Main, :update_state!)(firing_edge, current_time, current_values, xn, p)
+ else
+ if verbose println("No edge fired") end
+ break
+ end
+ if verbose
+ @show turns
+ @show edge_candidates
+ @show ind_edge, detected_event, nbr_candidates
+ @show current_loc
+ @show current_time
+ @show current_values
+ if turns == 500
+ @warn "We've reached 500 turns"
+ end
+ end
+ # For debug
+ #=
+ if turns > 100
+ println("Number of turns in next_state! is suspicious")
+ @show first_round, detected_event
+ @show length(edge_candidates)
+ @show tnplus1, tr_nplus1, xnplus1
+ @show edge_candidates
+ for edge in edge_candidates
@show getfield(edge, :check_constraints)(time_S, values_S, xn, p)
+ end
+ error("Unpredicted behavior automaton")
+ end
+ =#
end
- error("Unpredicted behavior automaton")
- end
- =#
- end
- if verbose
- println("Time flies with the flow...")
- end
- # Now time flies according to the flow
- for i in eachindex(current_values)
- @inbounds coeff_deriv = flow[current_loc][i]
- if coeff_deriv > 0
- @inbounds current_values[i] += coeff_deriv*(tnplus1 - current_time)
- end
- end
- current_time = tnplus1
- if verbose
- @show current_loc
- @show current_time
- @show current_values
- end
- # Now firing an edge according to the event
- while true
- turns += 1
- edges_from_current_loc = map_edges[current_loc]
- # Save all edges that satisfies transition predicate (synchronous ones)
- nbr_candidates = _find_edge_candidates!(edge_candidates, edges_from_current_loc, Λ,
- current_time, current_values, xnplus1, p, false)
- # Search the one we must chose
- ind_edge, detected_event = _get_edge_index(edge_candidates, nbr_candidates, detected_event, tr_nplus1)
- # Update the state with the chosen one (if it exists)
- if ind_edge > 0
- current_loc = getfield(edge_candidates[ind_edge], :update_state!)(current_time, current_values, xnplus1, p)
- end
- if ind_edge == 0 || detected_event
if verbose
- if detected_event
- println("Synchronized with $(tr_nplus1)")
+ println("Time flies with the flow...")
+ end
+ # Now time flies according to the flow
+ for i in eachindex(current_values)
+ @inbounds coeff_deriv = flow[current_loc][i]
+ if coeff_deriv > 0
+ @inbounds current_values[i] += coeff_deriv*(tnplus1 - current_time)
+ end
+ end
+ current_time = tnplus1
+ if verbose
+ @show current_loc
+ @show current_time
+ @show current_values
+ end
+ # Now firing an edge according to the event
+ while true
+ turns += 1
+ edges_from_current_loc = map_edges[current_loc]
+ # Save all edges that satisfies transition predicate (synchronous ones)
+ nbr_candidates = _find_edge_candidates!(edge_candidates, edges_from_current_loc, Λ,
+ current_time, current_values, xnplus1, p, false)
+ # Search the one we must chose
+ ind_edge, detected_event = _get_edge_index(edge_candidates, nbr_candidates, detected_event, tr_nplus1)
+ # Update the state with the chosen one (if it exists)
+ if ind_edge > 0
+ firing_edge = edge_candidates[ind_edge]
+ current_loc = getfield(Main, :update_state!)(firing_edge, current_time, current_values, xnplus1, p)
+ end
+ if ind_edge == 0 || detected_event
+ if verbose
+ if detected_event
+ println("Synchronized with $(tr_nplus1)")
+ @show turns
+ @show edge_candidates
+ @show ind_edge, detected_event, nbr_candidates
+ @show detected_event
+ @show current_loc
+ @show current_time
+ @show current_values
+ else
+ println("No edge fired")
+ end
+ end
+ break
+ end
+ if verbose
@show turns
@show edge_candidates
@show ind_edge, detected_event, nbr_candidates
@@ -256,44 +274,32 @@ function next_state!(Snplus1::StateLHA, A::LHA,
@show current_loc
@show current_time
@show current_values
- else
- println("No edge fired")
+ if turns == 500
+ @warn "We've reached 500 turns"
+ end
end
- end
- break
- end
- if verbose
- @show turns
- @show edge_candidates
- @show ind_edge, detected_event, nbr_candidates
- @show detected_event
- @show current_loc
- @show current_time
- @show current_values
- if turns == 500
- @warn "We've reached 500 turns"
- end
- end
- # For debug
- #=
- if turns > 100
- println("Number of turns in next_state! is suspicious")
- @show detected_event
- @show length(edge_candidates)
- @show tnplus1, tr_nplus1, xnplus1
- @show edge_candidates
- for edge in edge_candidates
+ # For debug
+ #=
+ if turns > 100
+ println("Number of turns in next_state! is suspicious")
+ @show detected_event
+ @show length(edge_candidates)
+ @show tnplus1, tr_nplus1, xnplus1
+ @show edge_candidates
+ for edge in edge_candidates
@show getfield(edge, :check_constraints)(time_S, values_S, x, p)
+ end
+ error("Unpredicted behavior automaton")
+ end
+ =#
+ end
+ setfield!(Snplus1, :loc, current_loc)
+ setfield!(Snplus1, :time, current_time)
+ if verbose
+ @show Snplus1
+ println("##### End next_state!")
end
- error("Unpredicted behavior automaton")
end
- =#
- end
- setfield!(Snplus1, :loc, current_loc)
- setfield!(Snplus1, :time, current_time)
- if verbose
- @show Snplus1
- println("##### End next_state!")
end
end
diff --git a/core/model.jl b/core/model.jl
index 3151968..5d2b022 100644
--- a/core/model.jl
+++ b/core/model.jl
@@ -23,308 +23,325 @@ function _update_values!(values::Vector{Vector{Int}}, times::Vector{Float64}, tr
@inbounds(transitions[idx] = tr_n)
end
-"""
- `simulate(m)`
+function generate_code_simulation(model_name::Symbol, f!::Symbol, isabsorbing::Symbol)
-Simulates a model. If `m::SynchronizedModel`, the simulation is synchronized with a
-Linear Hybrid Automaton.
-"""
-function simulate(m::ContinuousTimeModel; p::Union{Nothing,AbstractVector{Float64}} = nothing)
- p_sim = getfield(m, :p)
- if p != nothing
- p_sim = p
- end
- time_bound = getfield(m, :time_bound)
- buffer_size = getfield(m, :buffer_size)
- estim_min_states = getfield(m, :estim_min_states)
- # First alloc
- full_values = Vector{Vector{Int}}(undef, getfield(m, :dim_state))
- for i = eachindex(full_values) full_values[i] = zeros(Int, estim_min_states) end
- times = zeros(Float64, estim_min_states)
- transitions = Vector{Transition}(undef, estim_min_states)
- # Initial values
- for i = eachindex(full_values) full_values[i][1] = getfield(m, :x0)[i] end
- times[1] = getfield(m, :t0)
- transitions[1] = nothing
- # Values at time n
- n = 1
- xn = copy(getfield(m, :x0))
- tn = getfield(m, :t0)
- # at time n+1
- isabsorbing::Bool = getfield(m, :isabsorbing)(p_sim,xn)
- # If x0 is absorbing
- if isabsorbing
- _resize_trajectory!(full_values, times, transitions, 1)
- values = full_values[getfield(m, :_g_idx)]
- if isbounded(m)
- _finish_bounded_trajectory!(values, times, transitions, time_bound)
- end
- return Trajectory(m, values, times, transitions)
- end
- # Alloc of vectors where we stock n+1 values
- vec_x = zeros(Int, getfield(m, :dim_state))
- l_t = Float64[0.0]
- l_tr = Transition[nothing]
- # First we fill the allocated array
- for i = 2:estim_min_states
- getfield(m, :f!)(vec_x, l_t, l_tr, xn, tn, p_sim)
- tn = l_t[1]
- if tn > time_bound ||Â vec_x == xn
- isabsorbing = (vec_x == xn)
- break
- end
- n += 1
- copyto!(xn, vec_x)
- _update_values!(full_values, times, transitions, xn, tn, l_tr[1], i)
- end
- # If simulation ended before the estimation of states
- if n < estim_min_states
- _resize_trajectory!(full_values, times, transitions, n)
- values = full_values[getfield(m, :_g_idx)]
- if isbounded(m)
- _finish_bounded_trajectory!(values, times, transitions, time_bound)
- end
- return Trajectory(m, values, times, transitions)
- end
- # Otherwise, buffering system
- size_tmp = 0
- while !isabsorbing &&Â tn <= time_bound
- # Alloc buffer
- _resize_trajectory!(full_values, times, transitions, estim_min_states+size_tmp+buffer_size)
- i = 0
- while i < buffer_size
- i += 1
- getfield(m, :f!)(vec_x, l_t, l_tr, xn, tn, p_sim)
- tn = l_t[1]
- if tn > time_bound
- i -= 1
- break
+ return quote
+ import MarkovProcesses: simulate
+
+ """
+ `simulate(m)`
+
+ Simulates a model. If `m::SynchronizedModel`, the simulation is synchronized with a
+ Linear Hybrid Automaton.
+ """
+ function simulate(m::$(model_name); p::Union{Nothing,AbstractVector{Float64}} = nothing)
+ p_sim = getfield(m, :p)
+ if p != nothing
+ p_sim = p
end
- if vec_x == xn
- isabsorbing = true
- i -= 1
- break
+ time_bound = getfield(m, :time_bound)
+ buffer_size = getfield(m, :buffer_size)
+ estim_min_states = getfield(m, :estim_min_states)
+ # First alloc
+ full_values = Vector{Vector{Int}}(undef, getfield(m, :dim_state))
+ for i = eachindex(full_values) full_values[i] = zeros(Int, estim_min_states) end
+ times = zeros(Float64, estim_min_states)
+ transitions = Vector{Transition}(undef, estim_min_states)
+ # Initial values
+ for i = eachindex(full_values) full_values[i][1] = getfield(m, :x0)[i] end
+ times[1] = getfield(m, :t0)
+ transitions[1] = nothing
+ # Values at time n
+ n = 1
+ xn = copy(getfield(m, :x0))
+ tn = getfield(m, :t0)
+ # at time n+1
+ isabsorbing::Bool = $(isabsorbing)(p_sim,xn)
+ # If x0 is absorbing
+ if isabsorbing
+ MarkovProcesses._resize_trajectory!(full_values, times, transitions, 1)
+ values = full_values[getfield(m, :_g_idx)]
+ if isbounded(m)
+ MarkovProcesses._finish_bounded_trajectory!(values, times, transitions, time_bound)
+ end
+ return Trajectory(m, values, times, transitions)
end
- copyto!(xn, vec_x)
- _update_values!(full_values, times, transitions,
- xn, tn, l_tr[1], estim_min_states+size_tmp+i)
- end
- # If simulation ended before the end of buffer
- if i < buffer_size
- _resize_trajectory!(full_values, times, transitions, estim_min_states+size_tmp+i)
+ # Alloc of vectors where we stock n+1 values
+ vec_x = zeros(Int, getfield(m, :dim_state))
+ l_t = Float64[0.0]
+ l_tr = Transition[nothing]
+ # First we fill the allocated array
+ for i = 2:estim_min_states
+ $(f!)(vec_x, l_t, l_tr, xn, tn, p_sim)
+ tn = l_t[1]
+ if tn > time_bound ||Â vec_x == xn
+ isabsorbing = (vec_x == xn)
+ break
+ end
+ n += 1
+ copyto!(xn, vec_x)
+ MarkovProcesses._update_values!(full_values, times, transitions, xn, tn, l_tr[1], i)
+ end
+ # If simulation ended before the estimation of states
+ if n < estim_min_states
+ MarkovProcesses._resize_trajectory!(full_values, times, transitions, n)
+ values = full_values[getfield(m, :_g_idx)]
+ if isbounded(m)
+ MarkovProcesses._finish_bounded_trajectory!(values, times, transitions, time_bound)
+ end
+ return Trajectory(m, values, times, transitions)
+ end
+ # Otherwise, buffering system
+ size_tmp = 0
+ while !isabsorbing &&Â tn <= time_bound
+ # Alloc buffer
+ MarkovProcesses._resize_trajectory!(full_values, times, transitions, estim_min_states+size_tmp+buffer_size)
+ i = 0
+ while i < buffer_size
+ i += 1
+ $(f!)(vec_x, l_t, l_tr, xn, tn, p_sim)
+ tn = l_t[1]
+ if tn > time_bound
+ i -= 1
+ break
+ end
+ if vec_x == xn
+ isabsorbing = true
+ i -= 1
+ break
+ end
+ copyto!(xn, vec_x)
+ MarkovProcesses._update_values!(full_values, times, transitions,
+ xn, tn, l_tr[1], estim_min_states+size_tmp+i)
+ end
+ # If simulation ended before the end of buffer
+ if i < buffer_size
+ MarkovProcesses._resize_trajectory!(full_values, times, transitions, estim_min_states+size_tmp+i)
+ end
+ size_tmp += i
+ n += i
+ end
+ values = full_values[getfield(m, :_g_idx)]
+ if isbounded(m)
+ # Add last value: the convention is the last transition is nothing,
+ # the trajectory is bounded
+ MarkovProcesses._finish_bounded_trajectory!(values, times, transitions, time_bound)
+ end
+ return Trajectory(m, values, times, transitions)
end
- size_tmp += i
- n += i
end
- values = full_values[getfield(m, :_g_idx)]
- if isbounded(m)
- # Add last value: the convention is the last transition is nothing,
- # the trajectory is bounded
- _finish_bounded_trajectory!(values, times, transitions, time_bound)
- end
- return Trajectory(m, values, times, transitions)
end
-function simulate(product::SynchronizedModel; p::Union{Nothing,AbstractVector{Float64}} = nothing,
- verbose::Bool = false)
- m = getfield(product, :m)
- A = getfield(product, :automaton)
- p_sim = getfield(m, :p)
- if p != nothing
- p_sim = p
- end
- x0 = getfield(m, :x0)
- t0 = getfield(m, :t0)
- time_bound = getfield(m, :time_bound)
- S0 = init_state(A, x0, t0)
- buffer_size = getfield(m, :buffer_size)
- estim_min_states = getfield(m, :estim_min_states)
- edge_candidates = Vector{Edge}(undef, 2)
- # First alloc
- full_values = Vector{Vector{Int}}(undef, getfield(m, :dim_state))
- for i = eachindex(full_values) full_values[i] = zeros(Int, estim_min_states) end
- times = zeros(Float64, estim_min_states)
- transitions = Vector{Transition}(undef, estim_min_states)
- # Initial values
- for i = eachindex(full_values) full_values[i][1] = x0[i] end
- times[1] = t0
- transitions[1] = nothing
- # Values at time n
- n = 1
- xn = copy(x0)
- tn = copy(t0)
- Sn = copy(S0)
- next_state!(Sn, A, x0, t0, nothing, S0, x0, p_sim, edge_candidates; verbose = verbose)
- isabsorbing::Bool = getfield(m, :isabsorbing)(p_sim,xn)
- isacceptedLHA::Bool = isaccepted(Sn)
- # Alloc of vectors where we stock n+1 values
- vec_x = zeros(Int, getfield(m, :dim_state))
- l_t = Float64[0.0]
- l_tr = Transition[nothing]
- Snplus1 = deepcopy(Sn)
- # If x0 is absorbing
- if isabsorbing || isacceptedLHA
- _resize_trajectory!(full_values, times, transitions, 1)
- values = full_values[getfield(m, :_g_idx)]
- if isbounded(m)
- _finish_bounded_trajectory!(values, times, transitions, time_bound)
- end
- if isabsorbing && !isacceptedLHA
- next_state!(Snplus1, A, xn, time_bound, nothing, Sn, xn, p_sim, edge_candidates; verbose = verbose)
- copyto!(Sn, Snplus1)
- end
- return SynchronizedTrajectory(Sn, product, values, times, transitions)
- end
- # First we fill the allocated array
- for i = 2:estim_min_states
- getfield(m, :f!)(vec_x, l_t, l_tr, xn, tn, p_sim)
- if l_t[1] > time_bound ||Â vec_x == xn
- tn = l_t[1]
- isabsorbing = (vec_x == xn)
- break
- end
- n += 1
- next_state!(Snplus1, A, vec_x, l_t[1], l_tr[1], Sn, xn, p_sim, edge_candidates; verbose = verbose)
- copyto!(xn, vec_x)
- tn = l_t[1]
- tr_n = l_tr[1]
- _update_values!(full_values, times, transitions, xn, tn, tr_n, i)
- copyto!(Sn, Snplus1)
- isacceptedLHA = isaccepted(Snplus1)
- if isabsorbing ||Â isacceptedLHA
- break
- end
- end
- # If simulation ended before the estimation of states
- if n < estim_min_states
- _resize_trajectory!(full_values, times, transitions, n)
- values = full_values[getfield(m, :_g_idx)]
- if isbounded(m)
- _finish_bounded_trajectory!(values, times, transitions, time_bound)
- end
- if isabsorbing && !isacceptedLHA
- next_state!(Snplus1, A, xn, time_bound, nothing, Sn, xn, p_sim, edge_candidates; verbose = verbose)
- copyto!(Sn, Snplus1)
- end
- return SynchronizedTrajectory(Sn, product, values, times, transitions)
- end
- # Otherwise, buffering system
- size_tmp = 0
- while !isabsorbing &&Â tn <= time_bound && !isacceptedLHA
- # Alloc buffer
- _resize_trajectory!(full_values, times, transitions, estim_min_states+size_tmp+buffer_size)
- i = 0
- while i < buffer_size
- i += 1
- getfield(m, :f!)(vec_x, l_t, l_tr, xn, tn, p_sim)
- if l_t[1] > time_bound
+function generate_synchronized_simulation_code()
+
+ return quote
+ import MarkovProcesses: simulate
+
+ function simulate(product::SynchronizedModel; p::Union{Nothing,AbstractVector{Float64}} = nothing,
+ verbose::Bool = false)
+ m = getfield(product, :m)
+ A = getfield(product, :automaton)
+ p_sim = getfield(m, :p)
+ func_f! = getfield(m, :f!)
+ func_isabsorbing = getfield(m, :isabsorbing)
+ if p != nothing
+ p_sim = p
+ end
+ x0 = getfield(m, :x0)
+ t0 = getfield(m, :t0)
+ time_bound = getfield(m, :time_bound)
+ S0 = init_state(A, x0, t0)
+ buffer_size = getfield(m, :buffer_size)
+ estim_min_states = getfield(m, :estim_min_states)
+ edge_candidates = Vector{Edge}(undef, 2)
+ # First alloc
+ full_values = Vector{Vector{Int}}(undef, getfield(m, :dim_state))
+ for i = eachindex(full_values) full_values[i] = zeros(Int, estim_min_states) end
+ times = zeros(Float64, estim_min_states)
+ transitions = Vector{Transition}(undef, estim_min_states)
+ # Initial values
+ for i = eachindex(full_values) full_values[i][1] = x0[i] end
+ times[1] = t0
+ transitions[1] = nothing
+ # Values at time n
+ n = 1
+ xn = copy(x0)
+ tn = copy(t0)
+ Sn = copy(S0)
+ next_state!(Sn, A, x0, t0, nothing, S0, x0, p_sim, edge_candidates; verbose = verbose)
+ isabsorbing::Bool = func_isabsorbing(p_sim,xn)
+ isacceptedLHA::Bool = isaccepted(Sn)
+ # Alloc of vectors where we stock n+1 values
+ vec_x = zeros(Int, getfield(m, :dim_state))
+ l_t = Float64[0.0]
+ l_tr = Transition[nothing]
+ Snplus1 = deepcopy(Sn)
+ # If x0 is absorbing
+ if isabsorbing || isacceptedLHA
+ MarkovProcesses._resize_trajectory!(full_values, times, transitions, 1)
+ values = full_values[getfield(m, :_g_idx)]
+ if isbounded(m)
+ MarkovProcesses._finish_bounded_trajectory!(values, times, transitions, time_bound)
+ end
+ if isabsorbing && !isacceptedLHA
+ next_state!(Snplus1, A, xn, time_bound, nothing, Sn, xn, p_sim, edge_candidates; verbose = verbose)
+ copyto!(Sn, Snplus1)
+ end
+ return SynchronizedTrajectory(Sn, product, values, times, transitions)
+ end
+ # First we fill the allocated array
+ for i = 2:estim_min_states
+ func_f!(vec_x, l_t, l_tr, xn, tn, p_sim)
+ if l_t[1] > time_bound ||Â vec_x == xn
+ tn = l_t[1]
+ isabsorbing = (vec_x == xn)
+ break
+ end
+ n += 1
+ next_state!(Snplus1, A, vec_x, l_t[1], l_tr[1], Sn, xn, p_sim, edge_candidates; verbose = verbose)
+ copyto!(xn, vec_x)
tn = l_t[1]
- i -= 1
- break
+ tr_n = l_tr[1]
+ MarkovProcesses._update_values!(full_values, times, transitions, xn, tn, tr_n, i)
+ copyto!(Sn, Snplus1)
+ isacceptedLHA = isaccepted(Snplus1)
+ if isabsorbing ||Â isacceptedLHA
+ break
+ end
end
- if vec_x == xn
- isabsorbing = true
- i -= 1
- break
+ # If simulation ended before the estimation of states
+ if n < estim_min_states
+ MarkovProcesses._resize_trajectory!(full_values, times, transitions, n)
+ values = full_values[getfield(m, :_g_idx)]
+ if isbounded(m)
+ MarkovProcesses._finish_bounded_trajectory!(values, times, transitions, time_bound)
+ end
+ if isabsorbing && !isacceptedLHA
+ next_state!(Snplus1, A, xn, time_bound, nothing, Sn, xn, p_sim, edge_candidates; verbose = verbose)
+ copyto!(Sn, Snplus1)
+ end
+ return SynchronizedTrajectory(Sn, product, values, times, transitions)
end
- next_state!(Snplus1, A, vec_x, l_t[1], l_tr[1], Sn, xn, p_sim, edge_candidates; verbose = verbose)
- copyto!(xn, vec_x)
- tn = l_t[1]
- tr_n = l_tr[1]
- _update_values!(full_values, times, transitions,
- xn, tn, tr_n, estim_min_states+size_tmp+i)
- copyto!(Sn, Snplus1)
- isacceptedLHA = isaccepted(Snplus1)
- if isabsorbing || isacceptedLHA
- break
+ # Otherwise, buffering system
+ size_tmp = 0
+ while !isabsorbing &&Â tn <= time_bound && !isacceptedLHA
+ # Alloc buffer
+ MarkovProcesses._resize_trajectory!(full_values, times, transitions, estim_min_states+size_tmp+buffer_size)
+ i = 0
+ while i < buffer_size
+ i += 1
+ func_f!(vec_x, l_t, l_tr, xn, tn, p_sim)
+ if l_t[1] > time_bound
+ tn = l_t[1]
+ i -= 1
+ break
+ end
+ if vec_x == xn
+ isabsorbing = true
+ i -= 1
+ break
+ end
+ next_state!(Snplus1, A, vec_x, l_t[1], l_tr[1], Sn, xn, p_sim, edge_candidates; verbose = verbose)
+ copyto!(xn, vec_x)
+ tn = l_t[1]
+ tr_n = l_tr[1]
+ MarkovProcesses._update_values!(full_values, times, transitions,
+ xn, tn, tr_n, estim_min_states+size_tmp+i)
+ copyto!(Sn, Snplus1)
+ isacceptedLHA = isaccepted(Snplus1)
+ if isabsorbing || isacceptedLHA
+ break
+ end
+ end
+ # If simulation ended before the end of buffer
+ if i < buffer_size
+ MarkovProcesses._resize_trajectory!(full_values, times, transitions, estim_min_states+size_tmp+i)
+ end
+ size_tmp += i
+ n += i
end
+ values = full_values[getfield(m, :_g_idx)]
+ if isbounded(m) && !isaccepted(Sn)
+ # Add last value: the convention is the last transition is nothing,
+ # the trajectory is bounded
+ MarkovProcesses._finish_bounded_trajectory!(values, times, transitions, time_bound)
+ end
+ if isabsorbing && !isacceptedLHA
+ next_state!(Snplus1, A, xn, time_bound, nothing, Sn, xn, p_sim, edge_candidates; verbose = verbose)
+ copyto!(Sn, Snplus1)
+ end
+ return SynchronizedTrajectory(Sn, product, values, times, transitions)
end
- # If simulation ended before the end of buffer
- if i < buffer_size
- _resize_trajectory!(full_values, times, transitions, estim_min_states+size_tmp+i)
- end
- size_tmp += i
- n += i
- end
- values = full_values[getfield(m, :_g_idx)]
- if isbounded(m) && !isaccepted(Sn)
- # Add last value: the convention is the last transition is nothing,
- # the trajectory is bounded
- _finish_bounded_trajectory!(values, times, transitions, time_bound)
- end
- if isabsorbing && !isacceptedLHA
- next_state!(Snplus1, A, xn, time_bound, nothing, Sn, xn, p_sim, edge_candidates; verbose = verbose)
- copyto!(Sn, Snplus1)
- end
- return SynchronizedTrajectory(Sn, product, values, times, transitions)
-end
-"""
- `volatile_simulate(sm::SynchronizedModel; p, verbose)`
+ """
+ `volatile_simulate(sm::SynchronizedModel; p, verbose)`
-Simulates a model synchronized with an automaton but does not store the values of the simulation
-in order to improve performance.
-It returns the last state of the simulation `S::StateLHA` not a trajectory `σ::SynchronizedTrajectory`.
-"""
-function volatile_simulate(product::SynchronizedModel;
- p::Union{Nothing,AbstractVector{Float64}} = nothing, verbose::Bool = false)
- m = product.m
- A = product.automaton
- p_sim = getfield(m, :p)
- if p != nothing
- p_sim = p
- end
- x0 = getfield(m, :x0)
- t0 = getfield(m, :t0)
- time_bound = getfield(m, :time_bound)
- S0 = init_state(A, x0, t0)
- edge_candidates = Vector{Edge}(undef, 2)
- # Values at time n
- n = 1
- xn = copy(x0)
- tn = copy(t0)
- Sn = copy(S0)
- next_state!(Sn, A, x0, t0, nothing, S0, x0, p_sim, edge_candidates; verbose = verbose)
- isabsorbing::Bool = getfield(m, :isabsorbing)(p_sim,xn)
- isacceptedLHA::Bool = isaccepted(Sn)
- # Alloc of vectors where we stock n+1 values
- vec_x = zeros(Int, getfield(m, :dim_state))
- l_t = Float64[0.0]
- l_tr = Transition[nothing]
- Snplus1 = deepcopy(Sn)
- # If x0 is absorbing
- if isabsorbing || isacceptedLHA
- if !isacceptedLHA
- next_state!(Snplus1, A, xn, time_bound, nothing, Sn, xn, p_sim, edge_candidates; verbose = verbose)
- copyto!(Sn, Snplus1)
- end
- return Sn
- end
- while !isabsorbing &&Â tn <= time_bound && !isacceptedLHA
- getfield(m, :f!)(vec_x, l_t, l_tr, xn, tn, p_sim)
- if l_t[1] > time_bound
- tn = l_t[1]
- break
- end
- if vec_x == xn
- isabsorbing = true
- break
+ Simulates a model synchronized with an automaton but does not store the values of the simulation
+ in order to improve performance.
+ It returns the last state of the simulation `S::StateLHA` not a trajectory `σ::SynchronizedTrajectory`.
+ """
+ function volatile_simulate(product::SynchronizedModel;
+ p::Union{Nothing,AbstractVector{Float64}} = nothing, verbose::Bool = false)
+ m = product.m
+ A = product.automaton
+ p_sim = getfield(m, :p)
+ func_f! = getfield(m, :f!)
+ func_isabsorbing = getfield(m, :isabsorbing)
+ if p != nothing
+ p_sim = p
+ end
+ x0 = getfield(m, :x0)
+ t0 = getfield(m, :t0)
+ time_bound = getfield(m, :time_bound)
+ S0 = init_state(A, x0, t0)
+ edge_candidates = Vector{Edge}(undef, 2)
+ # Values at time n
+ n = 1
+ xn = copy(x0)
+ tn = copy(t0)
+ Sn = copy(S0)
+ next_state!(Sn, A, x0, t0, nothing, S0, x0, p_sim, edge_candidates; verbose = verbose)
+ isabsorbing::Bool = func_isabsorbing(p_sim,xn)
+ isacceptedLHA::Bool = isaccepted(Sn)
+ # Alloc of vectors where we stock n+1 values
+ vec_x = zeros(Int, getfield(m, :dim_state))
+ l_t = Float64[0.0]
+ l_tr = Transition[nothing]
+ Snplus1 = deepcopy(Sn)
+ # If x0 is absorbing
+ if isabsorbing || isacceptedLHA
+ if !isacceptedLHA
+ next_state!(Snplus1, A, xn, time_bound, nothing, Sn, xn, p_sim, edge_candidates; verbose = verbose)
+ copyto!(Sn, Snplus1)
+ end
+ return Sn
+ end
+ while !isabsorbing &&Â tn <= time_bound && !isacceptedLHA
+ func_f!(vec_x, l_t, l_tr, xn, tn, p_sim)
+ if l_t[1] > time_bound
+ tn = l_t[1]
+ break
+ end
+ if vec_x == xn
+ isabsorbing = true
+ break
+ end
+ next_state!(Snplus1, A, vec_x, l_t[1], l_tr[1], Sn, xn, p_sim, edge_candidates; verbose = verbose)
+ copyto!(xn, vec_x)
+ tn = l_t[1]
+ tr_n = l_tr[1]
+ copyto!(Sn, Snplus1)
+ isacceptedLHA = isaccepted(Snplus1)
+ n += 1
+ end
+ if isabsorbing && !isacceptedLHA
+ next_state!(Snplus1, A, xn, time_bound, nothing, Sn, xn, p_sim, edge_candidates; verbose = verbose)
+ copyto!(Sn, Snplus1)
+ end
+ return Sn
end
- next_state!(Snplus1, A, vec_x, l_t[1], l_tr[1], Sn, xn, p_sim, edge_candidates; verbose = verbose)
- copyto!(xn, vec_x)
- tn = l_t[1]
- tr_n = l_tr[1]
- copyto!(Sn, Snplus1)
- isacceptedLHA = isaccepted(Snplus1)
- n += 1
- end
- if isabsorbing && !isacceptedLHA
- next_state!(Snplus1, A, xn, time_bound, nothing, Sn, xn, p_sim, edge_candidates; verbose = verbose)
- copyto!(Sn, Snplus1)
end
- return Sn
end
-
"""
`simulate(pm::ParametricModel, p_prior::AbstractVector{Float64})
@@ -425,7 +442,7 @@ function distribute_prob_accept_lha(sm::SynchronizedModel, nbr_sim::Int)
end
function Base.show(io::IO, m::ContinuousTimeModel)
- print(io, "$(m.name) model (ContinuousTimeModel)\n")
+ print(io, "$(typeof(m)) <: ContinuousTimeModel model\n")
print(io, "- variables :\n")
for (var, idx) in m.map_var_idx
print(io, "* $var (index = $idx in state space)\n")
@@ -454,7 +471,7 @@ function check_consistency(m::ContinuousTimeModel)
@assert length(m.g) <= m.dim_state
@assert length(m._g_idx) == length(m.g)
@assert m.buffer_size >= 0
- @assert typeof(m.isabsorbing(m.p, m.x0)) == Bool
+ @assert typeof(getfield(Main, m.isabsorbing)(m.p, m.x0)) == Bool
return true
end
diff --git a/core/network_model.jl b/core/network_model.jl
index b55bac4..9f3a552 100644
--- a/core/network_model.jl
+++ b/core/network_model.jl
@@ -72,7 +72,6 @@ fill_params!(dict_params::Dict{ParameterModel,Int}, l_dim_params::Vector{Int},
propensity::Real, list_species::Vector) = nothing
macro network_model(expr_network,expr_name...)
- model_name = isempty(expr_name) ? "Unnamed macro generated" : expr_name[1]
transitions = Transition[]
dict_species = Dict{VariableModel,Int}()
dict_params = Dict{ParameterModel,Int}()
@@ -142,12 +141,18 @@ macro network_model(expr_network,expr_name...)
=#
end
dim_params = l_dim_params[1]
- # Let's write some lines that creates the function f! (step of a simulation) for this biochemical network
- nbr_rand = rand(1:1000)
+
+ # Creation of names variables
+ model_name = isempty(expr_name) ? "Network" : expr_name[1]
+ model_name = Symbol(replace(model_name, ' ' => '_') * "Model")
+ id = Dates.format(Dates.now(), "YmHMs")
nbr_reactions = length(list_expr_reactions)
- basename_func = "$(replace(model_name, ' '=>'_'))_$(nbr_rand)"
+ basename_func = "$(model_name)_$(id)"
basename_func = replace(basename_func, '-'=>'_')
- expr_model_f! = "function $(basename_func)_f!(xnplus1::Vector{Int}, l_t::Vector{Float64}, l_tr::Vector{Transition}, xn::Vector{Int}, tn::Float64, p::Vector{Float64})\n\t"
+
+ # Writing of f!
+ symbol_func_f! = Symbol("$(basename_func)_f!")
+ str_expr_model_f! = "function $(symbol_func_f!)(xnplus1::Vector{Int}, l_t::Vector{Float64}, l_tr::Vector{Transition}, xn::Vector{Int}, tn::Float64, p::Vector{Float64})\n\t"
# Computation of nu and propensity functions in f!
str_l_a = "l_a = ("
str_test_isabsorbing = "@inbounds("
@@ -156,7 +161,7 @@ macro network_model(expr_network,expr_name...)
local isreaction = @capture(expr_reaction, TR_: (reactants_ => products_, propensity_))
# Writing of propensities function
str_propensity = get_str_propensity(propensity, dict_species, dict_params)
- expr_model_f! *= "@inbounds a$(i) = " * str_propensity * "\n\t"
+ str_expr_model_f! *= "@inbounds a$(i) = " * str_propensity * "\n\t"
# Anticipating the write of the function isabsorbing
str_test_isabsorbing *= str_propensity * "+"
# Update the nu of the i-th reaction
@@ -187,52 +192,65 @@ macro network_model(expr_network,expr_name...)
nu[dict_species[sym_species]] += mult
end
end
- expr_model_f! *= "nu_$i = $(Tuple(nu))\n\t"
+ str_expr_model_f! *= "nu_$i = $(Tuple(nu))\n\t"
# Anticipating the line l_a = (..)
str_l_a *= "a$(i), "
end
str_test_isabsorbing = str_test_isabsorbing[1:(end-1)] * ")"
str_l_a = str_l_a[1:(end-2)] * ")\n\t"
- expr_model_f! *= str_l_a
- expr_model_f! *= "asum = sum(l_a)\n\t"
- expr_model_f! *= "if asum == 0.0\n\t\t"
- expr_model_f! *= "copyto!(xnplus1, xn)\n\t\t"
- expr_model_f! *= "return nothing\n\t"
- expr_model_f! *= "end\n\t"
+ str_expr_model_f! *= str_l_a
+ str_expr_model_f! *= "asum = sum(l_a)\n\t"
+ str_expr_model_f! *= "if asum == 0.0\n\t\t"
+ str_expr_model_f! *= "copyto!(xnplus1, xn)\n\t\t"
+ str_expr_model_f! *= "return nothing\n\t"
+ str_expr_model_f! *= "end\n\t"
# Computation of array of transitions
- expr_model_f! *= "l_nu = (" * reduce(*, ["nu_$i, " for i = 1:nbr_reactions])[1:(end-2)] * ")\n\t"
- expr_model_f! *= "l_sym_R = $(Tuple(transitions))\n\t"
+ str_expr_model_f! *= "l_nu = (" * reduce(*, ["nu_$i, " for i = 1:nbr_reactions])[1:(end-2)] * ")\n\t"
+ str_expr_model_f! *= "l_sym_R = $(Tuple(transitions))\n\t"
# Simulation of the reaction
- expr_model_f! *= "u1 = rand()\n\t"
- expr_model_f! *= "u2 = rand()\n\t"
- expr_model_f! *= "tau = - log(u1) / asum\n\t"
- expr_model_f! *= "b_inf = 0.0\n\t"
- expr_model_f! *= "b_sup = a1\n\t"
- expr_model_f! *= "reaction = 0\n\n\t"
- expr_model_f! *= "for i = 1:$(nbr_reactions)\n\t\t"
- expr_model_f! *= "if b_inf < asum*u2 < b_sup\n\t\t\t"
- expr_model_f! *= "reaction = i\n\t\t\t"
- expr_model_f! *= "break\n\t\t"
- expr_model_f! *= "end\n\t\t"
- expr_model_f! *= "@inbounds b_inf += l_a[i]\n\t\t"
- expr_model_f! *= "@inbounds b_sup += l_a[i+1]\n\t"
- expr_model_f! *= "end\n\t"
- expr_model_f! *= "nu = l_nu[reaction]\n\t"
- expr_model_f! *= "for i = 1:$(dim_state)\n\t\t"
- expr_model_f! *= "@inbounds xnplus1[i] = xn[i]+nu[i]\n\t"
- expr_model_f! *= "end\n\t"
- expr_model_f! *= "@inbounds l_t[1] = tn + tau\n\t"
- expr_model_f! *= "@inbounds l_tr[1] = l_sym_R[reaction]\n"
- expr_model_f! *= "end\n"
- expr_model_isabsorbing = "isabsorbing_$(basename_func)(p::Vector{Float64},xn::Vector{Int}) = $(str_test_isabsorbing) === 0.0"
- @everywhere eval(Meta.parse($expr_model_f!))
- @everywhere eval(Meta.parse($expr_model_isabsorbing))
- symbol_func_f! = Symbol("$(basename_func)_f!")
+ str_expr_model_f! *= "u1 = rand()\n\t"
+ str_expr_model_f! *= "u2 = rand()\n\t"
+ str_expr_model_f! *= "tau = - log(u1) / asum\n\t"
+ str_expr_model_f! *= "b_inf = 0.0\n\t"
+ str_expr_model_f! *= "b_sup = a1\n\t"
+ str_expr_model_f! *= "reaction = 0\n\n\t"
+ str_expr_model_f! *= "for i = 1:$(nbr_reactions)\n\t\t"
+ str_expr_model_f! *= "if b_inf < asum*u2 < b_sup\n\t\t\t"
+ str_expr_model_f! *= "reaction = i\n\t\t\t"
+ str_expr_model_f! *= "break\n\t\t"
+ str_expr_model_f! *= "end\n\t\t"
+ str_expr_model_f! *= "@inbounds b_inf += l_a[i]\n\t\t"
+ str_expr_model_f! *= "@inbounds b_sup += l_a[i+1]\n\t"
+ str_expr_model_f! *= "end\n\t"
+ str_expr_model_f! *= "nu = l_nu[reaction]\n\t"
+ str_expr_model_f! *= "for i = 1:$(dim_state)\n\t\t"
+ str_expr_model_f! *= "@inbounds xnplus1[i] = xn[i]+nu[i]\n\t"
+ str_expr_model_f! *= "end\n\t"
+ str_expr_model_f! *= "@inbounds l_t[1] = tn + tau\n\t"
+ str_expr_model_f! *= "@inbounds l_tr[1] = l_sym_R[reaction]\n"
+ str_expr_model_f! *= "end\n"
+
+ # Writing of isabsorbing
symbol_func_isabsorbing = Symbol("isabsorbing_$(basename_func)")
+ str_expr_model_isabsorbing = "$(symbol_func_isabsorbing)(p::Vector{Float64},xn::Vector{Int}) = $(str_test_isabsorbing) === 0.0"
+
+ # Creation of code
+ expr_model_f! = Meta.parse(str_expr_model_f!)
+ expr_model_isabsorbing = Meta.parse(str_expr_model_isabsorbing)
+
map_idx_var_model = Dict(value => key for (key, value) in dict_species)
model_g = [map_idx_var_model[i] for i = 1:length(list_species)]
- return :(ContinuousTimeModel($dim_state, $dim_params, $dict_species, $dict_params, $transitions,
- $(zeros(dim_params)), $(zeros(Int, dim_state)), 0.0, $(getfield(Main, symbol_func_f!)), $(getfield(Main, symbol_func_isabsorbing));
- g = $model_g, name=$model_name))
+
+ return quote
+ @everywhere @eval $(MarkovProcesses.generate_code_model_type_def(model_name))
+ @everywhere @eval $(MarkovProcesses.generate_code_model_type_constructor(model_name))
+ @everywhere @eval $(MarkovProcesses.generate_code_simulation(model_name, symbol_func_f!, symbol_func_isabsorbing))
+ @everywhere @eval $expr_model_f!
+ @everywhere @eval $expr_model_isabsorbing
+
+ getfield(Main, $(Meta.quot(model_name)))($dim_state, $dim_params, $dict_species, $dict_params, $transitions,
+ $(zeros(dim_params)), $(zeros(Int, dim_state)), 0.0,
+ $(Meta.quot(symbol_func_f!)), $(Meta.quot(symbol_func_isabsorbing)); g = $model_g)
+ end
end
diff --git a/core/trajectory.jl b/core/trajectory.jl
index f2ec60c..8a55b46 100644
--- a/core/trajectory.jl
+++ b/core/trajectory.jl
@@ -197,8 +197,7 @@ function check_consistency(σ::AbstractTrajectory)
end
function Base.show(io::IO, σ::Trajectory)
- print(io, "Trajectory\n")
- print(io, "- Model name: $((σ.m).name) \n")
+ print(io, "Trajectory of $(σ.m)\n")
print(io, "- Variable trajectories:\n")
for obs_var in σ.m.g
print(io, "* $obs_var: $(σ[obs_var])\n")
diff --git a/core/utils.jl b/core/utils.jl
index c2322e7..83bf39e 100644
--- a/core/utils.jl
+++ b/core/utils.jl
@@ -23,7 +23,7 @@ function cosmos_get_values(name_file::String)
return dict_values
end
-load_model(name_model::String) = include("$(get_module_path())/models/$(name_model).jl")
-load_automaton(automaton::String) = include("$(get_module_path())/automata/$(automaton).jl")
-load_plots() = include(get_module_path() * "/core/plots.jl")
+load_model(name_model::String) = Base.MainInclude.include("$(get_module_path())/models/$(name_model).jl")
+load_automaton(automaton::String) = Base.MainInclude.include("$(get_module_path())/automata/$(automaton).jl")
+load_plots() = Base.MainInclude.include(get_module_path() * "/core/plots.jl")
diff --git a/models/ER.jl b/models/ER.jl
index fffe144..bd5d23c 100644
--- a/models/ER.jl
+++ b/models/ER.jl
@@ -50,7 +50,12 @@ end
@inbounds(p[1]*xn[1]*xn[2] + (p[2]+p[3])*xn[3] === 0.0)
g_ER = [:P]
-ER = ContinuousTimeModel(d,k,dict_var_ER,dict_p_ER,l_tr_ER,p_ER,x0_ER,t0_ER, getfield(Main, :ER_f!), getfield(Main, :isabsorbing_ER); g=g_ER,name="ER SSA pkg")
+@everywhere @eval $(MarkovProcesses.generate_code_model_type_def(:ERModel))
+@everywhere @eval $(MarkovProcesses.generate_code_model_type_constructor(:ERModel))
+@everywhere @eval $(MarkovProcesses.generate_code_simulation(:ERModel, :ER_f!, :isabsorbing_ER))
+
+ER = ERModel(d, k, dict_var_ER, dict_p_ER, l_tr_ER, p_ER, x0_ER, t0_ER,
+ :ER_f!, :isabsorbing_ER; g=g_ER)
function create_ER(new_p::Vector{Float64})
ER_new = deepcopy(ER)
@@ -59,5 +64,3 @@ function create_ER(new_p::Vector{Float64})
return ER_new
end
-export ER, create_ER
-
diff --git a/models/SIR.jl b/models/SIR.jl
index 16ac093..240ab49 100644
--- a/models/SIR.jl
+++ b/models/SIR.jl
@@ -17,7 +17,6 @@ t0_SIR = 0.0
copyto!(xnplus1, xn)
return nothing
end
- # column-major order
nu_1 = (-1, 1, 0)
nu_2 = (0, -1, 1)
l_nu = (nu_1, nu_2)
@@ -48,7 +47,12 @@ end
@everywhere isabsorbing_SIR(p::Vector{Float64}, xn::Vector{Int}) = (p[1]*xn[1]*xn[2] + p[2]*xn[2]) === 0.0
g_SIR = [:I]
-SIR = ContinuousTimeModel(d,k,dict_var,dict_p,l_tr_SIR,p_SIR,x0_SIR,t0_SIR, getfield(Main, :SIR_f!), getfield(Main, :isabsorbing_SIR); g=g_SIR, name="SIR SSA pkg")
+@everywhere @eval $(MarkovProcesses.generate_code_model_type_def(:SIRModel))
+@everywhere @eval $(MarkovProcesses.generate_code_model_type_constructor(:SIRModel))
+@everywhere @eval $(MarkovProcesses.generate_code_simulation(:SIRModel, :SIR_f!, :isabsorbing_SIR))
+
+SIR = SIRModel(d, k, dict_var, dict_p, l_tr_SIR, p_SIR, x0_SIR, t0_SIR,
+ :SIR_f!, :isabsorbing_SIR; g=g_SIR)
function create_SIR(new_p::Vector{Float64})
SIR_new = deepcopy(SIR)
@@ -57,5 +61,3 @@ function create_SIR(new_p::Vector{Float64})
return SIR_new
end
-export SIR, create_SIR
-
diff --git a/models/SIR_tauleap.jl b/models/SIR_tauleap.jl
index 85dbdc1..01af738 100644
--- a/models/SIR_tauleap.jl
+++ b/models/SIR_tauleap.jl
@@ -9,7 +9,7 @@ l_tr_SIR_tauleap = [:U]
p_SIR_tauleap = [0.0012, 0.05, 5.0]
x0_SIR_tauleap = [95, 5, 0]
t0_SIR_tauleap = 0.0
-@everywhere function SIR_tauleap_f!(xnplus1::Vector{Int}, l_t::Vector{Float64}, l_tr::Vector{Transition},
+@everywhere function SIRTauleap_f!(xnplus1::Vector{Int}, l_t::Vector{Float64}, l_tr::Vector{Transition},
xn::Vector{Int}, tn::Float64, p::Vector{Float64})
@inbounds tau = p[3]
@inbounds a1 = p[1] * xn[1] * xn[2]
@@ -31,12 +31,17 @@ t0_SIR_tauleap = 0.0
l_t[1] = tn + tau
l_tr[1] = :U
end
-@everywhere isabsorbing_SIR_tauleap(p::Vector{Float64}, xn::Vector{Int}) = (p[1]*xn[1]*xn[2] + p[2]*xn[2]) === 0.0
+@everywhere isabsorbing_SIRTauleap(p::Vector{Float64}, xn::Vector{Int}) = (p[1]*xn[1]*xn[2] + p[2]*xn[2]) === 0.0
g_SIR_tauleap = [:I]
-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,
- getfield(Main, :SIR_tauleap_f!), getfield(Main, :isabsorbing_SIR_tauleap); g=g_SIR_tauleap, name="SIR tauleap pkg")
+@everywhere @eval $(MarkovProcesses.generate_code_model_type_def(:SIRTauleapModel))
+@everywhere @eval $(MarkovProcesses.generate_code_model_type_constructor(:SIRTauleapModel))
+@everywhere @eval $(MarkovProcesses.generate_code_simulation(:SIRTauleapModel, :SIRTauleap_f!, :isabsorbing_SIRTauleap))
+
+SIR_tauleap = SIRTauleapModel(d, k, dict_var_SIR_tauleap, dict_p_SIR_tauleap, l_tr_SIR_tauleap,
+ p_SIR_tauleap, x0_SIR_tauleap, t0_SIR_tauleap,
+ :SIRTauleap_f!, :isabsorbing_SIRTauleap; g=g_SIR_tauleap)
+
function create_SIR_tauleap(new_p::Vector{Float64})
SIR_tauleap_new = deepcopy(SIR_tauleap)
@assert length(SIR_tauleap_new.p) == length(new_p)
@@ -44,5 +49,3 @@ function create_SIR_tauleap(new_p::Vector{Float64})
return SIR_tauleap_new
end
-export SIR_tauleap, create_SIR_tauleap
-
diff --git a/models/doping_3way_oscillator.jl b/models/doping_3way_oscillator.jl
index 5e4be8e..a338897 100644
--- a/models/doping_3way_oscillator.jl
+++ b/models/doping_3way_oscillator.jl
@@ -12,5 +12,3 @@ set_x0!(doping_3way_oscillator, [:A,:B,:C,:DA,:DB,:DC], [333,333,333,10,10,10])
set_param!(doping_3way_oscillator, [:rA,:rB,:rC], [1.0, 1.0, 1.0])
set_time_bound!(doping_3way_oscillator, 0.2)
-export doping_3way_oscillator
-
diff --git a/models/intracellular_viral_infection.jl b/models/intracellular_viral_infection.jl
index 9de65fa..4984277 100644
--- a/models/intracellular_viral_infection.jl
+++ b/models/intracellular_viral_infection.jl
@@ -12,5 +12,3 @@ set_x0!(intracellular_viral_infection, [:N, :A, :G, :T, :S, :V], [10000, 10000,
set_param!(intracellular_viral_infection, [:cn, :ca, :k1, :k2, :k3, :k4, :k5, :k6], [1.0, 1.0, 1.0, 0.025, 100.0, 0.25, 0.2, 7.5E-6])
set_time_bound!(intracellular_viral_infection, 200.0)
-export intracellular_viral_infection
-
diff --git a/models/poisson.jl b/models/poisson.jl
index 02de095..604e033 100644
--- a/models/poisson.jl
+++ b/models/poisson.jl
@@ -10,7 +10,7 @@ l_tr_poisson = [:R]
p_poisson = [5.0]
x0_poisson = [0]
t0_poisson = 0.0
-@everywhere function poisson_f!(xnplus1::Vector{Int}, l_t::Vector{Float64}, l_tr::Vector{Transition},
+@everywhere function Poisson_f!(xnplus1::Vector{Int}, l_t::Vector{Float64}, l_tr::Vector{Transition},
xn::Vector{Int}, tn::Float64, p::Vector{Float64})
u1 = rand()
@@ -19,13 +19,17 @@ t0_poisson = 0.0
l_t[1] = tn + tau
l_tr[1] = :R
end
-@everywhere isabsorbing_poisson(p::Vector{Float64}, xn::Vector{Int}) = p[1] === 0.0
+@everywhere isabsorbing_Poisson(p::Vector{Float64}, xn::Vector{Int}) = p[1] === 0.0
g_poisson = [:N]
-poisson = ContinuousTimeModel(d,k,dict_var_poisson,dict_p_poisson,l_tr_poisson,
- p_poisson,x0_poisson,t0_poisson,
- getfield(Main, :poisson_f!), getfield(Main, :isabsorbing_poisson);
- g=g_poisson, time_bound=1.0, name="Poisson process pkg")
+@everywhere @eval $(MarkovProcesses.generate_code_model_type_def(:PoissonModel))
+@everywhere @eval $(MarkovProcesses.generate_code_model_type_constructor(:PoissonModel))
+@everywhere @eval $(MarkovProcesses.generate_code_simulation(:PoissonModel, :Poisson_f!, :isabsorbing_Poisson))
+
+poisson = PoissonModel(d, k, dict_var_poisson, dict_p_poisson, l_tr_poisson,
+ p_poisson, x0_poisson, t0_poisson,
+ :Poisson_f!, :isabsorbing_Poisson; g=g_poisson, time_bound=1.0)
+
function create_poisson(new_p::Vector{Float64})
poisson_new = deepcopy(poisson)
@assert length(poisson_new.p) == length(new_p)
@@ -33,5 +37,3 @@ function create_poisson(new_p::Vector{Float64})
return poisson_new
end
-export poisson, create_poisson
-
diff --git a/models/repressilator.jl b/models/repressilator.jl
index 371c267..b5e907f 100644
--- a/models/repressilator.jl
+++ b/models/repressilator.jl
@@ -18,8 +18,6 @@ set_observed_var!(repressilator, [:mRNA1, :mRNA2, :mRNA3, :P1, :P2, :P3])
set_x0!(repressilator, [:mRNA1, :mRNA2, :mRNA3], fill(0, 3))
set_x0!(repressilator, [:P1, :P2, :P3], [5, 0, 15])
set_param!(repressilator, :n, 2.0)
-set_param!(repressilator, [:α, :α0, :β], [4400.0, 2.0, 4.0])
-set_time_bound!(repressilator, 15.0)
-
-export repressilator
+set_param!(repressilator, [:α, :α0, :β, :n], [400.0, 0.0, 2.0, 2.0])
+set_time_bound!(repressilator, 200.0)
diff --git a/models/square_wave_oscillator.jl b/models/square_wave_oscillator.jl
index e4bb08e..223a8c7 100644
--- a/models/square_wave_oscillator.jl
+++ b/models/square_wave_oscillator.jl
@@ -7,7 +7,7 @@ l_tr_square = [:Tu, :t1, :t2, :t3]
p_square = [1, 0, 0, 0, 0, 0, 5.0]
x0_square = [1, 0, 1]
t0_square = 0.0
-function square_wave_f!(xnplus1::Vector{Int}, l_t::Vector{Float64}, l_tr::Vector{Transition},
+function SquareWave_f!(xnplus1::Vector{Int}, l_t::Vector{Float64}, l_tr::Vector{Transition},
xn::Vector{Int}, tn::Float64, p::Vector{Float64})
if p[1] == 0 && p[3] == 0 && p[5] == 0
copyto!(xnplus1, xn)
@@ -56,12 +56,15 @@ function square_wave_f!(xnplus1::Vector{Int}, l_t::Vector{Float64}, l_tr::Vector
l_tr[1] = possible_transitions[transition]
end
end
-isabsorbing_square_wave(p::Vector{Float64}, xn::Vector{Int}) = (p[1] == 0 && p[3] == 0 && p[5] == 0)
+isabsorbing_SquareWave(p::Vector{Float64}, xn::Vector{Int}) = (p[1] == 0 && p[3] == 0 && p[5] == 0)
g_square_wave = [:A, :HIGH, :LOW]
-square_wave_oscillator = ContinuousTimeModel(d, k, dict_var_square, dict_params_square, l_tr_square,
- p_square, x0_square, t0_square, square_wave_f!, isabsorbing_square_wave;
- g=g_square_wave, name="square wave oscillator pkg", time_bound = 105.0)
+@everywhere @eval $(MarkovProcesses.generate_code_model_type_def(:SquareWaveModel))
+@everywhere @eval $(MarkovProcesses.generate_code_model_type_constructor(:SquareWaveModel))
+@everywhere @eval $(MarkovProcesses.generate_code_simulation(:SquareWaveModel, :SquareWave_f!, :isabsorbing_SquareWave))
-export square_wave_oscillator
+square_wave_oscillator = SquareWaveModel(d, k, dict_var_square, dict_params_square, l_tr_square,
+ p_square, x0_square, t0_square,
+ :SquareWave_f!, :isabsorbing_SquareWave;
+ g=g_square_wave, time_bound = 105.0)
diff --git a/tests/unit/simulate_available_models.jl b/tests/unit/simulate_available_models.jl
index 8c6bbec..814e6f9 100644
--- a/tests/unit/simulate_available_models.jl
+++ b/tests/unit/simulate_available_models.jl
@@ -5,11 +5,19 @@ load_model("SIR")
load_model("SIR_tauleap")
load_model("ER")
load_model("poisson")
+load_model("intracellular_viral_infection")
+load_model("repressilator")
+load_model("doping_3way_oscillator")
+load_model("square_wave_oscillator")
simulate(SIR)
simulate(ER)
simulate(SIR_tauleap)
simulate(poisson)
+simulate(intracellular_viral_infection)
+simulate(repressilator)
+simulate(doping_3way_oscillator)
+simulate(square_wave_oscillator)
return true
--
GitLab