_tests_simulate.jl 8.54 KB
Newer Older
1

2
3
4
5
6
7
8
struct OldTrajectory <: AbstractTrajectory 
    m::ContinuousTimeModel
    values::Matrix{Int}
    times::Vector{Float64}
    transitions::Vector{Union{String,Nothing}}
end 

9
10
# File for benchmarking simulation and memory access of the package.

11
12
# Trajectories

13
_get_values_col(σ::OldTrajectory, var::String) = 
14
@view σ.values[(σ.m)._map_obs_var_idx[var],:] 
15
_get_values_row(σ::OldTrajectory, var::String) = 
16
@view σ.values[:,(σ.m)._map_obs_var_idx[var]] 
17

18
_get_state_col(σ::OldTrajectory, idx::Int) = 
19
@view σ.values[:,idx]
20
_get_state_row(σ::OldTrajectory, idx::Int) = 
21
@view σ.values[idx,:]
22

23
_get_value_col(σ::OldTrajectory, var::String, idx::Int) = 
24
σ.values[(σ.m)._map_obs_var_idx[var],idx] 
25
_get_value_row(σ::OldTrajectory, var::String, idx::Int) = 
26
27
28
σ.values[idx,(σ.m)._map_obs_var_idx[var]] 

# Model
29
30
31
32
33
34
35
36

function _simulate_col(m::ContinuousTimeModel)
    # trajectory fields
    full_values = zeros(m.d, 0)
    times = zeros(0)
    transitions = Vector{Union{String,Nothing}}(undef,0)
    # values at time n
    n = 0
37
    xn = @view m.x0[:]
38
39
40
41
42
    tn = m.t0 
    tr = [""]
    # at time n+1
    xnplus1 = zeros(Int, m.d)
    tnplus1 = zeros(Float64, 1)
43
44
    isabsorbing = (m.isabsorbing(m.p,xn))::Bool
    while !isabsorbing && (tn <= m.time_bound)
45
46
47
48
49
50
        m.f!(xnplus1, tnplus1, tr, xn, tn, m.p)
        full_values = hcat(full_values, xnplus1)
        push!(times, tnplus1[1])
        push!(transitions, tr[1])
        xn, tn = xnplus1, tnplus1[1]
        n += 1
51
        isabsorbing = m.isabsorbing(m.p,xn)::Bool
52
53
    end
    values = @view full_values[m._g_idx,:] 
54
    if isbounded(m)
55
56
57
58
59
60
        if times[end] > m.time_bound
            values[:, end] = values[:,end-1]
            times[end] = m.time_bound
            transitions[end] = nothing
        end
    end
61
    return OldTrajectory(m, values, times, transitions)
62
63
64
65
66
67
68
69
70
end

function _simulate_row(m::ContinuousTimeModel)
    # trajectory fields
    full_values = zeros(m.d, 0)
    times = zeros(0)
    transitions = Vector{Union{String,Nothing}}(undef,0)
    # values at time n
    n = 0
71
    xn = @view m.x0[:]
72
73
74
75
76
    tn = m.t0 
    tr = [""]
    # at time n+1
    xnplus1 = zeros(Int, m.d)
    tnplus1 = zeros(Float64, 1)
77
78
    isabsorbing = (m.isabsorbing(m.p,xn))::Bool
    while !isabsorbing && (tn <= m.time_bound)
79
80
81
82
83
84
        m.f!(xnplus1, tnplus1, tr, xn, tn, m.p)
        full_values = vcat(full_values, xnplus1)
        push!(times, tnplus1[1])
        push!(transitions, tr[1])
        xn, tn = xnplus1, tnplus1[1]
        n += 1
85
        isabsorbing = m.isabsorbing(m.p,xn)::Bool
86
87
    end
    values = @view full_values[m._g_idx,:] 
88
    if isbounded(m)
89
90
91
92
93
94
        if times[end] > m.time_bound
            values[:, end] = values[:,end-1]
            times[end] = m.time_bound
            transitions[end] = nothing
        end
    end
95
    return OldTrajectory(m, values, times, transitions)
96
97
98
99
100
101
102
103
104
105
end


function _simulate_col_buffer(m::ContinuousTimeModel; buffer_size::Int = 5)
    # trajectory fields
    full_values = zeros(m.d, 0)
    times = zeros(0)
    transitions = Vector{Union{String,Nothing}}(undef,0)
    # values at time n
    n = 0
106
    xn = @view m.x0[:]
107
108
109
110
111
    tn = m.t0 
    # at time n+1
    mat_x = zeros(Int, m.d, buffer_size)
    l_t = zeros(Float64, buffer_size)
    l_tr = Vector{String}(undef, buffer_size)
112
113
    isabsorbing = m.isabsorbing(m.p,xn)::Bool
    while !isabsorbing && (tn <= m.time_bound)
114
        i = 0
115
        while i < buffer_size && !isabsorbing && (tn <= m.time_bound)
116
117
118
119
            i += 1
            m.f!(mat_x, l_t, l_tr, i, xn, tn, m.p)
            xn = @view mat_x[:,i]
            tn = l_t[i]
120
            isabsorbing = m.isabsorbing(m.p,xn)::Bool
121
122
123
124
125
        end
        full_values = hcat(full_values, @view mat_x[:,1:i])
        append!(times, @view l_t[1:i])
        append!(transitions,  @view l_tr[1:i])
        n += i
126
        isabsorbing = m.isabsorbing(m.p,xn)::Bool
127
128
    end
    values = @view full_values[m._g_idx,:] 
129
    if isbounded(m)
130
131
132
133
134
135
        if times[end] > m.time_bound
            values[:, end] = values[:,end-1]
            times[end] = m.time_bound
            transitions[end] = nothing
        end
    end
136
    return OldTrajectory(m, values, times, transitions)
137
138
139
140
141
142
143
144
145
end

function _simulate_row_buffer(m::ContinuousTimeModel; buffer_size::Int = 5)
    # trajectory fields
    full_values = zeros(0, m.d)
    times = zeros(0)
    transitions = Vector{Union{String,Nothing}}(undef,0)
    # values at time n
    n = 0
146
    xn = @view m.x0[:]
147
148
149
150
151
    tn = m.t0 
    # at time n+1
    mat_x = zeros(Int, buffer_size, m.d)
    l_t = zeros(Float64, buffer_size)
    l_tr = Vector{String}(undef, buffer_size)
152
153
    isabsorbing = m.isabsorbing(m.p,xn)::Bool
    while !isabsorbing && (tn <= m.time_bound)
154
        i = 0
155
        while i < buffer_size && !isabsorbing && (tn <= m.time_bound)
156
157
            i += 1
            m.f!(mat_x, l_t, l_tr, i, xn, tn, m.p)
158
            xn = @view mat_x[i,:]
159
            tn = l_t[i]
160
            isabsorbing = m.isabsorbing(m.p,xn)::Bool
161
162
163
164
165
        end
        full_values = vcat(full_values, @view mat_x[1:i,:])
        append!(times, @view l_t[1:i])
        append!(transitions,  @view l_tr[1:i])
        n += i
166
        isabsorbing = m.isabsorbing(m.p,xn)::Bool
167
168
    end
    values = @view full_values[:,m._g_idx] 
169
    if isbounded(m)
170
171
172
173
174
175
        if times[end] > m.time_bound
            values[end,:] = values[end-1,:]
            times[end] = m.time_bound
            transitions[end] = nothing
        end
    end
176
    return OldTrajectory(m, values, times, transitions)
177
178
end

179
180
181
182
183
184
185
186
function _simulate_without_view(m::ContinuousTimeModel)
    # trajectory fields
    full_values = Matrix{Int}(undef, 1, m.d)
    full_values[1,:] = m.x0
    times = Float64[m.t0]
    transitions = Union{String,Nothing}[nothing]
    # values at time n
    n = 0
187
    xn = @view m.x0[:]
188
189
190
191
192
    tn = m.t0 
    # 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)
193
194
    isabsorbing = m.isabsorbing(m.p,xn)::Bool
    while !isabsorbing && (tn <= m.time_bound)
195
        i = 0
196
        while i < m.buffer_size && !isabsorbing && (tn <= m.time_bound)
197
            i += 1
198
            m.f!(mat_x, l_t, l_tr, i, @view(xn[:]), tn, m.p)
199
200
            xn = mat_x[i,:]
            tn = l_t[i]
201
            isabsorbing = m.isabsorbing(m.p,@view(xn[:]))::Bool
202
203
204
205
206
        end
        full_values = vcat(full_values, mat_x[1:i,:])
        append!(times, l_t[1:i])
        append!(transitions,  l_tr[1:i])
        n += i
207
        isabsorbing = m.isabsorbing(m.p,@view(xn[:]))::Bool
208
    end
209
    if isbounded(m)
210
211
212
213
214
215
216
217
218
219
220
        if times[end] > m.time_bound
            full_values[end,:] = full_values[end-1,:]
            times[end] = m.time_bound
            transitions[end] = nothing
        else
            full_values = vcat(full_values, reshape(full_values[end,:], 1, m.d))
            push!(times, m.time_bound)
            push!(transitions, nothing)
        end
    end
    values = full_values[:,m._g_idx] 
221
    return OldTrajectory(m, values, times, transitions)
222
223
end

224
225
226
227
228
229
230
231
232
233
234
235
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
# With trajectory values in Matrix type
function _simulate_27d56(m::ContinuousTimeModel)
    # trajectory fields
    full_values = Matrix{Int}(undef, 1, m.d)
    full_values[1,:] = m.x0
    times = Float64[m.t0]
    transitions = Union{String,Nothing}[nothing]
    # values at time n
    n = 0
    xn = view(reshape(m.x0, 1, m.d), 1, :) # View for type stability
    tn = m.t0 
    # 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)
    isabsorbing::Bool = m.isabsorbing(m.p,xn)
    # use sizehint! ?
    while !isabsorbing && tn <= m.time_bound
        i = 0
        while i < m.buffer_size && !isabsorbing && tn <= m.time_bound
            i += 1
            m.f!(mat_x, l_t, l_tr, i, xn, tn, m.p)
            tn = l_t[i]
            if tn > m.time_bound
                i -= 1 # 0 is an ok value, 1:0 is allowed
                break
            end
            xn = view(mat_x, i, :)
            isabsorbing = m.isabsorbing(m.p,xn)
        end
        full_values = vcat(full_values, view(mat_x, 1:i, :))
        append!(times, view(l_t, 1:i))
        append!(transitions,  view(l_tr, 1:i))
        n += i
    end
    if isbounded(m)
        #=
        if times[end] > m.time_bound
            full_values[end,:] = full_values[end-1,:]
            times[end] = m.time_bound
            transitions[end] = nothing
        else
        end
        =#
        full_values = vcat(full_values, reshape(full_values[end,:], 1, m.d))
        push!(times, m.time_bound)
        push!(transitions, nothing)
    end
    values = view(full_values, :, m._g_idx)
    return OldTrajectory(m, values, times, transitions)
end