model.jl 10.1 KB
Newer Older
1

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

4
5
6
7
8
9
10
11
12
13
14
15
16
17
function _resize_trajectory!(values::Vector{Vector{Int}}, times::Vector{Float64}, 
                             transitions::Vector{Transition}, size::Int)
    for i = eachindex(values) resize!(values[i], size) end
    resize!(times, size)
    resize!(transitions, size)
end


function _finish_bounded_trajectory!(values::Vector{Vector{Int}}, times::Vector{Float64}, 
                                    transitions::Vector{Transition}, time_bound::Float64)
    for i = eachindex(values) push!(values[i], values[i][end]) end
    push!(times, time_bound)
    push!(transitions, nothing)
end
18
19
20
21
22
23
24
25

function _update_values!(values::Vector{Vector{Int}}, times::Vector{Float64}, transitions::Vector{Transition},
                         xn::Vector{Int}, tn::Float64, tr_n::Transition, idx::Int)
    for k = eachindex(values) values[k][idx] = xn[k] end
    times[idx] = tn
    transitions[idx] = tr_n
end

26
27
"""
    `simulate(m)`
28

29
30
31
Simulates a model. If `m::SynchronizedModel`, the simulation is synchronized with a 
Linear Hybrid Automaton.
"""
32
function simulate(m::ContinuousTimeModel)
33
    # First alloc
34
    full_values = Vector{Vector{Int}}(undef, m.d)
35
    for i = eachindex(full_values) full_values[i] = zeros(Int, m.estim_min_states) end
36
37
38
    times = zeros(Float64, m.estim_min_states)
    transitions = Vector{Transition}(undef, m.estim_min_states)
    # Initial values
39
    for i = eachindex(full_values) full_values[i][1] = m.x0[i] end
40
41
42
43
    times[1] = m.t0
    transitions[1] = nothing
    # Values at time n
    n = 1
44
    xn = m.x0 # View for type stability
45
    tn = m.t0 
46
    # at time n+1
47
    isabsorbing::Bool = m.isabsorbing(m.p,xn)
48
    # If x0 is absorbing
49
50
51
52
53
54
55
56
    if isabsorbing
        _resize_trajectory!(full_values, times, transitions, 1)
        values = full_values[m._g_idx]
        if isbounded(m)
            _finish_bounded_trajectory!(values, times, transitions, m.time_bound)
        end
        return Trajectory(m, values, times, transitions)
    end
57
58
    # Alloc of vectors where we stock n+1 values
    vec_x = zeros(Int, m.d)
59
60
    l_t = Float64[0.0]
    l_tr = Transition[nothing]
61
    # First we fill the allocated array
62
    for i = 2:m.estim_min_states
63
        m.f!(vec_x, l_t, l_tr, xn, tn, m.p)
64
65
66
67
68
        tn = l_t[1]
        if tn > m.time_bound
            break
        end
        n += 1
69
70
        xn = vec_x
        _update_values!(full_values, times, transitions, xn, tn, l_tr[1], i)
71
72
73
74
75
        isabsorbing = m.isabsorbing(m.p,xn)
        if isabsorbing 
            break
        end
    end
76
    # If simulation ended before the estimation of states
77
78
79
80
81
82
83
84
85
    if n < m.estim_min_states
        _resize_trajectory!(full_values, times, transitions, n)
        values = full_values[m._g_idx]
        if isbounded(m)
            _finish_bounded_trajectory!(values, times, transitions, m.time_bound)
        end
        return Trajectory(m, values, times, transitions)
    end
    # Otherwise, buffering system
86
    size_tmp = 0
87
    while !isabsorbing && tn <= m.time_bound
88
89
        # Alloc buffer
        _resize_trajectory!(full_values, times, transitions, m.estim_min_states+size_tmp+m.buffer_size)
90
91
92
        i = 0
        while i < m.buffer_size
            i += 1
93
94
            m.f!(vec_x, l_t, l_tr, xn, tn, m.p)
            tn = l_t[1]
95
            if tn > m.time_bound
96
                i -= 1
97
98
                break
            end
99
100
101
            xn = vec_x
            _update_values!(full_values, times, transitions, 
                            xn, tn, l_tr[1], m.estim_min_states+size_tmp+i)
102
            isabsorbing = m.isabsorbing(m.p,xn)
103
104
105
            if isabsorbing 
                break
            end
106
        end
107
        # If simulation ended before the end of buffer
108
        if i < m.buffer_size
109
            _resize_trajectory!(full_values, times, transitions, m.estim_min_states+size_tmp+i)
110
        end
111
        size_tmp += i
112
        n += i
113
    end
114
    values = full_values[m._g_idx]
115
    if isbounded(m)
116
117
118
        # Add last value: the convention is the last transition is nothing,
        # the trajectory is bounded
        _finish_bounded_trajectory!(values, times, transitions, m.time_bound)
Bentriou Mahmoud's avatar
Bentriou Mahmoud committed
119
    end
120
121
    return Trajectory(m, values, times, transitions)
end
122
123
124
function simulate(product::SynchronizedModel)
    m = product.m
    A = product.automaton
125
    # First alloc
126
    full_values = Vector{Vector{Int}}(undef, m.d)
127
128
129
130
131
132
133
    for i = eachindex(full_values) full_values[i] = zeros(Int, m.estim_min_states) end
    times = zeros(Float64, m.estim_min_states)
    transitions = Vector{Transition}(undef, m.estim_min_states)
    # Initial values
    for i = eachindex(full_values) full_values[i][1] = m.x0[i] end
    times[1] = m.t0
    transitions[1] = nothing
134
    S0 = init_state(A, m.x0, m.t0)
135
136
    # Values at time n
    n = 1
137
    xn = m.x0 # View for type stability
138
    tn = m.t0 
139
    Sn = copy(S0)
140
141
    isabsorbing::Bool = m.isabsorbing(m.p,xn)
    isacceptedLHA::Bool = isaccepted(Sn)
142
143
    # If x0 is absorbing
    if isabsorbing || isacceptedLHA 
144
145
146
147
148
149
150
        _resize_trajectory!(full_values, times, transitions, 1)
        values = full_values[m._g_idx]
        if isbounded(m)
            _finish_bounded_trajectory!(values, times, transitions, m.time_bound)
        end
        return SynchronizedTrajectory(Sn, product, values, times, transitions)
    end
151
152
    # Alloc of vectors where we stock n+1 values
    vec_x = zeros(Int, m.d)
153
154
    l_t = Float64[0.0]
    l_tr = Transition[nothing]
155
    Snplus1 = copy(Sn)
156
    # First we fill the allocated array
157
    for i = 2:m.estim_min_states
158
        m.f!(vec_x, l_t, l_tr, xn, tn, m.p)
159
160
161
162
163
        tn = l_t[1]
        if tn > m.time_bound
            break
        end
        n += 1
164
        xn = vec_x
165
166
        tr_n = l_tr[1]
        next_state!(Snplus1, A, xn, tn, tr_n, Sn)
167
        _update_values!(full_values, times, transitions, xn, tn, tr_n, i)
168
169
170
171
172
173
174
        Sn = Snplus1
        isabsorbing = m.isabsorbing(m.p,xn)
        isacceptedLHA = isaccepted(Snplus1)
        if isabsorbing || isacceptedLHA 
            break
        end
    end
175
    # If simulation ended before the estimation of states
176
177
178
179
180
181
    if n < m.estim_min_states
        _resize_trajectory!(full_values, times, transitions, n)
        values = full_values[m._g_idx]
        if isbounded(m)
            _finish_bounded_trajectory!(values, times, transitions, m.time_bound)
        end
182
        return SynchronizedTrajectory(Sn, product, values, times, transitions)
183
184
    end
    # Otherwise, buffering system
185
    size_tmp = 0
186
    while !isabsorbing && tn <= m.time_bound && !isacceptedLHA
187
188
        # Alloc buffer
        _resize_trajectory!(full_values, times, transitions, m.estim_min_states+size_tmp+m.buffer_size)
189
        i = 0
190
        while i < m.buffer_size
191
            i += 1
192
193
            m.f!(vec_x, l_t, l_tr, xn, tn, m.p)
            tn = l_t[1]
194
            if tn > m.time_bound
195
                i -= 1
196
197
                break
            end
198
199
            xn = vec_x
            tr_n = l_tr[1]
200
            next_state!(Snplus1, A, xn, tn, tr_n, Sn)
201
202
            _update_values!(full_values, times, transitions, 
                            xn, tn, tr_n, m.estim_min_states+size_tmp+i)
203
            Sn = Snplus1
204
205
            isabsorbing = m.isabsorbing(m.p,xn)
            isacceptedLHA = isaccepted(Snplus1)
206
            if isabsorbing || isacceptedLHA
207
208
                break
            end
209
        end
210
        # If simulation ended before the end of buffer
211
        if i < m.buffer_size
212
            _resize_trajectory!(full_values, times, transitions, m.estim_min_states+size_tmp+i)
213
        end
214
        size_tmp += i
215
216
        n += i
    end
217
    values = full_values[m._g_idx]
218
    if isbounded(m) && !isaccepted(Sn)
219
220
221
        # Add last value: the convention is the last transition is nothing,
        # the trajectory is bounded
        _finish_bounded_trajectory!(values, times, transitions, m.time_bound)
222
    end
223
    return SynchronizedTrajectory(Sn, product, values, times, transitions)
224
225
end

226

227
228
229
function simulate(m::ContinuousTimeModel, n::Int)
end

230
function set_observed_var!(m::Model, g::Vector{String})
231
    dobs = length(g)
232
233
    _map_obs_var_idx = Dict{String}{Int}()
    _g_idx = zeros(Int, dobs)
234
235
236
237
238
239
240
241
242
    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

243
244
245
246
247
248
249
250
251
252
253
function observe_all!(m::Model)
    g = Vector{String}(undef, m.d)
    _g_idx = collect(1:m.d)
    for var in keys(m.map_var_idx)
        g[m.map_var_idx[var]] = var
    end
    m.g = g
    m._g_idx = _g_idx
    m._map_obs_var_idx = m.map_var_idx
end

254
isbounded(m::ContinuousTimeModel) = m.time_bound < Inf
255
256
257
258
259
260
261
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
262
    @assert typeof(m.isabsorbing(m.p, m.x0)) == Bool
263
    return true
264
end
265

266
267
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)
268
269
270
271
272
273
274
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

275
get_param(m::ContinuousTimeModel) = m.p
276
getindex(m::ContinuousTimeModel, name_p::String) = m.p[m.map_param_idx[name_p]]
277
set_time_bound!(m::ContinuousTimeModel, b::Float64) = (m.time_bound = b)
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
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
301