Support for recurrence networks based on LightGraphs (#112)

heliosdrm · web-flow · commit df0a9a9db8ad · 2021-01-11T15:04:04.000+01:00
* add SimpleGraph and SimpleDiGraph * add Matrix constructores * interface for SimpleGraphs from RP * remove option to keep diagonal in SimpleGraph * Revert "remove option to keep diagonal in SimpleGraph" This reverts commit d365ec1. * fix calculation of SimpleGraph * rna function, do not export LightGraph names * v1.5 * add LightGraphs in tests * fix diameter
diff --git a/Project.toml b/Project.toml
@@ -1,12 +1,13 @@
 name = "RecurrenceAnalysis"
 uuid = "639c3291-70d9-5ea2-8c5b-839eba1ee399"
 repo = "https://github.com/JuliaDynamics/RecurrenceAnalysis.jl.git"
-version = "1.4.0"
+version = "1.5.0"
 
 [deps]
 DelayEmbeddings = "5732040d-69e3-5649-938a-b6b4f237613f"
 DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab"
 Distances = "b4f34e82-e78d-54a5-968a-f98e89d6e8f7"
+LightGraphs = "093fc24a-ae57-5d10-9952-331d41423f4d"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
@@ -17,16 +18,18 @@ UnicodePlots = "b8865327-cd53-5732-bb35-84acbb429228"
 [compat]
 DelayEmbeddings = "1.17, 1.18, 1.19"
 Distances = "0.8, 0.9, 0.10"
+LightGraphs = "1"
 StaticArrays = "0.8,0.9,0.10,0.11,0.12, 1.0"
 UnicodePlots = "0.3,1"
 julia = "1.5"
 
 [extras]
 DynamicalSystemsBase = "6e36e845-645a-534a-86f2-f5d4aa5a06b4"
+LightGraphs = "093fc24a-ae57-5d10-9952-331d41423f4d"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Statistics", "Test", "SparseArrays", "Random", "DynamicalSystemsBase"]
+test = ["Statistics", "Test", "SparseArrays", "Random", "DynamicalSystemsBase", "LightGraphs"]
diff --git a/src/RecurrenceAnalysis.jl b/src/RecurrenceAnalysis.jl
@@ -45,15 +45,19 @@ export embed,
        nmprt,
        rqa,
        sorteddistances,
-       transitivity,
-       @windowed
+       @windowed,
+       rna,
+       # deprecated:
+       transitivity  
 
 include("distance_matrix.jl")
 include("matrices.jl")
 include("plot.jl")
 include("histograms.jl")
 include("rqa.jl")
 include("radius.jl")
+include("graphs.jl")
+include("rna.jl")
 include("windowed.jl")
 include("deprecate.jl")
 
diff --git a/src/deprecate.jl b/src/deprecate.jl
@@ -64,3 +64,20 @@ function rqa(::Type{NamedTuple}, R; onlydiagonal=false, kwargs...)
         )
     end
 end
+
+function transitivity(R::AbstractRecurrenceMatrix)
+    @warn "`transitivity(x::AbstractRecurrenceMatrix)` is deprecated`, use `rna` to analyse network parameters"
+    if size(R, 1) ≠ size(R, 2)
+        @warn "Computing network transitivity of a non-square adjacency matrix is impossible"
+        return NaN
+    end
+    R² = R.data * R.data
+    numerator = zero(eltype(R²))
+    for col = 1:size(R,2)
+        rows = view(rowvals(R), nzrange(R,col))
+        for r = rows
+            numerator += R²[r, col]
+        end
+    end
+    trans = numerator / (sum(R²) - LinearAlgebra.tr(R²))
+end
diff --git a/src/graphs.jl b/src/graphs.jl
@@ -0,0 +1,39 @@
+using LightGraphs
+
+"""
+    SimpleGraph(R::AbstractRecurrenceMatrix)
+
+Create a recurrent network as a `SimpleGraph`, from the symmetric recurrence matrix `R`.
+That matrix is typically the result of calculating a [`RecurrenceMatrix`](@ref)
+or a [`JointRecurrenceMatrix`](@ref).
+
+The recurrence structure of `R` is interpreted as the adjacency matrix of an
+undirected complex network, where two different vertices are connected if they are
+neighbors in the embedded phase space, i.e.
+
+```math
+A_{i,j} = R_{i,j} - \\delta_{i,j}
+```
+
+Following this definition, diagonal points of `R` are ommited, i.e.
+the graph does not contain self-connected nodes.
+
+See the package [LightGraphs.jl](https://juliagraphs.org/LightGraphs.jl/stable/)
+for further options to work with `SimpleGraph` objects, besides the functions
+for Recurrence Network Analysis provided in this package.
+
+# References
+
+[1] : R.V. Donner *et al.* "Recurrence networks — a novel paradigm for nonlinear time series analysis",
+*New Journal of Physics* 12, 033025 (2010).
+
+[2] : R.V. Donner *et al.* "Complex Network Analysis of Recurrences", in:
+Webber, C.L. & Marwan N. (eds.) *Recurrence Quantification Analysis.
+Theory and Best Practices*, Springer, pp. 101-165 (2015).
+"""
+function LightGraphs.SimpleGraphs.SimpleGraph(R::AbstractRecurrenceMatrix)
+    graph = SimpleGraph(R.data)
+    delta = SimpleGraphFromIterator(Edge(v,v) for v=1:size(graph, 1))
+    graph = difference(graph, delta)
+    return graph
+end
diff --git a/src/matrices.jl b/src/matrices.jl
@@ -59,6 +59,13 @@ function colvals(x::SparseMatrixCSC)
 end
 colvals(x::ARM) = colvals(x.data)
 
+# Convert to matrices
+Base.Array{T}(R::ARM) where T = Matrix{T}(R.data)
+Base.Matrix{T}(R::ARM) where T = Matrix{T}(R.data)
+Base.Array(R::ARM) = Matrix(R.data)
+Base.Matrix(R::ARM) = Matrix(R.data)
+SparseArrays.SparseMatrixCSC{T}(R::ARM) where T = SparseMatrixCSC{T}(R.data)
+SparseArrays.SparseMatrixCSC(R::ARM) = SparseMatrixCSC(R.data)
 
 """
     RecurrenceMatrix(x, ε; kwargs...)
@@ -209,6 +216,7 @@ JointRecurrenceMatrix(args...; kwargs...) = JointRecurrenceMatrix{WithinRange}(a
 
 """
     JointRecurrenceMatrix(R1, R2; kwargs...)
+    
 Create a joint recurrence matrix from given recurrence matrices `R1, R2`.
 """
 function JointRecurrenceMatrix(
diff --git a/src/rna.jl b/src/rna.jl
@@ -0,0 +1,46 @@
+"""
+    rna(R)
+    rna(args...; kwargs...)
+
+Calculate a set of Recurrence Network parameters.
+The input `R` can be a symmetric recurrence matrix that is
+interpreted as the adjacency matrix of an undirected complex network,
+such that linked vertices are neighboring points in the phase space.
+
+Alternatively, the inputs can be a graph object or any 
+valid inputs to the `SimpleGraph` constructor of the
+[LightGraphs](https://github.com/JuliaGraphs/LightGraphs.jl) package.
+
+## Return
+
+The returned value is a dictionary that contains the following entries,
+with the corresponding global network properties[^1][^2]:
+* `:density`: edge density, approximately equivalent to the global recurrence rate in the phase space.
+* `:transitivity`: network transitivity, which describes the
+global clustering of points following Barrat's and Weigt's formulation [^3].
+* `:averagepath`: mean value of the shortest path lengths taken over
+all pairs of connected vertices, related to the average separation
+between points in the phase.
+* `:diameter`: maximum value of the shortest path lengths between
+pairs of connected vertices, related to the phase space diameter.
+
+## References
+
+[^1]: R.V. Donner *et al.* ["Recurrence networks — a novel paradigm for nonlinear time series analysis",
+*New Journal of Physics* 12, 033025 (2010)](https://doi.org/10.1088/1367-2630/12/3/033025)
+
+[^2]: R.V. Donner *et al.*, [The geometry of chaotic dynamics — a complex network perspective,
+*Eur. Phys. J.* B 84, 653–672 (2011)](https://doi.org/10.1140/epjb/e2011-10899-1)
+
+[^3]: A. Barrat & M. Weight, ["On the properties of small-world network models",
+*The European Physical Journal B* 13, 547–560 (2000)](https://doi.org/10.1007/s100510050067)
+"""
+function rna(args...; kwargs...)
+    graph = SimpleGraph(args...; kwargs...)
+    return Dict{Symbol, Float64}(
+        :density => density(graph),
+        :transitivity => global_clustering_coefficient(graph),
+        :averagepath => mean(1 ./ closeness_centrality(graph)),
+        :diameter => diameter(graph)
+    )
+end
diff --git a/src/rqa.jl b/src/rqa.jl
@@ -76,59 +76,6 @@ function _rrdenominator(R::M; theiler=0, kwargs...) where
     return k*(k+1)
 end
 
-"""
-    transitivity(R::AbstractRecurrenceMatrix) → T
-Returns the network transitivity `T` of the ε-recurrence network `R`. Here the
-recurrence plot `R` is identified as the network adjacency matrix `A`.
-
-## Description
-
-We quote from [^Donner2011], where the authors provide a complete description:
-Transitivity is related to fundamental algebraic relationships between triples
-of discrete objects. Specifically, in graph-theoretical terms, we identify the
-set `X` with the set of vertices `V` , and the relation `R` with the mutual
-adjacency of pairs of vertices. Hence, for a given vertex `i ∈ V` , transitivity
-refers to the fact that for two other vertices `j, k ∈ V` with
-`A_ij = A_ik = 1, A_jk = 1` also holds. In a general network, this is typically
-not the case for all vertices. Consequently, characterising the degree of
-transitivity (or, alternatively, the relative frequency of closed 3-loops, which
-are commonly referred to as triangles) with respect to 
-the whole network provides important information on the structural graph
-properties, which may be related to important general features of the underlying
-system.
-
-The network transitivity measures the fraction of closed triangles with respect to
-the number of linked triples of vertices in the whole network:
-
-```math
-\\mathcal{T} = \\frac{3 \\times \\textrm{number of triangles in the network}}{\\textrm{number of linked triples of vertices}}
-```
-or given the adjacency matrix `A`
-```math
-\\mathcal{T} = \\frac{trace(A^3)}{\\sum(A^2) - trace(A^2)}
-```
-
-## References
-
-[^Donner2011]: R.V. Donner *et al.*, [The geometry of chaotic dynamics — a complex network perspective, Eur. Phys. J. B 84, 653–672 (2011)](https://doi.org/10.1140/epjb/e2011-10899-1)
-"""
-function transitivity(R::ARM)
-    if size(R, 1) ≠ size(R, 2)
-        @warn "Computing network transitivity of a non-square adjacency matrix is impossible"
-        return NaN
-    end
-    R² = R.data * R.data
-    numerator = zero(eltype(R²))
-    for col = 1:size(R,2)
-        rows = view(rowvals(R), nzrange(R,col))
-        for r = rows
-            numerator += R²[r, col]
-        end
-    end
-    trans = numerator / (sum(R²) - LinearAlgebra.tr(R²))
-end
-
-
 ###########################################################################################
 # 0. Histograms
 ###########################################################################################
@@ -444,7 +391,6 @@ The returned value contains the following entries,
 which can be retrieved as from a dictionary (e.g. `results[:RR]`, etc.):
 
 * `:RR`: recurrence rate (see [`recurrencerate`](@ref))
-* `:TRANS`: transitivity (see [`transitivity`](@ref))
 * `:DET`: determinsm (see [`determinism`](@ref))
 * `:L`: average length of diagonal structures (see [`dl_average`](@ref))
 * `:Lmax`: maximum length of diagonal structures (see [`dl_max`](@ref))
@@ -466,7 +412,6 @@ less than the keyword `lminvert`) the average and maximum values
 (`:L`, `:Lmax`, `:TT`, `:Vmax`, `:MRT`)
 are returned as `0.0` but their respective entropies (`:ENTR`, `:VENTR`, `:RTE`)
 are returned as `NaN`.
-`:TRANS` is only returned when the input is a square matrix.
 
 ## Keyword Arguments
 
@@ -567,9 +512,6 @@ function rqa(::Type{Dict}, R; onlydiagonal=false, kwargs...)
             :RTE   => _rt_entropy(rthist),
             :NMPRT => maximum(rthist)
         )
-        if size(R, 1) == size(R, 2)
-            rqa_dict[:TRANS] = transitivity(R)
-        end
         return rqa_dict
     end
 end
diff --git a/src/windowed.jl b/src/windowed.jl
@@ -1,6 +1,5 @@
 const rqa_funs = [
     :recurrencerate,
-    :transitivity,
     :determinism,
     :dl_average,
     :dl_max,
diff --git a/test/deprecations.jl b/test/deprecations.jl
@@ -21,4 +21,6 @@ v = Dataset(rand(10, 2))
     @test_logs (:warn, "`JointRecurrenceMatrix{WithinRange}(x::AbstractMatrix, y, ε; kwargs...)` is deprecated, use `JointRecurrenceMatrix{WithinRange}(Dataset(x), y, ε; kwargs...)`") JointRecurrenceMatrix(xm, v, 0.5)
     @test_logs (:warn, "`JointRecurrenceMatrix{WithinRange}(x, y::AbstractMatrix, ε; kwargs...)` is deprecated, use `JointRecurrenceMatrix{WithinRange}(x, Dataset(y), ε; kwargs...)`") JointRecurrenceMatrix(v, ym, 0.5)
     @test_logs (:warn, "`JointRecurrenceMatrix{WithinRange}(x::AbstractMatrix, y::AbstractMatrix, ε; kwargs...)` is deprecated, use `JointRecurrenceMatrix{WithinRange}(Dataset(x), Dataset(y), ε; kwargs...)`") JointRecurrenceMatrix(xm, ym, 0.5)
+    # transitivity
+    @test_logs (:warn, "`transitivity(x::AbstractRecurrenceMatrix)` is deprecated`, use `rna` to analyse network parameters") transitivity(R)
 end
diff --git a/test/dynamicalsystems.jl b/test/dynamicalsystems.jl
@@ -1,5 +1,6 @@
 using RecurrenceAnalysis
 using DynamicalSystemsBase, Random, Statistics, SparseArrays
+using LightGraphs, LinearAlgebra
 using Test
 
 RA = RecurrenceAnalysis
@@ -48,7 +49,7 @@ dict_keys = ["Sine wave","White noise","Hénon (chaotic)","Hénon (periodic)"]
     ε = rqa_threshold[k]
     dmat = distancematrix(xe, ye)
     crmat = CrossRecurrenceMatrix(xe, ye, ε)
-    @test Matrix(crmat.data) == (dmat .≤ ε)
+    @test Matrix(crmat) == (dmat .≤ ε)
     rmat = RecurrenceMatrix(xe, ε; parallel = false)
     rmat_p = RecurrenceMatrix(xe, ε; parallel = true)
     jrmat = JointRecurrenceMatrix(xe, ye, ε; parallel = false)
@@ -117,4 +118,10 @@ dict_keys = ["Sine wave","White noise","Hénon (chaotic)","Hénon (periodic)"]
     @windowed rqaw = rqa(rmatw) width=50 step=40
     @test rqaw[:RR] == rrw
 
+    # Recurrence networks
+    graph = SimpleGraph(rmat)
+    amat = adjacency_matrix(graph)
+    @test amat == Matrix(rmat) - I
+    rna_dict = rna(rmat)
+    rna_dict[:density] == recurrencerate(rmat)
 end
diff --git a/test/smallmatrix.jl b/test/smallmatrix.jl
@@ -1,5 +1,5 @@
 using RecurrenceAnalysis, SparseArrays
-using Test
+using Test, Statistics
 
 cells = [(1,1),(2,2),(3,2),(6,2),(7,2),(8,2),(6,3),(7,3),(10,3),(2,4),(4,4),
     (1,5),(2,5),(3,5),(5,5),(9,5),(10,5),(2,6),(3,6),(6,6),(7,6),(10,6),
@@ -124,12 +124,64 @@ end
 ### Recurrence network
 @testset "Recurrence networks" begin
     # 7 edges, 12 linked triples, 1 triangle (1-2-4)
-    mat =  [0 1 0 1 0 0
-            1 0 1 1 1 0
-            0 1 0 0 0 0
-            1 1 0 0 0 1
-            0 1 0 0 0 1
-            0 0 0 1 1 0]
-    adjmat = RecurrenceMatrix(mat)
-    @test transitivity(adjmat) == 0.25 # (3/12)
+    adjmat7 =  [0 1 0 1 0 0
+                1 0 1 1 1 0
+                0 1 0 0 0 0
+                1 1 0 0 0 1
+                0 1 0 0 0 1
+                0 0 0 1 1 0]
+    rna_dict = rna(RecurrenceMatrix(adjmat7))
+    @test rna_dict[:density] == 7/15
+    @test rna_dict[:transitivity] == 0.25 # (3/12)
+    # taken from Donner et al. https://doi.org/10.1088/1367-2630/12/3/033025
+    adjmat10 = [0 0 0 1 0 0 0 1 0 0
+                0 0 0 0 1 0 0 0 1 0
+                0 0 0 0 0 1 0 0 0 1
+                1 0 0 0 0 0 1 0 0 0
+                0 1 0 0 0 0 0 1 0 0
+                0 0 1 0 0 0 0 0 1 0
+                0 0 0 1 0 0 0 0 0 1
+                1 0 0 0 1 0 0 0 0 0
+                0 1 0 0 0 1 0 0 0 0
+                0 0 1 0 0 0 1 0 0 0]
+    dismat10 = [0 3 4 1 2 5 2 1 4 3
+                3 0 3 4 1 2 5 2 1 4
+                4 3 0 3 4 1 2 5 2 1
+                1 4 3 0 3 4 1 2 5 2
+                2 1 4 3 0 3 4 1 2 5
+                5 2 1 4 3 0 3 4 1 2
+                2 5 2 1 4 3 0 3 4 1
+                1 2 5 2 1 4 3 0 3 4
+                4 1 2 5 2 1 4 3 0 3
+                3 4 1 2 5 2 1 4 3 0]
+    rna_dict = rna(RecurrenceMatrix(adjmat10))
+    @test rna_dict[:density] == 2/9
+    @test rna_dict[:transitivity] == 0
+    @test rna_dict[:averagepath] ≈ sum(dismat10) / (10*9)
+    @test rna_dict[:diameter] == maximum(dismat10)
+    adjmat9 =  [0 0 0 0 1 0 0 1 1
+                0 0 0 0 0 1 0 0 1
+                0 0 0 1 0 0 1 0 0
+                0 0 1 0 1 0 1 1 0
+                1 0 0 1 0 0 1 1 1
+                0 1 0 0 0 0 0 0 0
+                0 0 1 1 1 0 0 0 0
+                1 0 0 1 1 0 0 0 1
+                1 1 0 0 1 0 0 1 0]
+    dismat9 =  [0 1 2 2 1 2 3 1 1
+                1 0 2 3 1 1 2 2 1
+                2 2 0 4 1 1 1 3 2
+                2 3 4 0 3 4 5 1 3
+                1 1 1 3 0 1 2 2 1
+                2 1 1 4 1 0 1 3 2
+                3 2 1 5 2 1 0 4 3
+                1 2 3 1 2 3 4 0 2
+                1 1 2 3 1 2 3 2 0]
+    rna_dict = rna(RecurrenceMatrix(adjmat9))
+    triples = [3, 1, 1, 6, 10, 0, 3, 6, 6]
+    triangles = [3, 0, 1, 3, 5, 0, 2, 4, 3]
+    @test rna_dict[:density] == 7/18
+    @test rna_dict[:transitivity] ≈ sum(triangles) / sum(triples)
+    @test rna_dict[:averagepath] ≈ sum(dismat9) / (9*8)
+    @test rna_dict[:diameter] == maximum(dismat9)
 end
