Skip to content

Commit

Permalink
Resolving merge conflicts
Browse files Browse the repository at this point in the history
  • Loading branch information
@rusandris
rusandris committed Sep 7, 2023
2 parents 30d6149 + 2ffbfd0 commit 4ce2c5b
Show file tree
Hide file tree
Showing 12 changed files with 651 additions and 87 deletions.
4 changes: 3 additions & 1 deletion Project.toml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -6,6 +6,7 @@ version = "0.1.0"
[deps]
Cairo = "159f3aea-2a34-519c-b102-8c37f9878175"
Compose = "a81c6b42-2e10-5240-aca2-a61377ecd94b"
DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
DelayEmbeddings = "5732040d-69e3-5649-938a-b6b4f237613f"
DiscreteMarkovChains = "8abcb7ef-b365-4f7b-ac38-56893fb62f9f"
DynamicalSystemsBase = "6e36e845-645a-534a-86f2-f5d4aa5a06b4"
Expand All @@ -17,6 +18,7 @@ MetaGraphsNext = "fa8bd995-216d-47f1-8a91-f3b68fbeb377"
Permutations = "2ae35dd2-176d-5d53-8349-f30d82d94d4f"
Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
StatsPlots = "f3b207a7-027a-5e70-b257-86293d7955fd"


9 changes: 6 additions & 3 deletions src/StateTransitionNetworks.jl
View file
Original file line number Diff line number Diff line change
@@ -1,21 +1,24 @@
module StateTransitionNetworks

import Graphs: DiGraph,add_edge!, inneighbors,outneighbors, nv, ne, weights, is_strongly_connected,strongly_connected_components,edges,degree_histogram
import Graphs: DiGraph,add_edge!, inneighbors,outneighbors, nv, ne, weights, is_strongly_connected,strongly_connected_components,edges,degree_histogram,vertices
import StatsBase: sample, mean, var, Weights
using MetaGraphsNext
using GraphPlot,Cairo,Compose
import SparseArrays: spzeros
using DynamicalSystemsBase: DynamicalSystemsBase,poincaresos,StateSpaceSet
import LinearAlgebra: eigen, Diagonal
import LinearAlgebra: eigen, Diagonal, I, nullspace
import DataStructures: OrderedDict

include("create_stn.jl")
include("network_measures.jl")
include("plot_stn.jl")
include("timeseries_analysis.jl")
include("operations.jl")

export timeseries_to_grid, create_stn, check_stn!, prob_matrix, weight_matrix, calculate_weight_matrix, random_walk_on_stn, randomwalk_step,isnormalized
export timeseries_to_grid, create_stn, check_stn!, prob_matrix, weight_matrix, state_distribution, calculate_weight_matrix, random_walk_on_stn, randomwalk_step,isnormalized
export network_measures, sinai_kolmogorov_entropy, measure_convergence,lyapunov_measure
export plot_stn
export stn_analysis,read_bin,ndensity
export add_timeseries,are_equal,timeseries_to_common_grid

end
96 changes: 62 additions & 34 deletions src/create_stn.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,48 +4,47 @@ Discretizes a 2D timeseries/trajectory on a grid. Returns
a discrete timeseries containing the the cell coordinates and the list
of vertices with cell coordinates.
"""
function timeseries_to_grid(timeseries, grid)
function timeseries_to_grid(timeseries, grid; grid_edges = [])
M = zeros(grid,grid)
T = length(timeseries[:,1])
x_min = minimum(timeseries[:, 1])
y_min = minimum(timeseries[:, 2])
x_max = maximum(timeseries[:, 1])
y_max = maximum(timeseries[:, 2])

#partitioning with little extra space on both ends

dx = 0.5*(x_max-x_min)/grid
dy = 0.5*(y_max-y_min)/grid
if isempty(grid_edges)
x_min = minimum(timeseries[:, 1])
y_min = minimum(timeseries[:, 2])
x_max = maximum(timeseries[:, 1])
y_max = maximum(timeseries[:, 2])
else
x_min,y_min,x_max,y_max = grid_edges
end

x_grid = range(x_min-dx, x_max+dx, grid);
#partitioning between [`x_min`, `x_max`] into `grid` number of cells
x_grid = range(x_min, nextfloat(x_max, 100), grid+1);
x_min = x_grid[1]
y_grid = range(y_min-dy, y_max+dy, grid);
y_grid = range(y_min, nextfloat(y_max, 100), grid+1);
y_min = y_grid[1]

#arrays for space-discrete timeseries
x_n = Vector{Int64}(undef, T)
y_n = Vector{Int64}(undef, T)
num_vertex = 0
vertex_names = [];
vertex_names = OrderedDict{Tuple{Int64,Int64},Int64}()
x_n, y_n = [], [];

for row in timeseries
y = floor(Int,(row[2]-y_min)/Float64(y_grid.step))+1
x = floor(Int,(row[1]-x_min)/Float64(x_grid.step))+1
if M[y,x] == 0
num_vertex += 1
push!(vertex_names, [num_vertex, x, y])
vertex_names[(x,y)] = num_vertex
M[y,x] = 1
end
push!(x_n, x)
push!(y_n, y)
end
vertex_names = reduce(hcat, vertex_names)'
d_timeseries = [x_n y_n]
d_timeseries = collect(zip(x_n,y_n))
return d_timeseries, vertex_names
end


"""
create_stn(ts,grid::Int64,plane,idxs;make_ergodic=false, verbose=false,kwargs...) -> stn
Higher level function that creates a state transition network (STN) directly from the timeseries of a continuous system `ts` using `DynamicalSystemsBase.poincaresos`. The returned `stn` is a directed metagraph object. \\
Expand All @@ -60,6 +59,7 @@ function create_stn(ts,grid::Int64,plane,idxs;make_ergodic=false, verbose=false,

#psos
psection = DynamicalSystemsBase.poincaresos(DynamicalSystemsBase.StateSpaceSet(ts), plane; save_idxs=idxs,warning=true,kwargs...)

#method for time-discrete trajectory
create_stn(psection,grid; make_ergodic=make_ergodic, verbose=verbose)

Expand Down Expand Up @@ -90,23 +90,29 @@ Creates a state transition network (STN) using the discrete timeseries and verte
stn[i,j] -> (:prob => P[i,j],:weight => Q[i,j])
"""
function create_stn(discrete_timeseries,vertex_names;make_ergodic=false,verbose=false)
nr_vertices = length(vertex_names[:,1])
states = zeros(Int32,length(discrete_timeseries[:,1])) #discrete timeseries but only with vertex indices

for v in 1:nr_vertices
timepoints_at_v = findall(x -> x == vertex_names[v,2:end],[eachrow(discrete_timeseries)...]) #find all times when state was v
states[timepoints_at_v] .= v
end
nr_vertices = length(vertex_names)

#weight and transition probability matrices
Q = spzeros(nr_vertices, nr_vertices)
P = spzeros(nr_vertices, nr_vertices)

#probability distribution of states
states_distrib = zeros(nr_vertices)

#count transitions
next_states = circshift(states,-1)
for i in eachindex(states[1:end-1])
Q[states[i],next_states[i]] += 1
for i in eachindex(discrete_timeseries[1:end-1])
state = discrete_timeseries[i]
next_state = discrete_timeseries[i+1]
Q[vertex_names[state],vertex_names[next_state]] += 1
states_distrib[vertex_names[state]] += 1
end

#update end (final) state
end_state = discrete_timeseries[end]
states_distrib[vertex_names[end_state]] += 1

#normalize state distribution
states_distrib = states_distrib ./ sum(states_distrib)

#normalize Q and fill P by normalizing rows
Q = Q./sum(Q)
Expand All @@ -116,7 +122,7 @@ function create_stn(discrete_timeseries,vertex_names;make_ergodic=false,verbose=
stn = MetaGraph(
DiGraph(),
Int64,
Dict{Symbol, Int64},
Dict{Symbol, Union{Int64,Float64}},
Dict{Symbol, Float64},
nothing,
edge_data -> 1.0,
Expand All @@ -126,13 +132,12 @@ function create_stn(discrete_timeseries,vertex_names;make_ergodic=false,verbose=
#add edges and properties
#Properties: vertices -> Dict{Symbol,Int64}, edges -> Dict{Symbol,Float64}

for v in 1:nr_vertices
x,y = vertex_names[v,2:end]
stn[v] = Dict(:x => x,:y => y)
for state in keys(vertex_names)
stn[vertex_names[state]] = Dict(:x => state[1],:y => state[2],:prob => states_distrib[vertex_names[state]])
end

for i in 1:nr_vertices
for j in 1:nr_vertices
for i in 1:length(vertex_names)
for j in 1:length(vertex_names)
if P[i, j] != 0
stn[i,j] = Dict(:prob => P[i,j],:weight => Q[i,j])
end
Expand Down Expand Up @@ -160,14 +165,16 @@ function create_stn(P::AbstractMatrix;make_ergodic=false,verbose=false)
stn = MetaGraph(
DiGraph(),
Int64,
Dict{Symbol, Int64},
Dict{Symbol, Union{Int64,Float64}},
Dict{Symbol, Float64},
nothing,
edge_data -> 1.0,
0.0)

# add dummy coordinates
for v in 1:nr_vertices
stn[v] = Dict{Symbol, Int64}()
#stn[v] = Dict{Symbol, Int64}()
stn[v] = Dict(:x => v, :y => 0, :prob => 0.)
end

Q = calculate_weight_matrix(P)
Expand All @@ -185,6 +192,27 @@ function create_stn(P::AbstractMatrix;make_ergodic=false,verbose=false)
return stn,retcode
end


"""
state_distribution(stn) -> prob_states,pos_states
Returns the probability distribution of states and the positions of the vertices that correspond to these states.
"""
function state_distribution(stn)
prob_states = zeros(nv(stn)) #probability of vertices(states)
pos_states = zeros(Int32,nv(stn),2) #position of vertices

for v in collect(vertices(stn))
# v is the code of the vertex
v_label = label_for(stn,v)
prob_states[v] = stn[v_label][:prob]
pos_states[v,:] .= [stn[v_label][:x],stn[v_label][:y]]
end

return prob_states,pos_states
end



function prob_matrix(stn)
nr_vertices = nv(stn)
P = spzeros(Float64,(nr_vertices,nr_vertices))
Expand Down
49 changes: 22 additions & 27 deletions src/network_measures.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -46,9 +46,9 @@ end
Calculates the Sinai-Kolmogorov Entropy and Lyapunov measure of a STN
by using the analytical definitions of both quantities
"""
function network_measures(P::AbstractMatrix)
entropy, ret_code_entr = sinai_kolmogorov_entropy(P)
lyapunov, variance, covariance, ret_code_lyap = lyapunov_measure(P)
function network_measures(P::AbstractMatrix;x=nothing)
entropy, ret_code_entr = sinai_kolmogorov_entropy(P;x=x)
lyapunov, variance, covariance, ret_code_lyap = lyapunov_measure(P;x=x)
return entropy, lyapunov
end

Expand Down Expand Up @@ -109,21 +109,19 @@ end
Calculates analytically the Sinai-Kolmogorov entropy given the P transition probability matrix of the STN.
"""
function sinai_kolmogorov_entropy(P::AbstractMatrix)
λ, V = eigen(Matrix(P))
λ, X = eigen(transpose(Matrix(P)))

if real(λ[end]) ≈ 1
x = transpose(X[:,end]./sum(X[:,end]))
v = V[:,end]./V[1,end]
else
return -1, :StochasticMatrixError
function sinai_kolmogorov_entropy(P::AbstractMatrix;x=nothing)

if isnothing(x)
x = nullspace(Matrix(P' - I))
x = x ./sum(x)
end

v = ones(length(x))

L = Matrix(-log.(P))
replace!(L, Inf=>0.0)
L = P.*L
entropy = x*L*v
entropy = (x'*L*v)[1]
return real(entropy), :Success
end

Expand All @@ -133,28 +131,25 @@ end
Calculates analytically the Lyapunov measure given the the P transition probability matrix of the STN.
"""
function lyapunov_measure(P::AbstractMatrix)
λ, V = eigen(Matrix(P))
λ, X = eigen(transpose(Matrix(P)))
function lyapunov_measure(P::AbstractMatrix;x=nothing)

if real(λ[end]) ≈ 1
x = transpose(X[:,end]./sum(X[:,end]))
v = V[:,end]./V[1,end]
else
return -1, -1, -1, :StochasticMatrixError
if isnothing(x)
x = nullspace(Matrix(P' - I))
x = x ./sum(x)
end

x = x'
v = ones(length(x))

L = Matrix(-log.(P))
replace!(L, Inf=>0.0)
L2 = P.*L.^2
L = P.*L
I = Diagonal(ones(length(x)))
S = (I-v*x)+(P-v*x)*inv(I-P+v*x)
covariance = x*L*S*L*v
variance = x*L2*v-(x*L*v)^2
lyapunov = variance + 2*real(covariance)
covariance = (x*L*S*L*v)[1]
variance = (x*L2*v)[1] -(x*L*v)[1]^2
lyapunov = variance + 2*covariance
if imag(covariance) < 1.0e-3
return real(lyapunov), real(variance), real(covariance), :Success
return lyapunov, variance, covariance, :Success
else
@show covariance
return real(lyapunov), real(variance), real(covariance), :ComplexCovariancveWarning
Expand Down
Loading

0 comments on commit 4ce2c5b

Please sign in to comment.