Implement of ModelPrior + test.

Add several convenient methods for parameters of model.

Next step: implementation of simulate with resize!.

bench/pygmalion/multiple_long_sim_er.jl

@@ -13,7 +13,7 @@ pygmalion.load_model(str_m)
 str_oml = "P,R,time"
 ll_om = split(str_oml, ",")
 x0 = State(95.0, 95.0, 0.0, 0.0, 0.0, 0.0)
-p_true = Parameters(1.0, 1.0, 1.0)
+p_true = Parameters(0.2, 40.0, 1.0)
 u = Control(20.0)
 tml = 1:400
 g_all = create_observation_function([ObserverModel(str_oml, tml)]) 
@@ -29,6 +29,9 @@ println("MarkovProcesses:")
 using MarkovProcesses
 MarkovProcesses.load_model("ER")
 ER.time_bound = 20.0
+ER.estim_min_states = 7000
+set_param!(ER, "k1", 0.2)
+set_param!(ER, "k2", 40.0)
 @timev MarkovProcesses.simulate(ER)
 
 println("Default buffer size=10")
@@ -40,6 +43,6 @@ println("Buffer size 100")
 ER.buffer_size = 100
 b1_100 = @benchmark for i=1:$(nb_sim) MarkovProcesses.simulate($ER) end
 b2_100 = @benchmark for i=1:$(nb_sim) MarkovProcesses.simulate(ER) end
-@show minimum(b1), mean(b1), maximum(b1)
-@show minimum(b2), mean(b2), maximum(b2)
+@show minimum(b1_100), mean(b1_100), maximum(b1_100)
+@show minimum(b2_100), mean(b2_100), maximum(b2_100)
 

core/MarkovProcesses.jl

@@ -4,10 +4,13 @@ import Base: +, -, *
 import Base: copy, getfield, getindex, lastindex, setindex!, getproperty, setproperty!
 
 import StaticArrays: SVector
+import Distributions: Distribution, Product, Uniform, Normal
+
+export Distribution, Product, Uniform, Normal
 
 # Common types and constructors
 export Observations, AbstractTrajectory, Trajectory, SynchronizedTrajectory
-export Model, ContinuousTimeModel, SynchronizedModel
+export Model, ContinuousTimeModel, SynchronizedModel, ModelPrior
 export LHA, StateLHA, Edge
 
 # Trajectory related methods
@@ -22,7 +25,7 @@ export load_automaton, get_index, get_value, length_var, isaccepted
 
 # Model related methods
 export simulate, set_param!, get_param, set_observed_var!
-export set_time_bound!, getproperty
+export set_time_bound!, getproperty, draw!
 export isbounded, isaccepted, check_consistency
 export load_model, get_module_path
 

core/common.jl

@@ -22,6 +22,7 @@ mutable struct ContinuousTimeModel <: Model
     isabsorbing::Function
     time_bound::Float64
     buffer_size::Int
+    estim_min_states::Int
 end
 
 struct Trajectory <: AbstractTrajectory
@@ -31,6 +32,11 @@ struct Trajectory <: AbstractTrajectory
     transitions::Vector{Transition}
 end
 
+struct ModelPrior
+    m::Model
+    map_l_param_dist::Dict{Vector{String},Distribution}
+end
+
 struct Edge
     transitions::Vector{Transition}
     check_constraints::Function
@@ -74,7 +80,7 @@ end
 function ContinuousTimeModel(d::Int, k::Int, map_var_idx::Dict, map_param_idx::Dict, l_transitions::Vector{String}, 
               p::Vector{Float64}, x0::Vector{Int}, t0::Float64, 
               f!::Function, isabsorbing::Function; 
-              g::Vector{String} = keys(map_var_idx), time_bound::Float64 = Inf, buffer_size::Int = 10)
+              g::Vector{String} = keys(map_var_idx), time_bound::Float64 = Inf, buffer_size::Int = 10, estim_min_states::Int = 50)
     dobs = length(g)
     _map_obs_var_idx = Dict()
     _g_idx = Vector{Int}(undef, dobs)
@@ -90,7 +96,7 @@ function ContinuousTimeModel(d::Int, k::Int, map_var_idx::Dict, map_param_idx::D
         @warn "You have possibly redefined a function Model.isabsorbing used in a previously instantiated model."
     end
 
-    return ContinuousTimeModel(d, k, map_var_idx, _map_obs_var_idx, map_param_idx, l_transitions, p, x0, t0, f!, g, _g_idx, isabsorbing, time_bound, buffer_size)
+    return ContinuousTimeModel(d, k, map_var_idx, _map_obs_var_idx, map_param_idx, l_transitions, p, x0, t0, f!, g, _g_idx, isabsorbing, time_bound, buffer_size, estim_min_states)
 end
 
 LHA(A::LHA, map_var::Dict{String,Int}) = LHA(A.l_transitions, A.l_loc, A.Λ, 
@@ -99,3 +105,11 @@ LHA(A::LHA, map_var::Dict{String,Int}) = LHA(A.l_transitions, A.l_loc, A.Λ,
 Base.:*(m::ContinuousTimeModel, A::LHA) = SynchronizedModel(m, A)
 Base.:*(A::LHA, m::ContinuousTimeModel) = SynchronizedModel(m, A)
 
+function ModelPrior(m::Model, priors::Tuple{Vector{String},Distribution}...)
+    map = Dict{Vector{String},Distribution}()
+    for prior in priors
+        map[prior[1]] = prior[2]
+    end
+    return ModelPrior(m, map)
+end
+

core/model.jl

@@ -4,10 +4,11 @@ load_model(name_model::String) = include(get_module_path() * "/models/" * name_m
 # Simulation
 function simulate(m::ContinuousTimeModel)
     # trajectory fields
-    full_values = Matrix{Int}(undef, 1, m.d)
-    full_values[1,:] = m.x0
+    full_values = Vector{Vector{Int}}(undef, m.d)
+    for i = 1:m.d full_values[i] = Int[m.x0[i]] end
+    for i = 1:m.d sizehint!(full_values[i], m.estim_min_states) end
     times = Float64[m.t0]
-    transitions = Union{String,Nothing}[nothing]
+    transitions = Transition[nothing]
     # values at time n
     n = 0
     xn = view(reshape(m.x0, 1, m.d), 1, :) # View for type stability
@@ -17,39 +18,47 @@ function simulate(m::ContinuousTimeModel)
     l_t = zeros(Float64, m.buffer_size)
     l_tr = Vector{String}(undef, m.buffer_size)
     isabsorbing::Bool = m.isabsorbing(m.p,xn)
+    end_idx = -1
+    # use sizehint! ?
     while !isabsorbing && tn <= m.time_bound
-        i = 0
-        while i < m.buffer_size && !isabsorbing && tn <= m.time_bound
-            i += 1
+        for i = 1:m.buffer_size
             m.f!(mat_x, l_t, l_tr, i, xn, tn, m.p)
             tn = l_t[i]
             if tn > m.time_bound
                 i -= 1 # 0 is an ok value, 1:0 is allowed
+                end_idx = i
                 break
             end
             xn = view(mat_x, i, :)
             isabsorbing = m.isabsorbing(m.p,xn)
+            if isabsorbing 
+                end_idx = i
+                break
+            end
         end
-        full_values = vcat(full_values, view(mat_x, 1:i, :))
-        append!(times, view(l_t, 1:i))
-        append!(transitions,  view(l_tr, 1:i))
-        n += i
-        #isabsorbing = m.isabsorbing(m.p,xn)
+        if end_idx != -1
+            break 
+        end
+        for k = 1:m.d
+            append!(full_values[k], view(mat_x, :, k))
+        end
+        append!(times, l_t)
+        append!(transitions,  l_tr)
+        n += m.buffer_size
+    end
+    for k = 1:m.d
+        append!(full_values[k], view(mat_x, 1:end_idx, k))
     end
+    append!(times, view(l_t, 1:end_idx))
+    append!(transitions,  view(l_tr, 1:end_idx))
     if isbounded(m)
-        #=
-        if times[end] > m.time_bound
-            full_values[end,:] = full_values[end-1,:]
-            times[end] = m.time_bound
-            transitions[end] = nothing
-        else
+        for k = 1:m.d
+            push!(full_values[k], full_values[k][end])
         end
-        =#
-        full_values = vcat(full_values, reshape(full_values[end,:], 1, m.d))
         push!(times, m.time_bound)
         push!(transitions, nothing)
     end
-    values = view(full_values, :, m._g_idx)
+    values = full_values[m._g_idx]
     return Trajectory(m, values, times, transitions)
 end
 
@@ -57,12 +66,13 @@ function simulate(product::SynchronizedModel)
     # trajectory fields
     m = product.m
     A = product.automaton
-    full_values = Matrix{Int}(undef, 1, m.d)
-    full_values[1,:] = m.x0
+    full_values = Vector{Vector{Int}}(undef, m.d)
+    for i = 1:m.d full_values[i] = Int[m.x0[i]] end
+    for i = 1:m.d sizehint!(full_values[i], m.estim_min_states) end
     times = Float64[m.t0]
     transitions = Union{String,Nothing}[nothing]
     reshaped_x0 = view(reshape(m.x0, 1, m.d), 1, :) # View for type stability
-    S0 = init_state(A, reshaped_x0, m.t0) 
+    S0 = init_state(A, m.x0, m.t0)
     # values at time n
     n = 0
     xn = reshaped_x0 
@@ -75,7 +85,6 @@ function simulate(product::SynchronizedModel)
     isabsorbing::Bool = m.isabsorbing(m.p,xn)
     isacceptedLHA::Bool = isaccepted(Sn)
     Snplus1 = copy(Sn)
-    Sn_dump = Sn
     while !isabsorbing && tn <= m.time_bound && !isacceptedLHA
         i = 0
         while i < m.buffer_size && !isabsorbing && tn <= m.time_bound && !isacceptedLHA
@@ -93,7 +102,9 @@ function simulate(product::SynchronizedModel)
             isabsorbing = m.isabsorbing(m.p,xn)
             isacceptedLHA = isaccepted(Snplus1)
         end
-        full_values = vcat(full_values, view(mat_x, 1:i, :))
+        for k = 1:m.d
+            append!(full_values[k], view(mat_x, 1:i, k))
+        end
         append!(times, view(l_t, 1:i))
         append!(transitions,  view(l_tr, 1:i))
         n += i
@@ -101,11 +112,13 @@ function simulate(product::SynchronizedModel)
     # When the trajectory is accepted, we should not add an end value
     if isbounded(m) && !isaccepted(Sn)
         @assert times[end] < m.time_bound
-        #full_values = vcat(full_values, reshape(full_values[end,:], 1, m.d))
-        #push!(times, m.time_bound)
-        #push!(transitions, nothing)
+        for k = 1:m.d
+            push!(full_values[k], full_values[k][end])
+        end
+        push!(times, m.time_bound)
+        push!(transitions, nothing)
     end
-    values = view(full_values, :, m._g_idx)
+    values = full_values[m._g_idx]
     return SynchronizedTrajectory(Snplus1, product, values, times, transitions)
 end
 
@@ -136,8 +149,40 @@ function check_consistency(m::ContinuousTimeModel)
     @assert typeof(m.isabsorbing(m.p, view(reshape(m.x0, 1, m.d), 1, :))) == Bool
     return true
 end
+
 set_param!(m::ContinuousTimeModel, p::Vector{Float64}) = (m.p = p)
 set_param!(m::ContinuousTimeModel, name_p::String, p_i::Float64) = (m.p[m.map_param_idx[name_p]] = p_i)
+function set_param!(m::ContinuousTimeModel, l_name_p::Vector{String}, p::Vector{Float64}) 
+    nb_param = length(l_name_p)
+    for i = 1:nb_param
+        set_param!(m, l_name_p[i], p[i])
+    end
+end
+
 get_param(m::ContinuousTimeModel) = m.p
+getindex(m::ContinuousTimeModel, name_p::String) = m.p[m.map_param_idx[name_p]]
 set_time_bound!(m::ContinuousTimeModel, b::Float64) = (m.time_bound = b)
+set_time_bound!(sm::SynchronizedModel, b::Float64) = set_time_bound!(sm.m, b)
+
+function getproperty(m::ContinuousTimeModel, sym::Symbol)
+    if sym == :dobs
+        return length(m.g)
+    else
+        return getfield(m, sym)
+    end
+end
+get_proba_model(m::ContinuousTimeModel) = m
+get_proba_model(sm::SynchronizedModel) = sm.m
+
+# Prior methods
+function draw!(Π::ModelPrior)
+    dict_dist = Π.map_l_param_dist
+    for l_name in keys(dict_dist)
+        if length(l_name) == 1
+            set_param!(get_proba_model(Π.m), l_name[1], rand(dict_dist[l_name]))
+        else
+            set_param!(get_proba_model(Π.m), l_name, rand(dict_dist[l_name]))
+        end
+    end
+end
 

tests/run_unit.jl

@@ -22,5 +22,6 @@ using Test
     @test include("unit/side_effects_1.jl")
     @test include("unit/create_models.jl")
     @test include("unit/create_automata.jl")
+    @test include("unit/model_prior.jl")
 end
 

tests/unit/model_prior.jl

+
+using MarkovProcesses
+load_model("ER")
+
+test_all = true
+
+load_model("ER")
+prior1 = ModelPrior(ER, (["k2"], Uniform(2.0, 4.0)))
+draw!(prior1)
+test_all = test_all && 2.0 <= ER["k2"] <= 4.0
+
+prior2 = ModelPrior(ER, (["k3","k2"], Product(Uniform.([2.5,6.0], [3.5,7.0]))))
+draw!(prior2)
+test_all = test_all && 2.5 <= ER["k3"] <= 3.5 && 6.0 <= ER["k2"] <= 7.0
+
+prior3 = ModelPrior(ER, (["k3"], Uniform(10.0, 11.0)), (["k2"], Uniform(13.0, 14.0)))
+draw!(prior3)
+test_all = test_all && 10.0 <= ER["k3"] <= 11.0 && 13.0 <= ER["k2"] <= 14.0
+
+return test_all
+

tests/unit/set_param.jl

@@ -22,6 +22,11 @@ k3 = 10.0
 set_param!(ER, "k3", 10.0)
 test_all = test_all && (ER.p == [0.5, 4.0, 10.0])
 
+l_k = [20.0, 10.0]
+l_name = ["k3", "k2"]
+set_param!(ER,  l_name, l_k)
+test_all = test_all && (ER.p == [0.5, 10.0, 20.0])
+
 new_param_SIR = [0.0013, 0.08]
 set_param!(SIR, new_param_SIR)
 test_all = test_all && (SIR.p == new_param_SIR)
@@ -34,5 +39,20 @@ ki = 0.011
 set_param!(SIR, "ki", ki)
 test_all = test_all && (SIR.p == [0.011, 0.06])
 
+l_name = ["kr", "ki"]
+l_p = [0.03, 0.015]
+set_param!(SIR, l_name, l_p)
+test_all = test_all && (SIR.p == [0.015, 0.03])
+
+l_name = ["ki", "kr"]
+l_p = [0.03, 0.015]
+set_param!(SIR, l_name, l_p)
+test_all = test_all && (SIR.p == [0.03, 0.015])
+
+l_name = ["kr"]
+l_p = [0.02]
+set_param!(SIR, l_name, l_p)
+test_all = test_all && (SIR.p == [0.03, 0.02])
+
 return test_all