Commit 5d886fc4 authored by Bentriou Mahmoud's avatar Bentriou Mahmoud
Browse files

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).
parent c809fe7e
......@@ -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
......
......@@ -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.Λ,
......
......@@ -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
......
This diff is collapsed.
......@@ -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"