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

249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
function volatile_simulate(product::SynchronizedModel; 
                           p::Union{Nothing,AbstractVector{Float64}} = nothing, verbose::Bool = false)
    m = product.m
    A = product.automaton
    p_sim = m.p
    if p != nothing
        p_sim = p
    end
    S0 = init_state(A, m.x0, m.t0)
    # Values at time n
    n = 1
    xn = copy(m.x0)
    tn = m.t0 
    Sn = copy(S0)
    isabsorbing::Bool = m.isabsorbing(p_sim,xn)
    isacceptedLHA::Bool = isaccepted(Sn)
    # Alloc of vectors where we stock n+1 values
    vec_x = zeros(Int, m.d)
    l_t = Float64[0.0]
    l_tr = Transition[nothing]
    Snplus1 = copy(Sn)
270
271
272
273
274
    # If x0 is absorbing
    if isabsorbing || isacceptedLHA 
        if !isacceptedLHA
            next_state!(Snplus1, A, xn, m.time_bound, nothing, Sn; verbose = verbose)
            copyto!(Sn, Snplus1)
275
        end
276
277
278
        return Sn
    end
    while !isabsorbing && tn <= m.time_bound && !isacceptedLHA
279
280
281
282
283
284
285
286
        m.f!(vec_x, l_t, l_tr, xn, tn, p_sim)
        tn = l_t[1]
        if tn > m.time_bound
            i -= 1
            break
        end
        copyto!(xn, vec_x)
        tr_n = l_tr[1]
287
        next_state!(Snplus1, A, xn, tn, tr_n, Sn; verbose = verbose)
288
289
290
291
292
        copyto!(Sn, Snplus1)
        isabsorbing = m.isabsorbing(p_sim,xn)
        isacceptedLHA = isaccepted(Snplus1)
        n += 1
    end
293
294
295
296
    if isabsorbing && !isacceptedLHA
        next_state!(Snplus1, A, xn, m.time_bound, nothing, Sn; verbose = verbose)
        copyto!(Sn, Snplus1)
    end
297
298
299
300
    return Sn
end


301
302
303
304
305
306
307
308
309
310
311
312
313
"""
    `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
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
"""
    `distribute_mean_value_lha(sm::SynchronizedModel, str_var::String, nbr_stim::Int)`

Distribute over workers the computation of the mean value 
of a LHA over `nbr_sim` simulations of the model.
"""
function distribute_mean_value_lha(sm::SynchronizedModel, str_var::String, nbr_sim::Int)
    sum_val = @distributed (+) for i = 1:nbr_sim volatile_simulate(sm)[str_var] end
    return sum_val / nbr_sim
end

function mean_value_lha(sm::SynchronizedModel, str_var::String, nbr_sim::Int)
    sum_val = 0.0 
    for i = 1:nbr_sim 
        sum_val += volatile_simulate(sm)[str_var] 
    end
    return sum_val / nbr_sim
end
332

333
function distribute_prob_accept_lha(sm::SynchronizedModel, nbr_sim::Int)
334
335
    sum_val = @distributed (+) for i = 1:nbr_sim Int(isaccepted(volatile_simulate(sm))) end
    return sum_val / nbr_sim
336
337
end

338
isbounded(m::Model) = get_proba_model(m).time_bound < Inf
339
340
341
342
343
344
345
346
347
348
349
350
351
352
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)
353
    dobs = length(g)
354
355
    _map_obs_var_idx = Dict{String}{Int}()
    _g_idx = zeros(Int, dobs)
356
357
358
359
360
361
362
363
    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
364
365
function observe_all!(am::Model)
    m = get_proba_model(am)
366
367
368
369
370
371
372
373
374
    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
375
376
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)
377
function set_param!(m::ContinuousTimeModel, l_name_p::Vector{String}, p::Vector{Float64}) 
378
    for i = eachindex(l_name_p)
379
380
381
382
        set_param!(m, l_name_p[i], p[i])
    end
end

383
get_param(m::ContinuousTimeModel) = m.p
384
getindex(m::ContinuousTimeModel, name_p::String) = m.p[m.map_param_idx[name_p]]
385
set_time_bound!(am::Model, b::Float64) = (get_proba_model(am).time_bound = b)
386
387
388
389
390
391
392
393

function getproperty(m::ContinuousTimeModel, sym::Symbol)
    if sym == :dobs
        return length(m.g)
    else
        return getfield(m, sym)
    end
end
394
395
396
397
398
399
400
function getproperty(pm::ParametricModel, sym::Symbol)
    if sym == :df
        return length(pm.l_param)
    else
        return getfield(pm, sym)
    end
end
401

402
403
get_proba_model(m::ContinuousTimeModel) = m
get_proba_model(sm::SynchronizedModel) = sm.m
404
get_proba_model(pm::ParametricModel) = get_proba_model(pm.m)
405

406
get_observed_var(m::Model) = get_proba_model(am).g
407
408

# Prior methods
409
410
411
412
413
414
"""
    `draw_model!(pm::ParametricModel)`

Draw a parameter from the prior disitribution defined in `pm::ParametricModel`
and save it in the model contained in `pm`.
"""
415
draw_model!(pm::ParametricModel) = set_param!(get_proba_model(pm), pm.l_param, rand(pm.dist))
416
417
418
419
420
"""
    `draw!(vec_p, pm)`

Draw a parameter from the prior distribution defined in pm and stores it in vec_p.
"""
421
draw!(vec_p::AbstractVector{Float64}, pm::ParametricModel) = rand!(pm.dist, vec_p)
422
423
424
425
426
427
"""
    `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
428
function draw!(mat_p::AbstractMatrix{Float64}, pm::ParametricModel)
429
430
    for i = 1:size(mat_p)[2]
        draw!(view(mat_p, :, i), pm)
431
432
    end
end
433
434
435
436
437
438
439
440
441
442
443
444
"""
    `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
445
function prior_pdf!(res_pdf::AbstractVector{Float64}, pm::ParametricModel, mat_p::AbstractMatrix{Float64})
446
447
    for i = eachindex(res_pdf)
        res_pdf[i] = prior_pdf(pm, view(mat_p, :, i))
448
    end
449
end
Bentriou Mahmoud's avatar
Bentriou Mahmoud committed
450
451
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)
452

453
454
# to do: simulate(pm), create_res_dir, check_bounds, ajouter un champ _param_idx pour pm.
#
455