Allow rqa and associated functions to allow SparseMatrixCSC. (#122)

* Allow `rqa` and associated functions to allow `SparseMatrixCSC`.

* `theiler=1` by default for `SparseMatrixCSC`

Author: rs7q5

Project.toml

@@ -1,7 +1,7 @@
 name = "RecurrenceAnalysis"
 uuid = "639c3291-70d9-5ea2-8c5b-839eba1ee399"
 repo = "https://github.com/JuliaDynamics/RecurrenceAnalysis.jl.git"
-version = "1.6.1"
+version = "1.6.2"
 
 [deps]
 DelayEmbeddings = "5732040d-69e3-5649-938a-b6b4f237613f"

src/matrices/plot.jl

@@ -7,7 +7,7 @@ export recurrenceplot, coordinates, grayscale
     coordinates(R) -> xs, ys
 Return the coordinates of the recurrence points of `R` (in indices).
 """
-function coordinates(R::ARM) # TODO: allow scan every ÷100 elements?
+function coordinates(R::SparseMatrixCSC) # TODO: allow scan every ÷100 elements?
    rows = rowvals(R)
    is = zeros(Int, nnz(R))
    js = zeros(Int, nnz(R))
@@ -22,8 +22,8 @@ function coordinates(R::ARM) # TODO: allow scan every ÷100 elements?
    end
    return is, js
 end
-
-recurrenceplot(R::ARM; kwargs...) = recurrenceplot(stdout, R::ARM; kwargs...)
+coordinates(R::ARM) = coordinates(R.data)
+recurrenceplot(R::Union{ARM,SparseMatrixCSC}; kwargs...) = recurrenceplot(stdout, R::Union{ARM,SparseMatrixCSC}; kwargs...)
 
 """
     recurrenceplot([io,] R; minh = 25, maxh = 0.5, ascii, kwargs...) -> u
@@ -44,7 +44,7 @@ Notice that the accuracy of this function drops drastically for matrices whose s
 is significantly bigger than the width and height of the display (assuming each
 index of the matrix is one character).
 """
-function recurrenceplot(io::IO, R::ARM; minh = 25, maxh = 0.5, ascii = nothing, kwargs...)
+function recurrenceplot(io::IO, R::Union{ARM,SparseMatrixCSC}; minh = 25, maxh = 0.5, ascii = nothing, kwargs...)
     @assert maxh ≤ 1
     h, w = displaysize(io)
     h = max(minh, round(Int, maxh * h)) # make matrix as long as half the screen (but not too short)
@@ -85,7 +85,7 @@ function recurrenceplot(io::IO, R::ARM; minh = 25, maxh = 0.5, ascii = nothing,
     )
 end
 
-function Base.show(io::IO, ::MIME"text/plain", R::ARM)
+function Base.show(io::IO, ::MIME"text/plain", R::Union{ARM,SparseMatrixCSC})
     a = recurrenceplot(io, R)
     show(io, a)
 end

src/rqa/histograms.jl

@@ -101,8 +101,8 @@ end
 
 deftheiler(x::Union{RecurrenceMatrix,JointRecurrenceMatrix}) = 1
 deftheiler(x::CrossRecurrenceMatrix) = 0
-
-function diagonalhistogram(x::ARM; lmin::Integer=2, theiler::Integer=deftheiler(x), kwargs...)
+deftheiler(x::SparseMatrixCSC) = 1
+function diagonalhistogram(x::Union{ARM,SparseMatrixCSC}; lmin::Integer=2, theiler::Integer=deftheiler(x), kwargs...)
     (theiler < 0) && throw(ErrorException(
         "Theiler window length must be greater than or equal to 0"))
     (lmin < 1) && throw(ErrorException("lmin must be 1 or greater"))
@@ -113,7 +113,13 @@ function diagonalhistogram(x::ARM; lmin::Integer=2, theiler::Integer=deftheiler(
     if issymmetric(x)
         # If theiler==0, the LOI is counted separately to avoid duplication
         if theiler == 0
-            loi_hist = verticalhistograms(CrossRecurrenceMatrix(hcat(diag(x,0))),lmin=lmin)[1]
+            if all(isequal(1),x.nzval)==false #CrossRecurrenceMatrix won't work if all values aren't equal to 1 and theiler=0
+                xtmp = deepcopy(x)
+                xtmp.nzval .= 1
+                loi_hist = verticalhistograms(CrossRecurrenceMatrix(hcat(diag(xtmp,0))),lmin=lmin)[1]
+            else
+                loi_hist = verticalhistograms(CrossRecurrenceMatrix(hcat(diag(x,0))),lmin=lmin)[1]
+            end
         end
         inside = (dv .< max(theiler,1))
         f = 2
@@ -140,7 +146,8 @@ function diagonalhistogram(x::ARM; lmin::Integer=2, theiler::Integer=deftheiler(
     return dh
 end
 
-function verticalhistograms(x::ARM;
+
+function verticalhistograms(x::Union{ARM,SparseMatrixCSC};
     lmin::Integer=2, theiler::Integer=deftheiler(x), distances=true, kwargs...)
     (theiler < 0) && throw(ErrorException(
         "Theiler window length must be greater than or equal to 0"))
@@ -199,7 +206,7 @@ estimator of the Poincaré recurrence times [2].
 recurrence quantifications", in: Webber, C.L. & N. Marwan (eds.), *Recurrence
 Quantification Analysis. Theory and Best Practices*, Springer, pp. 3-43 (2015).
 """
-function recurrencestructures(x::ARM;
+function recurrencestructures(x::Union{ARM,SparseMatrixCSC};
     diagonal=true, vertical=true, recurrencetimes=true,
     theiler=deftheiler(x), kwargs...)
 

src/rqa/rqa.jl

@@ -52,7 +52,7 @@ URL: https://www.nsf.gov/pubs/2005/nsf05057/nmbs/nmbs.pdf
 recurrence quantifications", in: Webber, C.L. & N. Marwan (eds.), *Recurrence
 Quantification Analysis. Theory and Best Practices*, Springer, pp. 3-43 (2015).
 """
-function recurrencerate(R::ARM; theiler::Integer=deftheiler(R), kwargs...)::Float64
+function recurrencerate(R::Union{ARM,SparseMatrixCSC}; theiler::Integer=deftheiler(R), kwargs...)::Float64
     (theiler < 0) && throw(ErrorException(
         "Theiler window length must be greater than or equal to 0"))
     if theiler == 0
@@ -69,6 +69,13 @@ end
 # Calculate the denominator for the recurrence rate
 _rrdenominator(R::ARM; theiler=0, kwargs...) = length(R)
 
+function _rrdenominator(R::SparseMatrixCSC; theiler=0, kwargs...)
+
+    (theiler == 0) && (return length(R))
+    k = size(R,1) - theiler
+    return k*(k+1)
+end
+
 function _rrdenominator(R::M; theiler=0, kwargs...) where
     M<:Union{RecurrenceMatrix,JointRecurrenceMatrix}
 
@@ -153,7 +160,7 @@ is the number of lines of length equal to ``l``.
 points inside the Theiler window (see [`rqa`](@ref) for the
 default values and usage of the keyword argument `theiler`).
 """
-function determinism(R::ARM; kwargs...)
+function determinism(R::Union{ARM,SparseMatrixCSC}; kwargs...)
     npoints = recurrencerate(R; kwargs...)*_rrdenominator(R; kwargs...)
     return _determinism(diagonalhistogram(R; kwargs...), npoints)
 end
@@ -171,8 +178,7 @@ end
 Calculate the divergence of the recurrence matrix `R`
 (actually the inverse of [`dl_max`](@ref)).
 """
-divergence(R::ARM; kwargs...) = ( 1.0/dl_max(R; kwargs...) )
-
+divergence(R::Union{ARM,SparseMatrixCSC}; kwargs...) = ( 1.0/dl_max(R; kwargs...) )
 
 """
     trend(R[; border=10, theiler])
@@ -217,10 +223,10 @@ Dynamical Systems", in: Riley MA & Van Orden GC, *Tutorials in Contemporary
 Nonlinear Methods for the Behavioral Sciences*, 2005, 26-94.
 https://www.nsf.gov/pubs/2005/nsf05057/nmbs/nmbs.pdf
 """
-trend(R::ARM; theiler=deftheiler(R), kwargs...) =
+trend(R::Union{ARM,SparseMatrixCSC}; theiler=deftheiler(R), kwargs...) =
     _trend(tau_recurrence(R); theiler=theiler, kwargs...)
-
-function tau_recurrence(R::ARM)
+    
+function tau_recurrence(R::Union{ARM,SparseMatrixCSC})
     n = minimum(size(R))
     rv = rowvals(R)
     rr_τ = zeros(n)
@@ -290,7 +296,7 @@ is the number of lines of length equal to ``v``.
 points inside the Theiler window (see [`rqa`](@ref) for the
 default values and usage of the keyword argument `theiler`).
 """
-function laminarity(R::ARM; kwargs...)
+function laminarity(R::Union{ARM,SparseMatrixCSC}; kwargs...)
     npoints = recurrencerate(R)*_rrdenominator(R; kwargs...)
     return _laminarity(verticalhistograms(R; kwargs...)[1], npoints)
 end
@@ -308,8 +314,7 @@ default values and usage of the keyword argument `theiler`).
 The trapping time is the average of the vertical line structures and thus equal
 to [`vl_average`](@ref).
 """
-trappingtime(R::ARM; kwargs...) = vl_average(R; kwargs...)
-
+trappingtime(R::Union{ARM,SparseMatrixCSC}; kwargs...) = vl_average(R; kwargs...)
 ###########################################################################################
 # 3. Based on recurrence times
 ###########################################################################################
@@ -327,8 +332,7 @@ default values and usage of the keyword argument `theiler`).
 
 Equivalent to [`rt_average`](@ref).
 """
-meanrecurrencetime(R::ARM; kwargs...) = rt_average(R; kwargs...)
-
+meanrecurrencetime(R::Union{ARM,SparseMatrixCSC}; kwargs...) = rt_average(R; kwargs...)
 
 """
     nmprt(R[; lmin=2, theiler])
@@ -349,8 +353,7 @@ of recurrence times [1].
 [DOI:10.1103/physreve.75.036222](https://doi.org/10.1103/physreve.75.036222)
 
 """
-nmprt(R::ARM; kwargs) = maximum(verticalhistograms(R; kwargs...)[2])
-
+nmprt(R::Union{ARM,SparseMatrixCSC}, kwargs...) = maximum(verticalhistograms(R;theiler=deftheiler(R), kwargs...)[2])
 ###########################################################################################
 # 4. All in one
 ###########################################################################################

test/dynamicalsystems.jl

@@ -124,4 +124,20 @@ dict_keys = ["Sine wave","White noise","Hénon (chaotic)","Hénon (periodic)"]
     @test amat == Matrix(rmat) - I
     rna_dict = rna(rmat)
     rna_dict[:density] == recurrencerate(rmat)
+
+    #test sparse matrices
+    rqapar_sparse = rqa(rmat.data, theiler=1, lmin=3, border=20)
+    rqadiag_sparse = rqa(rmat.data, theiler=1, lmin=3, border=20, onlydiagonal=true)
+    for p in keys(rqadiag_sparse)
+        @test rqapar_sparse[p]==rqadiag_sparse[p]
+    end
+
+    #check for erroring here with sparse matrices
+    @windowed(rrw2 = recurrencerate(rmat.data,theiler=1), width=50, step=40)
+    @windowed(rqaw2 = rqa(rmat.data,theiler=1), width=50, step =40)
+    @test rqaw2[:RR]==rrw2
+    rc1 = coordinates(rmat)
+    rc2 = coordinates(rmat.data)
+    @test isequal(rc1[1],rc2[1])
+    @test isequal(rc2[2],rc2[2])
 end