model.jl 12.3 KB
Newer Older
1

2
import Random: rand, rand!
3
import Distributions: insupport, pdf
4

5
6
7
8
9
10
11
12
13
14
15
16
17
18
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
19
20
21
22
23
24
25
26

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

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

30
31
32
Simulates a model. If `m::SynchronizedModel`, the simulation is synchronized with a 
Linear Hybrid Automaton.
"""
33
function simulate(m::ContinuousTimeModel; p::Union{Nothing,AbstractVector{Float64}} = nothing)
34
35
36
37
    p_sim = m.p
    if p != nothing
        p_sim = p
    end
38
    # First alloc
39
    full_values = Vector{Vector{Int}}(undef, m.d)
40
    for i = eachindex(full_values) full_values[i] = zeros(Int, m.estim_min_states) end
41
42
43
    times = zeros(Float64, m.estim_min_states)
    transitions = Vector{Transition}(undef, m.estim_min_states)
    # Initial values
44
    for i = eachindex(full_values) full_values[i][1] = m.x0[i] end
45
46
47
48
    times[1] = m.t0
    transitions[1] = nothing
    # Values at time n
    n = 1
49
    xn = m.x0 
50
    tn = m.t0 
51
    # at time n+1
52
    isabsorbing::Bool = m.isabsorbing(p_sim,xn)
53
    # If x0 is absorbing
54
55
56
57
58
59
60
61
    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
62
63
    # Alloc of vectors where we stock n+1 values
    vec_x = zeros(Int, m.d)
64
65
    l_t = Float64[0.0]
    l_tr = Transition[nothing]
66
    # First we fill the allocated array
67
    for i = 2:m.estim_min_states
68
        m.f!(vec_x, l_t, l_tr, xn, tn, p_sim)
69
70
71
72
73
        tn = l_t[1]
        if tn > m.time_bound
            break
        end
        n += 1
74
75
        xn = vec_x
        _update_values!(full_values, times, transitions, xn, tn, l_tr[1], i)
76
        isabsorbing = m.isabsorbing(p_sim,xn)
77
78
79
80
        if isabsorbing 
            break
        end
    end
81
    # If simulation ended before the estimation of states
82
83
84
85
86
87
88
89
90
    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
91
    size_tmp = 0
92
    while !isabsorbing && tn <= m.time_bound
93
94
        # Alloc buffer
        _resize_trajectory!(full_values, times, transitions, m.estim_min_states+size_tmp+m.buffer_size)
95
96
97
        i = 0
        while i < m.buffer_size
            i += 1
98
            m.f!(vec_x, l_t, l_tr, xn, tn, p_sim)
99
            tn = l_t[1]
100
            if tn > m.time_bound
101
                i -= 1
102
103
                break
            end
104
105
106
            xn = vec_x
            _update_values!(full_values, times, transitions, 
                            xn, tn, l_tr[1], m.estim_min_states+size_tmp+i)
107
            isabsorbing = m.isabsorbing(p_sim,xn)
108
109
110
            if isabsorbing 
                break
            end
111
        end
112
        # If simulation ended before the end of buffer
113
        if i < m.buffer_size
114
            _resize_trajectory!(full_values, times, transitions, m.estim_min_states+size_tmp+i)
115
        end
116
        size_tmp += i
117
        n += i
118
    end
119
    values = full_values[m._g_idx]
120
    if isbounded(m)
121
122
123
        # 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
124
    end
125
126
    return Trajectory(m, values, times, transitions)
end
127
128
129
130

function simulate(product::SynchronizedModel; p::Union{Nothing,AbstractVector{Float64}} = nothing)
    m = product.m
    A = product.automaton
131
132
133
134
    p_sim = m.p
    if p != nothing
        p_sim = p
    end
135
    # First alloc
136
    full_values = Vector{Vector{Int}}(undef, m.d)
137
138
139
140
141
142
143
    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
144
    S0 = init_state(A, m.x0, m.t0)
145
146
    # Values at time n
    n = 1
147
    xn = m.x0 
148
    tn = m.t0 
149
    Sn = copy(S0)
150
    isabsorbing::Bool = m.isabsorbing(p_sim,xn)
151
    isacceptedLHA::Bool = isaccepted(Sn)
152
153
    # If x0 is absorbing
    if isabsorbing || isacceptedLHA 
154
155
156
157
158
159
160
        _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
161
162
    # Alloc of vectors where we stock n+1 values
    vec_x = zeros(Int, m.d)
163
164
    l_t = Float64[0.0]
    l_tr = Transition[nothing]
165
    Snplus1 = copy(Sn)
166
    # First we fill the allocated array
167
    for i = 2:m.estim_min_states
168
        m.f!(vec_x, l_t, l_tr, xn, tn, p_sim)
169
170
171
172
173
        tn = l_t[1]
        if tn > m.time_bound
            break
        end
        n += 1
174
        xn = vec_x
175
176
        tr_n = l_tr[1]
        next_state!(Snplus1, A, xn, tn, tr_n, Sn)
177
        _update_values!(full_values, times, transitions, xn, tn, tr_n, i)
178
        Sn = Snplus1
179
        isabsorbing = m.isabsorbing(p_sim,xn)
180
181
182
183
184
        isacceptedLHA = isaccepted(Snplus1)
        if isabsorbing || isacceptedLHA 
            break
        end
    end
185
    # If simulation ended before the estimation of states
186
187
188
189
190
191
    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
192
        return SynchronizedTrajectory(Sn, product, values, times, transitions)
193
194
    end
    # Otherwise, buffering system
195
    size_tmp = 0
196
    while !isabsorbing && tn <= m.time_bound && !isacceptedLHA
197
198
        # Alloc buffer
        _resize_trajectory!(full_values, times, transitions, m.estim_min_states+size_tmp+m.buffer_size)
199
        i = 0
200
        while i < m.buffer_size
201
            i += 1
202
            m.f!(vec_x, l_t, l_tr, xn, tn, p_sim)
203
            tn = l_t[1]
204
            if tn > m.time_bound
205
                i -= 1
206
207
                break
            end
208
209
            xn = vec_x
            tr_n = l_tr[1]
210
            next_state!(Snplus1, A, xn, tn, tr_n, Sn)
211
212
            _update_values!(full_values, times, transitions, 
                            xn, tn, tr_n, m.estim_min_states+size_tmp+i)
213
            Sn = Snplus1
214
            isabsorbing = m.isabsorbing(p_sim,xn)
215
            isacceptedLHA = isaccepted(Snplus1)
216
            if isabsorbing || isacceptedLHA
217
218
                break
            end
219
        end
220
        # If simulation ended before the end of buffer
221
        if i < m.buffer_size
222
            _resize_trajectory!(full_values, times, transitions, m.estim_min_states+size_tmp+i)
223
        end
224
        size_tmp += i
225
226
        n += i
    end
227
    values = full_values[m._g_idx]
228
    if isbounded(m) && !isaccepted(Sn)
229
230
231
        # Add last value: the convention is the last transition is nothing,
        # the trajectory is bounded
        _finish_bounded_trajectory!(values, times, transitions, m.time_bound)
232
    end
233
    return SynchronizedTrajectory(Sn, product, values, times, transitions)
234
235
end

236
237
238
239
240
241
242
243
244
245
246
247
248
249
"""
    `simulate(pm::ParametricModel, p_prior::AbstractVector{Float64})

Simulates the model contained in pm with p_prior values.
It simulates the model by taking the parameters contained in get_proba_model(pm).p and
replace the 1D parameters pm.l_param with p_prior.
"""
function simulate(pm::ParametricModel, p_prior::AbstractVector{Float64})
    full_p = copy(get_proba_model(pm).p)
    full_p[pm._param_idx] = p_prior
    
    return simulate(pm.m; p = full_p) 
end

250
251
252
function simulate(m::ContinuousTimeModel, n::Int)
end

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

# Set and get Model fields
function set_observed_var!(am::Model, g::Vector{String})
    m = get_proba_model(am)
268
    dobs = length(g)
269
270
    _map_obs_var_idx = Dict{String}{Int}()
    _g_idx = zeros(Int, dobs)
271
272
273
274
275
276
277
278
    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
279
280
function observe_all!(am::Model)
    m = get_proba_model(am)
281
282
283
284
285
286
287
288
289
    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
290
291
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)
292
function set_param!(m::ContinuousTimeModel, l_name_p::Vector{String}, p::Vector{Float64}) 
293
    for i = eachindex(l_name_p)
294
295
296
297
        set_param!(m, l_name_p[i], p[i])
    end
end

298
get_param(m::ContinuousTimeModel) = m.p
299
getindex(m::ContinuousTimeModel, name_p::String) = m.p[m.map_param_idx[name_p]]
300
set_time_bound!(am::Model, b::Float64) = (get_proba_model(am).time_bound = b)
301
302
303
304
305
306
307
308

function getproperty(m::ContinuousTimeModel, sym::Symbol)
    if sym == :dobs
        return length(m.g)
    else
        return getfield(m, sym)
    end
end
309
310
311
312
313
314
315
function getproperty(pm::ParametricModel, sym::Symbol)
    if sym == :df
        return length(pm.l_param)
    else
        return getfield(pm, sym)
    end
end
316

317
318
get_proba_model(m::ContinuousTimeModel) = m
get_proba_model(sm::SynchronizedModel) = sm.m
319
get_proba_model(pm::ParametricModel) = get_proba_model(pm.m)
320

321
get_observed_var(m::Model) = get_proba_model(am).g
322
323

# Prior methods
324
325
326
327
328
329
"""
    `draw_model!(pm::ParametricModel)`

Draw a parameter from the prior disitribution defined in `pm::ParametricModel`
and save it in the model contained in `pm`.
"""
330
draw_model!(pm::ParametricModel) = set_param!(get_proba_model(pm), pm.l_param, rand(pm.dist))
331
332
333
334
335
"""
    `draw!(vec_p, pm)`

Draw a parameter from the prior distribution defined in pm and stores it in vec_p.
"""
336
draw!(vec_p::AbstractVector{Float64}, pm::ParametricModel) = rand!(pm.dist, vec_p)
337
338
339
340
341
342
"""
    `draw!(mat_p, pm)`

Draw `size(mat_p)[2]` (number of columns of mat_p) parameters from the prior distribution 
defined in pm and stores it in mat_p.
"""
Bentriou Mahmoud's avatar
Bentriou Mahmoud committed
343
function draw!(mat_p::AbstractMatrix{Float64}, pm::ParametricModel)
344
345
    for i = 1:size(mat_p)[2]
        draw!(view(mat_p, :, i), pm)
346
347
    end
end
348
349
350
351
352
353
354
355
356
357
358
359
"""
    `prior_pdf(p_prior, pm)`
 
Computes the density of the prior distribution defined in pm on argument p_prior.
"""
prior_pdf(pm::ParametricModel, p_prior::AbstractVector{Float64}) = pdf(pm.dist, p_prior)
"""
    `prior_pdf(vec_res, mat_p, pm)`
 
Computes the density of the prior distribution defined in pm on each column
ov mat_p. Stores it in mat_p. (`length(vec_res) == size(mat_p)[2]`)
"""
Bentriou Mahmoud's avatar
Bentriou Mahmoud committed
360
function prior_pdf!(res_pdf::AbstractVector{Float64}, pm::ParametricModel, mat_p::AbstractMatrix{Float64})
361
362
    for i = eachindex(res_pdf)
        res_pdf[i] = prior_pdf(pm, view(mat_p, :, i))
363
    end
364
end
Bentriou Mahmoud's avatar
Bentriou Mahmoud committed
365
366
fill!(pm::ParametricModel, p_prior::AbstractVector{Float64}) = get_proba_model(pm).p[pm._param_idx] = p_prior
insupport(pm::ParametricModel, p_prior::AbstractVector{Float64}) = insupport(pm.dist, p_prior)
367

368
369
# to do: simulate(pm), create_res_dir, check_bounds, ajouter un champ _param_idx pour pm.
#
370