model.jl 4.83 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
20
    isabsorbing::Bool = m.isabsorbing(m.p,xn)
    while !isabsorbing && tn <= m.time_bound
21
        i = 0
22
        while i < m.buffer_size && !isabsorbing && tn <= m.time_bound
23
24
25
            i += 1
            m.f!(mat_x, l_t, l_tr, i, xn, tn, m.p)
            tn = l_t[i]
26
27
28
29
30
31
            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)
32
        end
33
34
35
        full_values = vcat(full_values, view(mat_x, 1:i, :))
        append!(times, view(l_t, 1:i))
        append!(transitions,  view(l_tr, 1:i))
36
        n += i
37
        #isabsorbing = m.isabsorbing(m.p,xn)
38
    end
39
40
    if isbounded(m)
        #=
Bentriou Mahmoud's avatar
Bentriou Mahmoud committed
41
        if times[end] > m.time_bound
42
            full_values[end,:] = full_values[end-1,:]
Bentriou Mahmoud's avatar
Bentriou Mahmoud committed
43
44
            times[end] = m.time_bound
            transitions[end] = nothing
45
        else
Bentriou Mahmoud's avatar
Bentriou Mahmoud committed
46
        end
47
48
49
50
        =#
        full_values = vcat(full_values, reshape(full_values[end,:], 1, m.d))
        push!(times, m.time_bound)
        push!(transitions, nothing)
Bentriou Mahmoud's avatar
Bentriou Mahmoud committed
51
    end
52
    values = view(full_values, :, m._g_idx)
53
54
55
    return Trajectory(m, values, times, transitions)
end

56
57
58
59
60
61
62
63
64
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
65
    S0 = init_state(A, reshaped_x0, m.t0) 
66
67
68
69
70
71
72
73
74
    # 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)
75
76
    isabsorbing::Bool = m.isabsorbing(m.p,xn)
    isacceptedLHA::Bool = isaccepted(Sn)
77
    Snplus1 = copy(Sn)
78
79
    Sn_dump = Sn
    while !isabsorbing && tn <= m.time_bound && !isacceptedLHA
80
        i = 0
81
        while i < m.buffer_size && !isabsorbing && tn <= m.time_bound && !isacceptedLHA
82
83
84
            i += 1
            m.f!(mat_x, l_t, l_tr, i, xn, tn, m.p)
            tn = l_t[i]
85
86
87
88
89
90
91
            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)
92
            Sn = Snplus1
93
94
            isabsorbing = m.isabsorbing(m.p,xn)
            isacceptedLHA = isaccepted(Snplus1)
95
96
97
98
99
100
        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
101
102
103
104
105
106
    # When the trajectory is accepted, we should not add an end value
    if isbounded(m) && !isaccepted(Sn)
        @assert times[end] < m.time_bound
        #full_values = vcat(full_values, reshape(full_values[end,:], 1, m.d))
        #push!(times, m.time_bound)
        #push!(transitions, nothing)
107
108
    end
    values = view(full_values, :, m._g_idx)
109
    return SynchronizedTrajectory(Snplus1, product, values, times, transitions)
110
111
end

112
113
114
function simulate(m::ContinuousTimeModel, n::Int)
end

115
116
117
118
119
120
121
122
123
124
125
126
127
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

128
isbounded(m::ContinuousTimeModel) = m.time_bound < Inf
129
130
131
132
133
134
135
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
136
    @assert typeof(m.isabsorbing(m.p, view(reshape(m.x0, 1, m.d), 1, :))) == Bool
137
    return true
138
end
139
140
141
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
142
set_time_bound!(m::ContinuousTimeModel, b::Float64) = (m.time_bound = b)
143