bugfixes of RQA functions and sorteddistances

Author: heliosdrm

.travis.yml

@@ -5,15 +5,14 @@ os:
   - osx
 julia:
   - 0.6
+  - 0.7
   - nightly
 branches:
   - only:
     - master
-    - julia_0.1_support
 matrix:
   allow_failures:
     - julia: nightly
-    - branches: julia_0.1_support
 notifications:
   email: false
 # uncomment the following lines to override the default test script

NEWS.md

@@ -1,5 +1,11 @@
 # RecurrenceAnalysis.jl News
 
+## 31-07-2018 - v0.2.1
+
+* Bugfixes:
+    - `sorteddistances`
+    - RQA functions (bugs in parameters of diagonal structures when `theiler==0`)
+
 ## 23-06-2018 - v0.2.0
 
 * Update to support Julia 0.6-0.7 (breaking compatibility with 0.5).

README.md

@@ -141,7 +141,7 @@ In some cases, specially with very long time series, it may be suitable to perfo
 ```julia
 @windowed determinism(rmat, theiler=2, lmin=3) width=1000 step=100
 ```
-will return a 91-element vector, such that each value is the determinism associated to a 1000-point fragment, starting at every 100 points (i.e. at `1`, `101`, &ldots; `9001`).
+will return a 91-element vector, such that each value is the determinism associated to a 1000-point fragment, starting at every 100 points (i.e. at `1`, `101`, … `9001`).
 
 The general syntax of that macro is:
 ```julia
@@ -154,7 +154,7 @@ where:
  * `w` is the width of the window for relevant data around each point.
  * `s` is the step or distance between points where the calculations are done (starting in the first point).
 
-To prevent syntax failures in the expansion of the macro, identify the RQA function (`rqa`, `recurrencerate`, `determinism`,&ldots;) directly by its name (avoid aliases), and use simple variable names (not complex expressions) for the arguments. On the other hand, the windowing options `w` and `s` can be given in any order. If `s` is ommitted, the calculations are done at every point, and the keyword `width` may be ommitted. (However, using `step=1` may be computationally very expensive, and that will provide just overly redundant results around each point, so it is advisable to set `step` a relatively big fraction of the window `width`.)
+To prevent syntax failures in the expansion of the macro, identify the RQA function (`rqa`, `recurrencerate`, `determinism`,…) directly by its name (avoid aliases), and use simple variable names (not complex expressions) for the arguments. On the other hand, the windowing options `w` and `s` can be given in any order. If `s` is ommitted, the calculations are done at every point, and the keyword `width` may be ommitted. (However, using `step=1` may be computationally very expensive, and that will provide just overly redundant results around each point, so it is advisable to set `step` a relatively big fraction of the window `width`.)
 
 The value returned by the macro will normally be a vector with the same type of numbers as expected by `expr`. In the case of `@windowed rqa(...) ...`, it will return a dictionary with a similar structure as in the default `rqa` function, but replacing scalar values by vectors.
 

src/radius.jl

@@ -24,7 +24,7 @@ The keyword arguments are the same that should be passed to the functions
 `rmat = recurrencematrix(x, d[i]; kwargs...)`, then
 `recurrencerate(rmat; kwargs...) == r[i]`.
 """
-function sorteddistances(x; theiler::Integer=0, scale=maximum, kwargs...)
+function sorteddistances(x; theiler::Integer=0, scale=1, kwargs...)
     # Create distance matrix
     kwargs = Dict(kwargs)
     argsdm = haskey(kwargs,:metric) ? (x, kwargs[:metric]) : (x,)
@@ -37,9 +37,9 @@ function sorteddistances(x; theiler::Integer=0, scale=maximum, kwargs...)
     ntheiler = theiler*n - theiler*(theiler-1)/2
     distarray = zeros(round(Integer,nd-ntheiler))
     pos = 0
-    for d = theiler:n-1
-        tmp = dm[n*d+1 : n+1 : n^2]
-        distarray[pos .+ (1:length(tmp))] = tmp
+    for d = 1:n
+        tmp = dm[d+theiler:n,d]
+        distarray[pos .+ (1:length(tmp))] .= tmp
         pos += length(tmp)
     end
     sort(distarray), (1.:length(distarray))/length(distarray)

src/rqa.jl

@@ -7,7 +7,7 @@ Calculate the recurrence rate (RR) of a recurrence matrix, ruling out
 the points within the Theiler window.
 """
 function recurrencerate(x::AbstractMatrix; theiler::Integer=0, kwargs...)
-    theiler < 0 && error("Theiler window length must be greater than 0")
+    theiler < 0 && error("Theiler window length must be greater than or equal to 0")
     if theiler == 0
         return typeof(0.0)( count(!iszero, x)/prod(size(x)) )
     end
@@ -23,13 +23,13 @@ end
 
 function tau_recurrence(x::AbstractMatrix{Bool})
     n = minimum(size(x))
-    [count(!iszero, diag(x,d))/(n-d) for d in (1:n-1)]
+    [count(!iszero, diag(x,d))/(n-d) for d in (0:n-1)]
 end
 
 # Based on diagonal lines
 
 function diagonalhistogram(x::AbstractMatrix{Bool}; theiler::Integer=1, kwargs...)
-    theiler < 0 && error("Theiler window length must be greater than 0")
+    theiler < 0 && error("Theiler window length must be greater than or equal to 0")
     bins = [0]
     nbins = 1
     current_diag = 0
@@ -70,19 +70,35 @@ function diagonalhistogram(x::AbstractMatrix{Bool}; theiler::Integer=1, kwargs..
 end
 
 function diagonalhistogram(x::SparseMatrixCSC{Bool}; theiler::Integer=1, kwargs...)
-    theiler < 0 && error("Theiler window length must be greater than 0")
+    theiler < 0 && error("Theiler window length must be greater than or equal to 0")
     m,n=size(x)
     rv = rowvals(x)
     dv = colvals(x) - rowvals(x)
+    loi_hist = Int[]
     if issymmetric(x)
-        valid = (dv .>= theiler)
+        valid = (dv .>= max(theiler,1))
         f = 2
+        # If theiler==0, the LOI is counted separately to avoid duplication
+        if theiler == 0
+            loi_hist = verticalhistogram(hcat(diag(x,0)))
+        end
     else
         valid = (abs.(dv) .>= theiler)
         f = 1
     end
     vmat = sparse(rv[valid], dv[valid] .+ (m+1), true)
-    f .* verticalhistogram(vmat)
+    dh = f .* verticalhistogram(vmat)
+    # Add frequencies of LOI if suitable
+    if (nbins_loi = length(loi_hist)) > 0
+        nbins = length(dh)
+        if nbins_loi > nbins
+            loi_hist[1:nbins] .+= dh
+            dh = loi_hist
+        else
+            dh[1:nbins_loi] .+= loi_hist
+        end
+    end
+    dh
 end
 
 """
@@ -93,7 +109,7 @@ the points within the Theiler window and diagonals shorter than a minimum value.
 """
 function determinism(diag_hist::Vector; lmin=2, kwargs...)
     if lmin < 2
-        error("lmin must be 2 or higher")
+        error("lmin must be 2 or greater")
     end
     nbins = length(diag_hist)
     diag_points = collect(1:nbins) .* diag_hist
@@ -112,7 +128,7 @@ the points within the Theiler window and diagonals shorter than a minimum value.
 """
 function avgdiag(diag_hist::Vector; lmin=2, kwargs...)
     if lmin < 2
-        error("lmin must be 2 or higher")
+        error("lmin must be 2 or greater")
     end
     nbins = length(diag_hist)
     diag_points = collect(1:nbins) .* diag_hist
@@ -148,7 +164,7 @@ the points within the Theiler window and diagonals shorter than a minimum value.
 """
 function entropy(diag_hist::Vector; lmin=2, kwargs...)
     if lmin < 2
-        error("lmin must be 2 or higher")
+        error("lmin must be 2 or greater")
     end
     nbins = length(diag_hist)
     if lmin <= nbins
@@ -173,9 +189,10 @@ the points within the Theiler window and in the outermost diagonals.
 function trend(npoints::Vector; theiler=1, border=10, kwargs...)
     nmax = length(npoints)
     rrk = npoints./collect(nmax:-1:1)
-    m = nmax-border
-    w = collect(theiler:m) .- m/2
-    typeof(0.0)( (w'*(rrk[theiler:m] .- mean(rrk[theiler:m])) ./ (w'*w))[1] )
+    a = 1+theiler
+    b = nmax-border
+    w = collect(a:b) .- b/2
+    typeof(0.0)( (w'*(rrk[a:b] .- mean(rrk[a:b])) ./ (w'*w))[1] )
 end
 
 function trend(x::AbstractMatrix; kwargs...)
@@ -294,6 +311,11 @@ maxvert(x::AbstractMatrix) = maxvert(verticalhistogram(x))
 
 """
     rqa(x; <keyword arguments>)
+
+Calculate RQA parameters of a recurrence matrix. See the functions
+`recurrencerate`, `determinism`, `avgdiag`, `maxdiag`, `divergence`, `entropy`,
+`trend`, `laminarity`, `trappingtime` and `maxvert` for the definition of
+the different parameters and the default values of the arguments.
 """
 function rqa(x; kwargs...)
     dhist = diagonalhistogram(x; kwargs...)