model.jl 14.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 = 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
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 = copy(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
        copyto!(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
        copyto!(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
            copyto!(xn, vec_x)
209
            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
            copyto!(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
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
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)
    # If x0 is absorbing
    if isabsorbing || isacceptedLHA 
        return Sn
    end
    # 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)
    while !isabsorbing && tn <= m.time_bound && !isacceptedLHA
        if verbose
            @show vec_x
            @show Sn
        end
        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]
        next_state!(Snplus1, A, xn, tn, tr_n, Sn)
        copyto!(Sn, Snplus1)
        isabsorbing = m.isabsorbing(p_sim,xn)
        isacceptedLHA = isaccepted(Snplus1)
        n += 1
    end
    return Sn
end


284
285
286
287
288
289
290
291
292
293
294
295
296
"""
    `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
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
"""
    `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
315

316
317
318
function distribute_prob_accept_lha(sm::SynchronizedModel, str_var::String, nbr_sim::Int)
    sum_val = @distributed (+) for i = 1:nbr_sim Int(isaccepted(volatile_simulate(sm))) end
    return sum_val / nbr_sim
319
320
end

321
isbounded(m::Model) = get_proba_model(m).time_bound < Inf
322
323
324
325
326
327
328
329
330
331
332
333
334
335
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)
336
    dobs = length(g)
337
338
    _map_obs_var_idx = Dict{String}{Int}()
    _g_idx = zeros(Int, dobs)
339
340
341
342
343
344
345
346
    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
347
348
function observe_all!(am::Model)
    m = get_proba_model(am)
349
350
351
352
353
354
355
356
357
    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
358
359
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)
360
function set_param!(m::ContinuousTimeModel, l_name_p::Vector{String}, p::Vector{Float64}) 
361
    for i = eachindex(l_name_p)
362
363
364
365
        set_param!(m, l_name_p[i], p[i])
    end
end

366
get_param(m::ContinuousTimeModel) = m.p
367
getindex(m::ContinuousTimeModel, name_p::String) = m.p[m.map_param_idx[name_p]]
368
set_time_bound!(am::Model, b::Float64) = (get_proba_model(am).time_bound = b)
369
370
371
372
373
374
375
376

function getproperty(m::ContinuousTimeModel, sym::Symbol)
    if sym == :dobs
        return length(m.g)
    else
        return getfield(m, sym)
    end
end
377
378
379
380
381
382
383
function getproperty(pm::ParametricModel, sym::Symbol)
    if sym == :df
        return length(pm.l_param)
    else
        return getfield(pm, sym)
    end
end
384

385
386
get_proba_model(m::ContinuousTimeModel) = m
get_proba_model(sm::SynchronizedModel) = sm.m
387
get_proba_model(pm::ParametricModel) = get_proba_model(pm.m)
388

389
get_observed_var(m::Model) = get_proba_model(am).g
390
391

# Prior methods
392
393
394
395
396
397
"""
    `draw_model!(pm::ParametricModel)`

Draw a parameter from the prior disitribution defined in `pm::ParametricModel`
and save it in the model contained in `pm`.
"""
398
draw_model!(pm::ParametricModel) = set_param!(get_proba_model(pm), pm.l_param, rand(pm.dist))
399
400
401
402
403
"""
    `draw!(vec_p, pm)`

Draw a parameter from the prior distribution defined in pm and stores it in vec_p.
"""
404
draw!(vec_p::AbstractVector{Float64}, pm::ParametricModel) = rand!(pm.dist, vec_p)
405
406
407
408
409
410
"""
    `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
411
function draw!(mat_p::AbstractMatrix{Float64}, pm::ParametricModel)
412
413
    for i = 1:size(mat_p)[2]
        draw!(view(mat_p, :, i), pm)
414
415
    end
end
416
417
418
419
420
421
422
423
424
425
426
427
"""
    `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
428
function prior_pdf!(res_pdf::AbstractVector{Float64}, pm::ParametricModel, mat_p::AbstractMatrix{Float64})
429
430
    for i = eachindex(res_pdf)
        res_pdf[i] = prior_pdf(pm, view(mat_p, :, i))
431
    end
432
end
Bentriou Mahmoud's avatar
Bentriou Mahmoud committed
433
434
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)
435

436
437
# to do: simulate(pm), create_res_dir, check_bounds, ajouter un champ _param_idx pour pm.
#
438