Fix broken argument in windowed function (#150)

* Fix broken argument in windowed function

Hopefully this one fixes both the function and the associated docs

* Revert formatting

* Replace all SparseMatrixCSC with AbstractMatrix

* Remove more SparseMatrixCSC

* Added condition for recurrence counts depending on a type of argument

* Removed views from windowed function

* SparseMatrix to SparseMatrixCSC

* A small fix to how windowing works

* Added some argtype declarations

* A couple of sanity checks

* Some more argtype declaration

* Tests for windowed function

* 2.0.4 => 2.0.5

Project.toml

@@ -1,7 +1,7 @@
 name = "RecurrenceAnalysis"
 uuid = "639c3291-70d9-5ea2-8c5b-839eba1ee399"
 repo = "https://github.com/JuliaDynamics/RecurrenceAnalysis.jl.git"
-version = "2.0.4"
+version = "2.0.5"
 
 [deps]
 DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab"

src/rqa/histograms.jl

@@ -107,8 +107,8 @@ end
 
 deftheiler(x::Union{RecurrenceMatrix,JointRecurrenceMatrix}) = 1
 deftheiler(x::CrossRecurrenceMatrix) = 0
-deftheiler(x::SparseMatrixCSC) = 1
-function diagonalhistogram(x::Union{ARM,SparseMatrixCSC}; lmin::Integer=2, theiler::Integer=deftheiler(x), kwargs...)
+deftheiler(x::AbstractMatrix) = 1
+function diagonalhistogram(x::Union{ARM,AbstractMatrix}; 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"))
@@ -153,7 +153,7 @@ function diagonalhistogram(x::Union{ARM,SparseMatrixCSC}; lmin::Integer=2, theil
 end
 
 
-function verticalhistograms(x::Union{ARM,SparseMatrixCSC};
+function verticalhistograms(x::Union{ARM,AbstractMatrix};
     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"))
@@ -212,7 +212,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::Union{ARM,SparseMatrixCSC};
+function recurrencestructures(x::Union{ARM,AbstractMatrix};
     diagonal=true, vertical=true, recurrencetimes=true,
     theiler=deftheiler(x), kwargs...)
 

src/rqa/rqa.jl

@@ -52,24 +52,29 @@ 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::Union{ARM,SparseMatrixCSC}; theiler::Integer=deftheiler(R), kwargs...)::Float64
+function recurrencerate(R::Union{ARM,AbstractMatrix}; theiler::Integer=deftheiler(R), kwargs...)::Float64
     (theiler < 0) && throw(ErrorException(
         "Theiler window length must be greater than or equal to 0"))
+    if R isa Union{ARM, SparseMatrixCSC}
+        n_recs = nnz(R)
+    else
+        n_recs = count(R)
+    end
     if theiler == 0
-        return nnz(R)/length(R)
+        return n_recs/length(R)
     end
     diags_remove = -(theiler-1):(theiler-1)
     theiler_nz = 0
     for d in diags_remove
         theiler_nz += nnz(diag(R,d))
     end
-    return (nnz(R) - theiler_nz)/_rrdenominator(R; theiler=theiler, kwargs...)
+    return (n_recs - theiler_nz)/_rrdenominator(R; theiler=theiler, kwargs...)
 end
 
 # Calculate the denominator for the recurrence rate
 _rrdenominator(R::ARM; theiler=0, kwargs...) = length(R)
 
-function _rrdenominator(R::SparseMatrixCSC; theiler=0, kwargs...)
+function _rrdenominator(R::AbstractMatrix; theiler=0, kwargs...)
 
     (theiler == 0) && (return length(R))
     k = size(R,1) - theiler
@@ -143,7 +148,7 @@ macro histogram_params(keyword, description, hist_fun)
         end
         push!(ret.args, quote
             @doc $doc ->
-            $fname(R::ARM; kwargs...) = $_fname($hist_fun(R; kwargs...))
+            $fname(R::Union{ARM,AbstractMatrix}; kwargs...) = $_fname($hist_fun(R::Union{ARM,AbstractMatrix}; kwargs...))
 
             function $_fname(hist::Vector{<:Integer})
                 $fbody
@@ -179,7 +184,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::Union{ARM,SparseMatrixCSC}; kwargs...)
+function determinism(R::Union{ARM,AbstractMatrix}; kwargs...)
     npoints = recurrencerate(R; kwargs...)*_rrdenominator(R; kwargs...)
     return _determinism(diagonalhistogram(R; kwargs...), npoints)
 end
@@ -197,7 +202,7 @@ end
 Calculate the divergence of the recurrence matrix `R`
 (actually the inverse of [`dl_max`](@ref)).
 """
-divergence(R::Union{ARM,SparseMatrixCSC}; kwargs...) = ( 1.0/dl_max(R; kwargs...) )
+divergence(R::Union{ARM,AbstractMatrix}; kwargs...) = ( 1.0/dl_max(R::Union{ARM,AbstractMatrix}; kwargs...) )
 
 """
     trend(R[; border=10, theiler])
@@ -242,7 +247,7 @@ 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::Union{ARM,SparseMatrixCSC}; theiler=deftheiler(R), kwargs...) =
+trend(R::Union{ARM,AbstractMatrix}; theiler=deftheiler(R), kwargs...) =
     _trend(tau_recurrence(R); theiler=theiler, kwargs...)
 
 """
@@ -251,7 +256,7 @@ trend(R::Union{ARM,SparseMatrixCSC}; theiler=deftheiler(R), kwargs...) =
 Compute the diagonalwise recurrence rate `τ-RR` from a recurrence matrix `R`.
 Note that this only works for squared recurrence matrices.
 """
-function tau_recurrence(R::Union{ARM,SparseMatrixCSC})
+function tau_recurrence(R::Union{ARM,AbstractMatrix})
     n = minimum(size(R))
     rv = rowvals(R)
     rr_τ = zeros(n)
@@ -321,7 +326,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::Union{ARM,SparseMatrixCSC}; kwargs...)
+function laminarity(R::Union{ARM,AbstractMatrix}; kwargs...)
     npoints = recurrencerate(R)*_rrdenominator(R; kwargs...)
     return _laminarity(verticalhistograms(R; kwargs...)[1], npoints)
 end
@@ -339,7 +344,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::Union{ARM,SparseMatrixCSC}; kwargs...) = vl_average(R; kwargs...)
+trappingtime(R::Union{ARM,AbstractMatrix}; kwargs...) = vl_average(R::Union{ARM,AbstractMatrix}; kwargs...)
 ###########################################################################################
 # 3. Based on recurrence times
 ###########################################################################################
@@ -357,7 +362,7 @@ default values and usage of the keyword argument `theiler`).
 
 Equivalent to [`rt_average`](@ref).
 """
-meanrecurrencetime(R::Union{ARM,SparseMatrixCSC}; kwargs...) = rt_average(R; kwargs...)
+meanrecurrencetime(R::Union{ARM,AbstractMatrix}; kwargs...) = rt_average(R::Union{ARM,AbstractMatrix}; kwargs...)
 
 """
     nmprt(R[; lmin=2, theiler])
@@ -378,7 +383,7 @@ of recurrence times [1].
 [DOI:10.1103/physreve.75.036222](https://doi.org/10.1103/physreve.75.036222)
 
 """
-nmprt(R::Union{ARM,SparseMatrixCSC}, kwargs...) = maximum(verticalhistograms(R;theiler=deftheiler(R), kwargs...)[2])
+nmprt(R::Union{ARM,AbstractMatrix}, kwargs...) = maximum(verticalhistograms(R;theiler=deftheiler(R), kwargs...)[2])
 ###########################################################################################
 # 4. All in one
 ###########################################################################################

src/rqa/windowed_function.jl

@@ -1,19 +1,17 @@
 """
-windowed(rmat, f, windth, step = 1; kwargs...)
+    windowed(rmat, f, width, step = 1; kwargs...)
 
 A convenience function that applies the RQA function `f`, such as `determinism`,
 to windowed views of the given recurrence matrix `rmat` with given window `width`
 and `step`. The `kwargs...` are propagated to the call `f(rmat_view; kwargs...)`.
 """
-function windowed(rmat, f, windth, step = 1; kwargs...)
-    windows = 1:step:(size(rmat, 1)-width)
-    map(1:length(windows)) do i
-        rmat_view = view(
-          rmat,
-          windows[i]:(windows[i]+width),
-          windows[i]:(windows[i]+width)
-        )
-        f(rmat_view; kwargs...)
-    end
-end
-
+function windowed(rmat::Union{ARM,AbstractMatrix}, f::Function, width::Integer, step=1::Integer; kwargs...)
+  (width < 2) && throw(ErrorException(
+        "Window width must be must be greater than or equal to 2"))
+  (step < 1) && throw(ErrorException(
+        "Step size must be greater than or equal to 1"))
+  windows = 1:step:(size(rmat, 1)-width+1)
+  map(1:length(windows)) do i
+    f(rmat[windows[i]:(windows[i]+width-1), windows[i]:(windows[i]+width-1)]; kwargs...)
+  end
+end

test/runtests.jl

@@ -12,4 +12,5 @@ defaultname(file) = uppercasefirst(replace(splitext(basename(file))[1], '_' => '
     testfile("rmatrix_analytic.jl")
     testfile("rqa_rna_analytic.jl")
     testfile("skeletontest.jl")
+    testfile("windowtest.jl")
 end

test/windowtest.jl

@@ -0,0 +1,43 @@
+using RecurrenceAnalysis, Test
+
+t = range(0, 2π; length = 300) # length is crucial and decides distance and thresholds
+c = cos.(t)
+s = sin.(t)
+X = StateSpaceSet(c, s)
+
+dmat = distancematrix(X)
+threshold = dmat[1, 2] + 0.01 # distance between two points
+
+@testset "window function" begin
+
+    rec = RecurrenceMatrix(X, threshold)
+
+    #classical
+
+    @test length(windowed(rec, recurrencerate, 30, 30)) == 10
+    @test length(windowed(rec, determinism, 30, 30)) == 10
+    @test length(windowed(rec, dl_average, 30, 30)) == 10
+    @test length(windowed(rec, dl_max, 30, 30)) == 10
+    @test length(windowed(rec, dl_entropy, 30, 30)) == 10
+    @test length(windowed(rec, divergence, 30, 30)) == 10
+    @test length(windowed(rec, trend, 30, 30)) == 10
+
+    #extended
+
+    @test length(windowed(rec, laminarity, 30, 30)) == 10
+    @test length(windowed(rec, trappingtime, 30, 30)) == 10
+    @test length(windowed(rec, vl_average, 30, 30)) == 10
+    @test length(windowed(rec, vl_max, 30, 30)) == 10
+    @test length(windowed(rec, vl_entropy, 30, 30)) == 10
+
+    #recurrence time
+
+    @test length(windowed(rec, meanrecurrencetime, 30, 30)) == 10
+    @test length(windowed(rec, nmprt, 30, 30)) == 10
+    @test length(windowed(rec, rt_entropy, 30, 30)) == 10
+    @test length(windowed(rec, rt_average, 30, 30)) == 10
+
+    #check if reacts to changes in width and step
+    @test length(windowed(rec, recurrencerate, 30, 1)) == 271
+    @test length(windowed(rec, recurrencerate, 10, 1)) == 291
+end