model.jl 6.03 KB
Newer Older
1

2
3
load_model(name_model::String) = include(get_module_path() * "/models/" * name_model * ".jl")

4
# Simulation
5
function simulate(m::ContinuousTimeModel)
6
    # trajectory fields
7
8
9
    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
10
    times = Float64[m.t0]
11
    transitions = Transition[nothing]
12
13
    # values at time n
    n = 0
14
    xn = view(reshape(m.x0, 1, m.d), 1, :) # View for type stability
15
    tn = m.t0 
16
    # at time n+1
17
18
19
    mat_x = zeros(Int, m.buffer_size, m.d)
    l_t = zeros(Float64, m.buffer_size)
    l_tr = Vector{String}(undef, m.buffer_size)
20
    isabsorbing::Bool = m.isabsorbing(m.p,xn)
21
22
    end_idx = -1
    # use sizehint! ?
23
    while !isabsorbing && tn <= m.time_bound
24
        for i = 1:m.buffer_size
25
26
            m.f!(mat_x, l_t, l_tr, i, xn, tn, m.p)
            tn = l_t[i]
27
28
            if tn > m.time_bound
                i -= 1 # 0 is an ok value, 1:0 is allowed
29
                end_idx = i
30
31
32
33
                break
            end
            xn = view(mat_x, i, :)
            isabsorbing = m.isabsorbing(m.p,xn)
34
35
36
37
            if isabsorbing 
                end_idx = i
                break
            end
38
        end
39
40
41
42
43
44
45
46
47
48
49
50
        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))
51
    end
52
53
    append!(times, view(l_t, 1:end_idx))
    append!(transitions,  view(l_tr, 1:end_idx))
54
    if isbounded(m)
55
56
        for k = 1:m.d
            push!(full_values[k], full_values[k][end])
Bentriou Mahmoud's avatar
Bentriou Mahmoud committed
57
        end
58
59
        push!(times, m.time_bound)
        push!(transitions, nothing)
Bentriou Mahmoud's avatar
Bentriou Mahmoud committed
60
    end
61
    values = full_values[m._g_idx]
62
63
64
    return Trajectory(m, values, times, transitions)
end

65
66
67
68
function simulate(product::SynchronizedModel)
    # trajectory fields
    m = product.m
    A = product.automaton
69
70
71
    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
72
73
74
    times = Float64[m.t0]
    transitions = Union{String,Nothing}[nothing]
    reshaped_x0 = view(reshape(m.x0, 1, m.d), 1, :) # View for type stability
75
    S0 = init_state(A, m.x0, m.t0)
76
77
78
79
80
81
82
83
84
    # values at time n
    n = 0
    xn = reshaped_x0 
    tn = m.t0
    Sn = copy(S0)
    # at time n+1
    mat_x = zeros(Int, m.buffer_size, m.d)
    l_t = zeros(Float64, m.buffer_size)
    l_tr = Vector{String}(undef, m.buffer_size)
85
86
    isabsorbing::Bool = m.isabsorbing(m.p,xn)
    isacceptedLHA::Bool = isaccepted(Sn)
87
    Snplus1 = copy(Sn)
88
    while !isabsorbing && tn <= m.time_bound && !isacceptedLHA
89
        i = 0
90
        while i < m.buffer_size && !isabsorbing && tn <= m.time_bound && !isacceptedLHA
91
92
93
            i += 1
            m.f!(mat_x, l_t, l_tr, i, xn, tn, m.p)
            tn = l_t[i]
94
95
96
97
98
99
100
            if tn > m.time_bound
                i -= 1 # 0 is an ok value, 1:0 is allowed
                break
            end
            xn = view(mat_x, i, :)
            tr_n = l_tr[i]
            next_state!(Snplus1, A, xn, tn, tr_n, Sn)
101
            Sn = Snplus1
102
103
            isabsorbing = m.isabsorbing(m.p,xn)
            isacceptedLHA = isaccepted(Snplus1)
104
        end
105
106
107
        for k = 1:m.d
            append!(full_values[k], view(mat_x, 1:i, k))
        end
108
109
110
111
        append!(times, view(l_t, 1:i))
        append!(transitions,  view(l_tr, 1:i))
        n += i
    end
112
113
114
    # When the trajectory is accepted, we should not add an end value
    if isbounded(m) && !isaccepted(Sn)
        @assert times[end] < m.time_bound
115
116
117
118
119
        for k = 1:m.d
            push!(full_values[k], full_values[k][end])
        end
        push!(times, m.time_bound)
        push!(transitions, nothing)
120
    end
121
    values = full_values[m._g_idx]
122
    return SynchronizedTrajectory(Snplus1, product, values, times, transitions)
123
124
end

125
126
127
function simulate(m::ContinuousTimeModel, n::Int)
end

128
129
130
131
132
133
134
135
136
137
138
139
140
function set_observed_var!(m::Model,g::Vector{String})
    dobs = length(g)
    _map_obs_var_idx = Dict()
    _g_idx = Vector{Int}(undef, dobs)
    for i = 1:dobs
        _g_idx[i] = m.map_var_idx[g[i]] # = ( (g[i] = i-th obs var)::String => idx in state space )
        _map_obs_var_idx[g[i]] = i
    end
    m.g = g
    m._g_idx = _g_idx
    m._map_obs_var_idx = _map_obs_var_idx
end

141
isbounded(m::ContinuousTimeModel) = m.time_bound < Inf
142
143
144
145
146
147
148
function check_consistency(m::ContinuousTimeModel) 
    @assert m.d == length(m.map_var_idx) 
    @assert m.k == length(m.map_param_idx)
    @assert m.k == length(m.p)
    @assert length(m.g) <= m.d
    @assert length(m._g_idx) == length(m.g)
    @assert m.buffer_size >= 0
149
    @assert typeof(m.isabsorbing(m.p, view(reshape(m.x0, 1, m.d), 1, :))) == Bool
150
    return true
151
end
152

153
154
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)
155
156
157
158
159
160
161
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

162
get_param(m::ContinuousTimeModel) = m.p
163
getindex(m::ContinuousTimeModel, name_p::String) = m.p[m.map_param_idx[name_p]]
164
set_time_bound!(m::ContinuousTimeModel, b::Float64) = (m.time_bound = b)
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
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
188