model.jl 4.58 KB
Newer Older
1

2
3
load_model(name_model::String) = include(get_module_path() * "/models/" * name_model * ".jl")

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

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
102
103
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

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

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

125
126
127
128
129
130
131
132
133
134
is_bounded(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.is_absorbing(m.p, view(reshape(m.x0, 1, m.d), 1, :))) == Bool
    return true
135
end
136
137
138
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)
get_param(m::ContinuousTimeModel) = m.p
139