model.jl 10.9 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
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
    
"""
    `simulate(m)`
22

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

252

253
254
255
function simulate(m::ContinuousTimeModel, n::Int)
end

256
257
258
259
260
261
262
263
264
265
266
267
268
function set_observed_var!(m::Model,g::Vector{String})
    dobs = length(g)
    _map_obs_var_idx = Dict()
    _g_idx = Vector{Int}(undef, dobs)
    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

269
isbounded(m::ContinuousTimeModel) = m.time_bound < Inf
270
271
272
273
274
275
276
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
277
    @assert typeof(m.isabsorbing(m.p, view(reshape(m.x0, 1, m.d), 1, :))) == Bool
278
    return true
279
end
280

281
282
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)
283
284
285
286
287
288
289
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

290
get_param(m::ContinuousTimeModel) = m.p
291
getindex(m::ContinuousTimeModel, name_p::String) = m.p[m.map_param_idx[name_p]]
292
set_time_bound!(m::ContinuousTimeModel, b::Float64) = (m.time_bound = b)
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
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
316