model.jl 4.29 KB
Newer Older
1

2
# Simulation
3
function simulate(m::ContinuousTimeModel)
4
    # trajectory fields
5
    full_values = Matrix{Int}(undef, 1, m.d)
6
7
8
    full_values[1,:] = m.x0
    times = Float64[m.t0]
    transitions = Union{String,Nothing}[nothing]
9
10
    # values at time n
    n = 0
11
    xn = view(reshape(m.x0, 1, m.d), 1, :) # View for type stability
12
    tn = m.t0 
13
    # at time n+1
14
15
16
    mat_x = zeros(Int, m.buffer_size, m.d)
    l_t = zeros(Float64, m.buffer_size)
    l_tr = Vector{String}(undef, m.buffer_size)
17
    is_absorbing::Bool = m.is_absorbing(m.p,xn)
18
    while !is_absorbing && (tn <= m.time_bound)
19
20
21
22
        i = 0
        while i < m.buffer_size && !is_absorbing && (tn <= m.time_bound)
            i += 1
            m.f!(mat_x, l_t, l_tr, i, xn, tn, m.p)
23
            xn = view(mat_x, i, :)
24
            tn = l_t[i]
25
            is_absorbing = m.is_absorbing(m.p,xn)
26
        end
27
28
29
        full_values = vcat(full_values, view(mat_x, 1:i, :))
        append!(times, view(l_t, 1:i))
        append!(transitions,  view(l_tr, 1:i))
30
        n += i
31
        is_absorbing = m.is_absorbing(m.p,xn)
32
    end
Bentriou Mahmoud's avatar
Bentriou Mahmoud committed
33
34
    if is_bounded(m)
        if times[end] > m.time_bound
35
            full_values[end,:] = full_values[end-1,:]
Bentriou Mahmoud's avatar
Bentriou Mahmoud committed
36
37
            times[end] = m.time_bound
            transitions[end] = nothing
38
        else
39
            full_values = vcat(full_values, reshape(full_values[end,:], 1, m.d))
40
41
            push!(times, m.time_bound)
            push!(transitions, nothing)
Bentriou Mahmoud's avatar
Bentriou Mahmoud committed
42
43
        end
    end
44
    values = view(full_values, :, m._g_idx)
45
46
47
    return Trajectory(m, values, times, transitions)
end

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
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
function simulate(product::SynchronizedModel)
    # trajectory fields
    m = product.m
    A = product.automaton
    full_values = Matrix{Int}(undef, 1, m.d)
    full_values[1,:] = m.x0
    times = Float64[m.t0]
    transitions = Union{String,Nothing}[nothing]
    reshaped_x0 = view(reshape(m.x0, 1, m.d), 1, :) # View for type stability
    S0 = init_state(A, reshaped_x0, t0) 
    # 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)
    is_absorbing::Bool = m.is_absorbing(m.p,xn)
    Snplus1 = copy(Sn)
    while !is_absorbing && (tn < m.time_bound)
        i = 0
        while i < m.buffer_size && !is_absorbing && (tn < m.time_bound)
            i += 1
            m.f!(mat_x, l_t, l_tr, i, xn, tn, m.p)
            xn = view(mat_x, i, :)
            tn = l_t[i]
            tr_n = l_tr[n]
            A.next_state!(Snplus1, A, xn, tn, tr_n, Sn)
            Sn = Snplus1
            is_absorbing = m.is_absorbing(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
        is_absorbing = m.is_absorbing(m.p,xn)
    end
    if is_bounded(m)
        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 = view(full_values, :, m._g_idx)
    return Trajectory(m, values, times, transitions)
end

102
103
104
105
106
107
108
109
function simulate(m::ContinuousTimeModel, n::Int)
    obs = ContinuousObservations(undef, n)
    for i = 1:n
        obs[i] = simulate(m)
    end
    return obs
end

110
111
112
113
114
115
116
117
118
119
120
121
122
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

Bentriou Mahmoud's avatar
Bentriou Mahmoud committed
123
is_bounded(m::Model) = m.time_bound < Inf
124
function check_consistency(m::Model) end
125
function simulate(m::Model, n::Int; bound::Float64 = Inf)::AbstractObservations end
126
127
function set_param!(m::Model, p::Vector{Float64})::Nothing end
function get_param(m::Model)::Vector{Float64} end
128

129
get_module_path() = dirname(dirname(pathof(@__MODULE__)))
130
function load_model(name_model::String)
131
    include(get_module_path() * "/models/" * name_model * ".jl")
132
133
end