model.jl 11 KB
Newer Older
1

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

5
6
load_model(name_model::String) = include(get_module_path() * "/models/" * name_model * ".jl")

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

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

29
30
"""
    `simulate(m)`
31

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

237
238
239
function simulate(m::ContinuousTimeModel, n::Int)
end

240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
isbounded(m::ContinuousTimeModel) = m.time_bound < Inf
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)
255
    dobs = length(g)
256
257
    _map_obs_var_idx = Dict{String}{Int}()
    _g_idx = zeros(Int, dobs)
258
259
260
261
262
263
264
265
    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
266
267
function observe_all!(am::Model)
    m = get_proba_model(am)
268
269
270
271
272
273
274
275
276
    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
277
278
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)
279
function set_param!(m::ContinuousTimeModel, l_name_p::Vector{String}, p::Vector{Float64}) 
280
    for i = eachindex(l_name_p)
281
282
283
284
        set_param!(m, l_name_p[i], p[i])
    end
end

285
get_param(m::ContinuousTimeModel) = m.p
286
getindex(m::ContinuousTimeModel, name_p::String) = m.p[m.map_param_idx[name_p]]
287
set_time_bound!(am::Model, b::Float64) = (get_proba_model(am).time_bound = b)
288
289
290
291
292
293
294
295

function getproperty(m::ContinuousTimeModel, sym::Symbol)
    if sym == :dobs
        return length(m.g)
    else
        return getfield(m, sym)
    end
end
296
297
298
299
300
301
302
303

function getproperty(pm::ParametricModel, sym::Symbol)
    if sym == :df
        return length(pm.l_param)
    else
        return getfield(pm, sym)
    end
end
304
305
get_proba_model(m::ContinuousTimeModel) = m
get_proba_model(sm::SynchronizedModel) = sm.m
306
307
get_proba_model(pm::ParametricModel) = get_proba_model(pm.m)
get_observed_var(m::Model) = get_proba_model(am).g
308
309

# Prior methods
310
311
312
313
314
draw_model!(pm::ParametricModel) = set_param!(get_proba_model(pm), pm.l_param, rand(pm.dist))
draw!(vec_p::AbstractVector{Float64}, pm::ParametricModel) = rand!(pm.dist, vec_p)
function draw!(mat_p::Matrix{Float64}, pm::ParametricModel)
    for i = 1:size(mat_p)[2]
        draw!(view(mat_p, :, i), pm)
315
316
    end
end
317

318
319
320
321
322
prior_density(vec_p::AbstractVector{Float64}, pm::ParametricModel) = pdf(pm.dist, vec_p)
function prior_density!(res_density::Vector{Float64}, mat_p::Matrix{Float64}, pm::ParametricModel)
    for i = eachindex(res_density)
        res_density[i] = prior_density(view(mat_p, :, i), pm)
    end
323
324
end

325
function check_bounds() end
326

327
328
# to do: simulate(pm), create_res_dir, check_bounds, ajouter un champ _param_idx pour pm.
#
329