trajectory.jl 7.98 KB
Newer Older
1
2
3
4
5
6

abstract type AbstractTrajectory end
ContinuousObservations = AbstractVector{AbstractTrajectory}

struct Trajectory <: AbstractTrajectory
    m::ContinuousTimeModel
7
8
9
    values::Matrix{Int}
    times::Vector{Float64}
    transitions::Vector{Union{String,Nothing}}
10
11
end

12
# Operations
13
14
function +(σ1::AbstractTrajectory,σ2::AbstractTrajectory) end
function -(σ1::AbstractTrajectory,σ2::AbstractTrajectory) end
15
16
17

# Top-level Lp distance function
"""
18
    `dist_lp(l_σ1, l_σ2; verbose, p, str_stat_list, str_stat_trajectory)`   
19
20
21
22
23
24
25
26
27
28

Function that computes Lp distance between two set of any dimensional trajectories.
...
# Arguments
- `l_σ1::Vector{AbstractTrajectory}` is the first set of trajectories. l_σ2 is the second.
- `verbose::Bool` If `true`, launch a verbose execution of the computation. 
- `str_stat_list::String` allows to set the statistic to apply on the distances of the trajectories. 
It is salso called linkage function in clustering field. Only "mean" is available for now on.
...
"""
29
function dist_lp(l_σ1::Vector{<:AbstractTrajectory}, l_σ2::Vector{<:AbstractTrajectory};  
30
31
32
33
34
35
                 verbose::Bool = false, p::Int = 1, str_stat_list::String = "mean", str_stat_trajectory::String = "mean")
    if str_stat_list == "mean"
        return dist_lp_mean(l_σ1, l_σ2; verbose = verbose, p = p, str_stat_trajectory = str_stat_trajectory)
    end
    error("Unrecognized statistic in dist_lp")
end
36

37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
"""
    `dist_lp(σ1, σ2; verbose, p, str_stat)`   

Function that computes Lp distance between two trajectories of any dimension.
It computes Lp distances for each observed variable (contained in `get_obs_var(σ)`). and then apply a statistic on these distances.
Requires `get_obs_var(σ1) == get_obs_var(σ2)`, it is verified if they are simulated from the same model.
...
# Arguments
- `σ1::AbstractTrajectory` is the first trajectory. σ2 is the second.
- `verbose::Bool` If `true`, launch a verbose execution of the computation. 
- `str_stat::String` allows to set the statistic to apply on the distances computed  of the trajectories. 
Only "mean" is available for now on.
...
"""
function dist_lp(σ1::AbstractTrajectory, σ2::AbstractTrajectory;  
                 verbose::Bool = false, p::Int = 1, str_stat::String = "mean")
    if str_stat == "mean"
        return dist_lp_mean(σ1, σ2; verbose = verbose, p = p)
    end
    error("Unrecognized statistic in dist_lp")
end
58

59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
"""
    `dist_lp(σ1, σ2, var; verbose, p, str_stat)`   

Function that computes Lp distance between two trajectories of any dimension.
It computes Lp distances for each observed variable (contained in `get_obs_var(σ)`). and then apply a statistic on these distances.
Requires `get_obs_var(σ1) == get_obs_var(σ2)`, it is verified if they are simulated from the same model.
...
# Arguments
- `σ1::AbstractTrajectory` is the first trajectory. σ2 is the second.
- `var::String` is an observed variable. Have to be contained in `get_obs_var(σ1)` and `get_obs_var(σ2)`.
- `verbose::Bool` If `true`, launch a verbose execution of the computation. 
...
"""
function dist_lp(σ1::AbstractTrajectory, σ2::AbstractTrajectory, var::String; 
                 verbose::Bool = false, p::Int = 1)
    if !is_bounded(σ1) || !is_bounded(σ2)
        @warn "Lp distance computed on unbounded trajectories. Result should be wrong"
    end
77
    return dist_lp(σ1[var], σ1["times"], σ2[var], σ2["times"]; verbose = false, p = p)
78
79
80
end

# Distance function. Vectorized version
81
function dist_lp(x_obs::AbstractVector{Int}, t_x::Vector{Float64}, y_obs::AbstractVector{Int}, t_y::Vector{Float64}; 
82
83
84
85
86
87
88
89
90
91
92
93
94
95
                 verbose::Bool = false, p::Int = 1)
    current_y_obs = y_obs[1]
    current_t_y = t_y[2]
    idx = 1
    res = 0.0
    for i = 1:(length(x_obs)-1)
        if verbose
            @show i
            @show (t_x[i], x_obs[i])
            @show current_t_y
            @show t_x[i+1]
        end
        last_t_y = t_x[i]
        while current_t_y < t_x[i+1]
96
            rect =  abs(current_y_obs - x_obs[i])^p * (current_t_y - last_t_y)
97
98
99
            res += rect
            if verbose
                println("-- in loop :")
100
                println("-- add : $rect abs($current_y_obs - $(x_obs[i]))^p * ($current_t_y - $(last_t_y)) / $(abs(current_y_obs - x_obs[i])^p) * $(current_t_y - last_t_y)")
101
102
103
104
105
106
107
108
109
                print("-- ")
                @show current_y_obs, current_t_y
            end
            idx += 1
            last_t_y = t_y[idx]
            current_y_obs = y_obs[idx]
            current_t_y = t_y[idx+1]
        end
        last_t_read = max(last_t_y, t_x[i])
110
        rect = abs(current_y_obs - x_obs[i])^p * (t_x[i+1] - last_t_read)
111
112
        res += rect
        if verbose
113
            println("add : $rect abs($current_y_obs - $(x_obs[i]))^p * ($(t_x[i+1]) - $last_t_read) / $(abs(current_y_obs - x_obs[i])^p) * $(t_x[i+1] - last_t_read)")
114
115
116
117
118
119
            @show t_x[i+1], t_y[idx+1]
            @show y_obs[idx], x_obs[i]
            @show last_t_read
            @show current_y_obs
        end
    end
120
    return res^(1/p)
121
122
123
124
125
126
127
end
# For all the observed variables
function dist_lp_mean(σ1::AbstractTrajectory, σ2::AbstractTrajectory;  
                      verbose::Bool = false, p::Int = 1)
    if get_obs_var(σ1) != get_obs_var(σ2) error("Lp distances should be computed with the same observed variables") end
    res = 0.0
    for var in get_obs_var(σ1)
128
        res += dist_lp(σ1, σ2, var; verbose = verbose, p = p)
129
130
131
132
    end
    return res / length_obs_var(σ1)
end
# For a list of trajectories
133
function dist_lp_mean(l_σ1::Vector{<:AbstractTrajectory}, l_σ2::Vector{<:AbstractTrajectory};  
134
135
136
137
138
139
140
141
142
                 verbose::Bool = false, p::Int = 1, str_stat_trajectory::String = "mean")
    res = 0.0
    for σ1 in l_σ1
        for σ2 in l_σ2
            res += dist_lp(σ1, σ2; verbose = verbose, p = p, str_stat = str_stat_trajectory)
        end
    end
    return res / (length(l_σ1)*length(l_σ2))
end
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
# Get the value at any time
function _f_step(l_x::AbstractArray, l_t::AbstractArray, t::Float64)
    @assert length(l_x) == length(l_t)
    for i = 1:(length(l_t)-1)
        if l_t[i] <= t < l_t[i+1]
            return l_x[i]
        end
    end
    if l_t[end] == t
        return l_x[end]
    end
    return missing
end
# Riemann sum
function _riemann_sum(f::Function, t_begin::Real, t_end::Real, step::Float64)
    res = 0.0
    l_t = collect(t_begin:step:t_end)
    for i in 2:length(l_t)
        res += f(l_t[i]) * (l_t[i] - l_t[i-1])
    end
    return res
end

function check_consistency(σ::AbstractTrajectory)
    test_sizes = length(σ.times) == length(σ.transitions) == size(σ.values)[1]
    test_obs_var = length_obs_var(σ) == length(σ.m.g) == size(σ.values)[2]
    if !test_sizes error("Issue between sizes of values/times/transitions in this trajectory.") end
    if !test_obs_var error("Issue between sizes of observed variables in model and values") end
    return true
end
173

174
# Properties of the trajectory
175
176
177
178
length_states(σ::AbstractTrajectory) = length(σ.times)
length_obs_var(σ::AbstractTrajectory) = size(σ.values)[2]
get_obs_var(σ::AbstractTrajectory) = (σ.m).g
is_bounded(σ::AbstractTrajectory) = σ.transitions[end] == nothing 
Bentriou Mahmoud's avatar
Bentriou Mahmoud committed
179

180
181
182
183
184
185
# Access to trajectory values
get_var_values(σ::AbstractTrajectory, var::String) = 
@view σ.values[:,(σ.m)._map_obs_var_idx[var]] 
get_state(σ::AbstractTrajectory, idx::Int) = @view σ.values[idx,:]
get_value(σ::AbstractTrajectory, var::String, idx::Int) = 
σ.values[idx,(σ.m)._map_obs_var_idx[var]] 
186
function δ(σ1::AbstractTrajectory,t::Float64) end
187

188
189
190
# Get var values ["I"]
function getindex(σ::AbstractTrajectory, var::String)
    if var  == "times"
191
        return σ.times
192
    elseif var == "transitions"
193
194
        return σ.transitions
    else
195
196
197
198
        return get_var_values(σ, var)
    end
end
# Get i-th state [i]
199
getindex(σ::AbstractTrajectory, idx::Int) = get_state(σ, idx)
200
# Get i-th value of var ["I", idx]
201
function getindex(σ::AbstractTrajectory, var::String, idx::Int)
202
203
204
205
206
207
    if var  == "times"
        return σ.times[idx]
    elseif var == "transitions"
        return σ.transitions[idx]
    else
        return get_value(σ, var, idx)
208
209
    end
end
210