model.jl 7.8 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
19
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

20
# Simulation
21
function simulate(m::ContinuousTimeModel)
22
    # First alloc
23
    full_values = Vector{Vector{Int}}(undef, m.d)
24
25
26
27
28
29
30
31
32
    for i = 1:m.d 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 = 1:m.d full_values[i][1] = m.x0[i] end
    times[1] = m.t0
    transitions[1] = nothing
    # Values at time n
    n = 1
33
    xn = view(reshape(m.x0, 1, m.d), 1, :) # View for type stability
34
    tn = m.t0 
35
    # at time n+1
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
    isabsorbing::Bool = m.isabsorbing(m.p,xn)
    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
        for k = 1:m.d full_values[k][n] = xn[k] end
        times[n] = tn
        transitions[n] = l_tr[1]
        isabsorbing = m.isabsorbing(m.p,xn)
        if isabsorbing 
            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 Trajectory(m, values, times, transitions)
    end
    # Otherwise, buffering system
76
77
    mat_x = zeros(Int, m.buffer_size, m.d)
    l_t = zeros(Float64, m.buffer_size)
78
    l_tr = Vector{Transition}(undef, m.buffer_size)
79
    while !isabsorbing && tn <= m.time_bound
80
81
82
        i = 0
        while i < m.buffer_size
            i += 1
83
84
            m.f!(mat_x, l_t, l_tr, i, xn, tn, m.p)
            tn = l_t[i]
85
            if tn > m.time_bound
86
                i -= 1
87
88
89
90
                break
            end
            xn = view(mat_x, i, :)
            isabsorbing = m.isabsorbing(m.p,xn)
91
92
93
            if isabsorbing 
                break
            end
94
        end
95
96
97
98
99
100
101
102
        _resize_trajectory!(full_values, times, transitions, n+i)
        # Update values
        rng_traj = (n+1):(n+i) 
        rng_buffer = 1:i
        for k = 1:m.d full_values[k][rng_traj] = view(mat_x, rng_buffer, k) end
        times[rng_traj] = l_t[rng_buffer]
        transitions[rng_traj] = l_tr[rng_buffer]
        n += i
103
    end
104
    values = full_values[m._g_idx]
105
    if isbounded(m)
106
107
108
        # 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
109
    end
110
111
112
    return Trajectory(m, values, times, transitions)
end

113
114
115
116
function simulate(product::SynchronizedModel)
    # trajectory fields
    m = product.m
    A = product.automaton
117
118
119
    full_values = Vector{Vector{Int}}(undef, m.d)
    for i = 1:m.d full_values[i] = Int[m.x0[i]] end
    for i = 1:m.d sizehint!(full_values[i], m.estim_min_states) end
120
121
122
    times = Float64[m.t0]
    transitions = Union{String,Nothing}[nothing]
    reshaped_x0 = view(reshape(m.x0, 1, m.d), 1, :) # View for type stability
123
    S0 = init_state(A, m.x0, m.t0)
124
125
126
127
128
129
130
131
132
    # values at time n
    n = 0
    xn = reshaped_x0 
    tn = m.t0
    Sn = copy(S0)
    # at time n+1
    mat_x = zeros(Int, m.buffer_size, m.d)
    l_t = zeros(Float64, m.buffer_size)
    l_tr = Vector{String}(undef, m.buffer_size)
133
134
    isabsorbing::Bool = m.isabsorbing(m.p,xn)
    isacceptedLHA::Bool = isaccepted(Sn)
135
    Snplus1 = copy(Sn)
136
    while !isabsorbing && tn <= m.time_bound && !isacceptedLHA
137
        i = 0
138
        while i < m.buffer_size && !isabsorbing && tn <= m.time_bound && !isacceptedLHA
139
140
141
            i += 1
            m.f!(mat_x, l_t, l_tr, i, xn, tn, m.p)
            tn = l_t[i]
142
143
144
145
146
147
148
            if tn > m.time_bound
                i -= 1 # 0 is an ok value, 1:0 is allowed
                break
            end
            xn = view(mat_x, i, :)
            tr_n = l_tr[i]
            next_state!(Snplus1, A, xn, tn, tr_n, Sn)
149
            Sn = Snplus1
150
151
            isabsorbing = m.isabsorbing(m.p,xn)
            isacceptedLHA = isaccepted(Snplus1)
152
        end
153
154
155
        for k = 1:m.d
            append!(full_values[k], view(mat_x, 1:i, k))
        end
156
157
158
159
        append!(times, view(l_t, 1:i))
        append!(transitions,  view(l_tr, 1:i))
        n += i
    end
160
161
162
    # When the trajectory is accepted, we should not add an end value
    if isbounded(m) && !isaccepted(Sn)
        @assert times[end] < m.time_bound
163
164
165
166
167
        for k = 1:m.d
            push!(full_values[k], full_values[k][end])
        end
        push!(times, m.time_bound)
        push!(transitions, nothing)
168
    end
169
    values = full_values[m._g_idx]
170
    return SynchronizedTrajectory(Snplus1, product, values, times, transitions)
171
172
end

173
174
175
function simulate(m::ContinuousTimeModel, n::Int)
end

176
177
178
179
180
181
182
183
184
185
186
187
188
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

189
isbounded(m::ContinuousTimeModel) = m.time_bound < Inf
190
191
192
193
194
195
196
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
197
    @assert typeof(m.isabsorbing(m.p, view(reshape(m.x0, 1, m.d), 1, :))) == Bool
198
    return true
199
end
200

201
202
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)
203
204
205
206
207
208
209
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

210
get_param(m::ContinuousTimeModel) = m.p
211
getindex(m::ContinuousTimeModel, name_p::String) = m.p[m.map_param_idx[name_p]]
212
set_time_bound!(m::ContinuousTimeModel, b::Float64) = (m.time_bound = b)
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
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
236