# Recurrence Plots

## Recurrence Matrices

A Recurrence plot (which refers to the plot of a matrix) is a way to quantify *recurrences* that occur in a trajectory. A recurrence happens when a trajectory visits the same neighborhood on the phase space that it was at some previous time.

The central structure used in these recurrences is the (cross-) recurrence matrix:

\[R[i, j] = \begin{cases} 1 \quad \text{if}\quad d(x[i], y[j]) \le \varepsilon\\ 0 \quad \text{else} \end{cases}\]

where $d(x[i], y[j])$ stands for the *distance* between trajectory $x$ at point $i$ and trajectory $y$ at point $j$. Both $x, y$ can be single timeseries, full trajectories or embedded timeseries (which are also trajectories).

If $x\equiv y$ then $R$ is called recurrence matrix, otherwise it is called cross-recurrence matrix. There is also the joint-recurrence variant, see below. With `RecurrenceAnalysis`

you can use the following functions to access these matrices

`RecurrenceAnalysis.RecurrenceMatrix`

— Type```
RecurrenceMatrix(x, ε; kwargs...)
RecurrenceMatrix{FAN}(x, ε; kwargs...)
```

Create a recurrence matrix from trajectory `x`

(either a `Dataset`

or a `Vector`

). Objects of type `<:AbstractRecurrenceMatrix`

are displayed as a `recurrenceplot`

.

**Description**

The recurrence matrix is a numeric representation of a "recurrence plot" [1, 2], in the form of a sparse square matrix of Boolean values.

`x`

must be a `Vector`

or an `AbstractDataset`

(possibly representing an embedded phase space; see `embed`

). If `d(x[i], x[j]) ≤ ε`

(with `d`

the distance function), then the cell `(i, j)`

of the matrix will have a `true`

value. The criteria to evaluate distances between data points are defined by the following keyword arguments:

`scale=1`

: a function of the distance matrix (see`distancematrix`

), or a fixed number, used to scale the value of`ε`

. Typical choices are`maximum`

or`mean`

, such that the threshold`ε`

is defined as a ratio of the maximum or the mean distance between data points, respectively (using`mean`

or`maximum`

calls specialized versions that are faster than the naive approach). Use`1`

to keep the distances unscaled (default).`fixedrate::Bool=false`

: a flag that indicates if`ε`

should be taken as a target fixed recurrence rate (see`recurrencerate`

). If`fixedrate`

is set to`true`

,`ε`

must be a value between 0 and 1, and`scale`

is ignored.`metric="euclidean"`

: metric of the distances, either`Metric`

or a string, as in`distancematrix`

.`parallel::Bool=false`

: whether to parallelize the computation of the recurrence matrix. This will split the computation of the matrix across multiple threads.

The parametrized constructor `RecurrenceMatrix{NeighborNumber}`

creates the recurrence matrix with a fixed number of neighbors for each point in the phase space, i.e. the number of recurrences is the same for all columns of the recurrence matrix. In such case, `ε`

is taken as the target fixed local recurrence rate, defined as a value between 0 and 1, and `scale`

and `fixedrate`

are ignored. This is often referred to in the literature as the method of "Fixed Amount of Nearest Neighbors" (or FAN for short); `RecurrenceMatrix{FAN}`

can be used as a convenient alias for `RecurrenceMatrix{NeighborNumber}`

.

If no parameter is specified, `RecurrenceMatrix`

returns a `RecurrenceMatrix{WithinRange}`

object, meaning that recurrences will be taken for pairs of data points whose distance is within the range determined by the input arguments. Note that while recurrence matrices with neighbors defined within a given range are always symmetric, those defined by a fixed amount of neighbors can be non-symmetric.

See also: `CrossRecurrenceMatrix`

, `JointRecurrenceMatrix`

and use `recurrenceplot`

to turn the result of these functions into a plottable format.

**References**

[1] : N. Marwan *et al.*, "Recurrence plots for the analysis of complex systems", *Phys. Reports 438*(5-6), 237-329 (2007). DOI:10.1016/j.physrep.2006.11.001

[2] : N. Marwan & C.L. Webber, "Mathematical and computational foundations of recurrence quantifications", in: Webber, C.L. & N. Marwan (eds.), *Recurrence Quantification Analysis. Theory and Best Practices*, Springer, pp. 3-43 (2015).

`RecurrenceAnalysis.CrossRecurrenceMatrix`

— Type```
CrossRecurrenceMatrix(x, y, ε; kwargs...)
CrossRecurrenceMatrix{FAN}(x, y, ε; kwargs...)
```

Create a cross recurrence matrix from trajectories `x`

and `y`

.

The cross recurrence matrix is a bivariate extension of the recurrence matrix. For the time series `x`

, `y`

, of length `n`

and `m`

, respectively, it is a sparse `n×m`

matrix of Boolean values, such that if `d(x[i], y[j]) ≤ ε`

, then the cell `(i, j)`

of the matrix will have a `true`

value.

Note that, unlike univariate recurrence matrices, cross recurrence matrices are not generally symmetric, regardless of the method used to make them.

See `RecurrenceMatrix`

for details, references and keywords. See also: `JointRecurrenceMatrix`

.

`RecurrenceAnalysis.JointRecurrenceMatrix`

— Type```
JointRecurrenceMatrix(x, y, ε; kwargs...)
JointRecurrenceMatrix{FAN}(x, y, ε; kwargs...)
```

Create a joint recurrence matrix from `x`

and `y`

.

The joint recurrence matrix considers the recurrences of the trajectories of `x`

and `y`

separately, and looks for points where both recur simultaneously. It is calculated by the element-wise multiplication of the recurrence matrices of `x`

and `y`

. If `x`

and `y`

are of different length, the recurrences are only calculated until the length of the shortest one.

See `RecurrenceMatrix`

for details, references and keywords. See also: `CrossRecurrenceMatrix`

.

`JointRecurrenceMatrix(R1, R2; kwargs...)`

Create a joint recurrence matrix from given recurrence matrices `R1, R2`

.

## Simple Recurrence Plots

The recurrence matrices are internally stored as sparse matrices with boolean values. Typically in the literature one does not "see" the matrices themselves but instead a plot of them (hence "Recurrence Plots"). By default, when a Recurrence Matrix is created we "show" a mini plot of it which is a text-based scatterplot.

Here is an example recurrence plot/matrix of a full trajectory of the Roessler system:

```
using DynamicalSystems
ro = Systems.roessler(ones(3), a=0.15, b=0.20, c=10.0)
N = 2000; Δt = 0.05
tr = trajectory(ro, N*Δt; Δt, Ttr = 10.0)
R = RecurrenceMatrix(tr, 5.0; metric = "euclidean")
recurrenceplot(R; ascii = true)
```

```
RecurrenceMatrix of size (2001, 2001) with 383597 entries
+------------------------------------------------------------+
| ..:'.::' .::'.:'..:' ..:'.:: .. .:'|
| .. .::'.:'' .:' .:'..''.::' . .::'.:' .:'..:' |
| .:':.:':::' . .:'.::'.::'.:' .:'.::'.::'.. ::'.::'.:' |
|:'.::'.::' :.:::.:' .. '.::'.::' :.:::.:' |
| ''' ''' ' .. .:''''' ..:'.. .:''''' ' .. ..:' ''' ..|
| .:'.::' .::'.:'.::' .::'.::' . .:''|
| .:' .:'..: ..:'..:'.:''.:' .:'..:'..: ..:'.::|
| ..:'.::'.::' .::'.::' :'.::' .::'.::'.:'' ' .::'.:'' |
|'.::'.::'.:' .:::.:'.. .:'' .:::.:::.:' .:':::' |
|:::::':::' '.::'.:'' ' ':::'.::' '.::'.. |
|'::' ::' .. .::'::' .. . .::'.::' .. .::'::' .|
| .:'.::' .::'.:' .:' ..:'.::' .::'|
| .::'.::'.: .:''.:'.::'.. .::'.:''.: . .:':.:|
| :'' :'.::'.: ::' ::'.:':.:'.. :' ::'.::'.:' ::'.::'.|
| . .:'..:' .:'.::'.:''.:' . ..:'.::' ..:'.::'|
| ..:'.::'.::' .::'.::'..'.::' ..:'.::'.:'' .::'.:'' |
|.::'.:::.:' :::.:'.::' ::' :::.:::.:' :':::'.:: |
|::::':::' .::'.::' . .::'.::' .::'.:'' |
| '' '' '.:' .:'' '' .:'..: .:'' ''' .:' .:'' '' .:|
| .:'..:'. .::'.:''.:'.. .:'.::'. . .::'.|
| .::'.::'.:' .:''.:'.::'.:' .::'.:':.:'. .:'..:'|
|. .:':.:'..:' '.::'.::' ''.:' ..:':::'..'' '.::'.::' |
|:::'.::'.:'' :::.::'.. ::' ::'.::'.:' :::.:'' |
|::::::::' .::'.::' .::::::' .::'.:: |
|::::::'.. ::.::' :::::'. ::.:' |
+------------------------------------------------------------+
```

`typeof(R)`

`RecurrenceMatrix{WithinRange}`

`summary(R)`

`"RecurrenceMatrix of size (2001, 2001) with 383597 entries"`

The above simple plotting functionality is possible through the package `UnicodePlots`

. The following function creates the plot:

`RecurrenceAnalysis.recurrenceplot`

— Function`recurrenceplot([io,] R; minh = 25, maxh = 0.5, ascii, kwargs...) -> u`

Create a text-based scatterplot representation of a recurrence matrix `R`

to be displayed in `io`

(by default `stdout`

) using `UnicodePlots`

. The matrix spans at minimum `minh`

rows and at maximum `maxh*displaysize(io)[1]`

(i.e. by default half the display). As we always try to plot in equal aspect ratio, if the width of the plot is even less, the minimum height is dictated by the width.

The keyword `ascii::Bool`

can ensure that all elements of the plot are ASCII characters (`true`

) or Unicode (`false`

).

The rest of the `kwargs`

are propagated into `UnicodePlots.scatterplot`

.

Notice that the accuracy of this function drops drastically for matrices whose size is significantly bigger than the width and height of the display (assuming each index of the matrix is one character).

Here is the same plot but using Unicode Braille characters

`recurrenceplot(R; ascii = false)`

```
⠀⠀⠀RecurrenceMatrix of size (2001, 2001) with 383597 entries⠀⠀
┌────────────────────────────────────────────────────────────┐
│⠀⠀⠀⠀⠀⠀⠀⠀⠀⢀⣠⠞⠃⢀⣴⠟⠁⠀⠀⠀⠀⠀⠀⠀⢀⣴⠟⠁⣠⡶⠋⢀⣤⠞⠁⠀⠀⠀⠀⠀⠀⠀⢀⣤⠞⠁⣀⡴⠟⠀⠀⠀⠀⢀⣤⠀⠀⣀⡴⠛│
│⠀⠀⠀⢀⡀⠀⠀⣀⣴⠟⠁⣠⡶⠋⠁⠀⠀⠀⣀⡴⠋⠀⣠⡾⠋⢀⣠⠘⠉⣀⡴⠟⠁⠀⠀⠀⢀⠀⠀⠀⣠⡶⠟⠁⣠⡾⠋⠀⠀⠀⠀⣠⡶⠋⢀⣠⡾⠋⠀⠀│
│⠀⣠⡶⠛⢁⣠⡾⠋⢁⣴⠾⠋⠀⡀⠀⠀⢠⡾⠋⢀⣴⠟⠋⣀⡴⠟⠁⢠⡾⠋⠀⠀⠀⠀⢠⡾⠋⢀⣴⡾⠋⢀⣴⠞⠉⢀⡀⠀⠀⣴⡾⠋⣀⣴⠟⠉⣀⡴⠋⠀│
│⡿⠋⣠⣴⡿⠋⣠⣶⠟⠁⠀⠀⠀⠀⠀⠀⢈⣠⣾⠟⢁⣠⡾⠋⠀⠀⠀⠀⢀⡀⠀⠀⠀⠀⠈⣠⣶⠟⠋⣠⣾⠟⠁⠀⠀⠀⠀⠀⠀⢁⣤⡾⠟⢁⣠⡾⠋⠀⠀⠀│
│⠀⠚⠛⠁⠀⠚⠛⠁⠀⠒⠀⣠⡄⠀⠀⣠⠞⠛⠁⠐⠛⠋⠀⠀⠀⢀⣠⠞⠋⢀⡄⠀⠀⣠⡼⠛⠁⠐⠛⠋⠀⠀⠒⠀⣠⡄⠀⢀⣤⠟⠋⠀⠐⠛⠉⠀⠀⠀⢀⣤│
│⠀⠀⠀⠀⠀⠀⠀⠀⠀⣤⠞⠉⢀⣴⠟⠁⠀⠀⠀⠀⠀⠀⠀⢀⣴⠟⠁⣠⡾⠋⢀⣴⠞⠉⠀⠀⠀⠀⠀⠀⠀⢀⣴⠞⠁⣀⡴⠟⠁⠀⠀⠀⠀⡀⠀⠀⣠⡴⠛⠁│
│⠀⠀⠀⠀⠀⠀⣀⡴⠂⠀⣠⡾⠋⢀⣠⠞⠀⠀⠀⠀⢀⣠⡾⠋⢀⣤⠞⠁⣠⡴⠛⠁⣠⡶⠃⠀⠀⠀⠀⣠⡴⠋⢀⣠⡾⠋⢀⣤⠞⠀⠀⠀⠀⢀⣤⠾⠋⢀⣴⠞│
│⠀⠀⠀⢀⣤⡾⠋⢀⣴⠟⠋⢀⡴⠟⠁⠀⠀⠀⣀⣴⠟⠁⣠⡴⠟⠁⠀⠸⠋⢀⣴⠞⠁⠀⠀⠀⢀⣴⠾⠋⣀⣴⠟⠁⣀⡴⠋⠁⠀⠒⠀⣠⣴⠟⠁⣠⡶⠋⠁⠀│
│⠃⣠⣶⠟⠉⣠⣾⠟⠁⣠⡾⠋⠀⠀⠀⠀⣠⣾⠟⢁⣤⡾⠋⢀⡀⠀⠀⢠⡴⠋⠁⠀⠀⠀⢠⣶⠟⢁⣠⡾⠟⢁⣠⠾⠋⠀⠀⠀⠀⣤⡾⠛⢁⣴⡾⠋⠀⠀⠀⠀│
│⣿⠟⣁⣴⡿⠋⣁⣴⠟⠋⠀⠀⠀⠀⠀⠀⠉⣠⣶⠿⠋⣠⡴⠛⠁⠀⠀⠈⠀⠀⠀⠀⠀⠀⠈⣡⣴⡿⠋⣠⣴⠟⠁⠀⠀⠀⠀⠀⠀⠉⣠⣾⠟⠁⣠⡄⠀⠀⠀⠀│
│⠁⠾⠿⠋⠀⠾⠟⠁⠀⠀⠀⢀⡀⠀⠀⣠⡾⠟⠁⠰⠾⠋⠀⠀⠀⠀⣠⡄⠀⢀⠀⠀⠀⣀⡼⠿⠋⠠⠾⠟⠁⠀⠀⠀⣀⡀⠀⠀⣠⠿⠟⠁⠰⠞⠋⠀⠀⠀⠀⣠│
│⠀⠀⠀⠀⠀⠀⠀⠀⠀⣠⡾⠋⢀⣴⠞⠉⠀⠀⠀⠀⠀⠀⠀⢀⣴⠞⠁⢠⡴⠋⠀⣠⡾⠋⠀⠀⠀⠀⠀⠀⠀⢀⣠⠞⠋⢀⣴⠟⠁⠀⠀⠀⠀⠀⠀⠀⢀⣴⠟⠁│
│⠀⠀⠀⠀⠀⠀⢀⣴⠞⠁⣠⡴⠟⠁⣠⠶⠀⠀⠀⠀⠀⣠⡴⠋⠁⣠⡾⠋⢀⣴⠟⠁⣀⡤⠀⠀⠀⠀⠀⣀⣴⠟⠁⣠⡶⠋⠁⣠⠖⠀⢀⠀⠀⠀⣠⡾⠋⢁⣠⠞│
│⠀⠀⠀⠀⠀⠶⠋⠁⠀⠾⠋⢀⣴⠞⠁⣠⡴⠀⠀⠰⠾⠋⠀⠰⠟⠁⣠⡼⠛⢁⣤⠞⠋⢀⡄⠀⠀⠀⠾⠋⠀⠰⠾⠋⣀⡴⠟⠁⣠⠾⠋⠀⠰⠟⠋⢀⡰⠟⠁⣠│
│⠀⠀⠀⠀⠀⢀⠀⠀⠀⣠⡾⠋⢀⣤⠞⠋⠀⠀⠀⠀⣠⡶⠃⢀⣴⠞⠋⢠⡴⠛⠁⣠⡾⠋⠀⠀⠀⠀⣀⠀⠀⢀⣠⡾⠋⢀⣴⠞⠁⠀⠀⠀⢀⣠⠞⠋⢀⣴⠞⠁│
│⠀⠀⢀⣠⡾⠋⢀⣴⠟⠋⢀⣴⠟⠁⠀⠀⠀⣀⣴⠟⠁⣠⣴⠟⠁⣠⠄⠈⢀⣴⠞⠉⠀⠀⠀⢀⣤⠞⠉⣀⣴⠟⠁⣠⡴⠛⠁⠀⠀⠀⣀⣴⠟⠁⣠⡶⠛⠁⠀⠀│
│⣠⣴⠟⠁⣠⣾⠟⢁⣠⡾⠋⠀⠀⠀⠀⠀⠸⠟⢁⣤⡾⠋⢀⣴⠞⠁⠀⠰⠟⠁⠀⠀⠀⠀⠰⠟⢁⣠⡾⠟⢁⣤⡾⠋⠀⠀⠀⠀⠀⠾⠋⢁⣴⡾⠋⢀⣴⠆⠀⠀│
│⠟⣁⣴⡿⠛⣁⣴⠿⠋⠀⠀⠀⠀⠀⠀⠀⢠⣶⡿⠋⣠⣶⠟⠁⠀⠀⠀⢀⠀⠀⠀⠀⠀⠀⢠⣴⡿⠋⣠⣴⠟⠉⠀⠀⠀⠀⠀⠀⠀⣠⣶⠟⠉⣠⡾⠋⠁⠀⠀⠀│
│⠀⠉⠉⠀⠀⠉⠁⠀⠀⠁⢀⡴⠂⠀⣠⡶⠋⠁⠀⠈⠉⠀⠀⠀⠀⣠⡶⠋⢀⣤⠆⠀⣀⡴⠋⠉⠀⠈⠉⠁⠀⠀⠀⣀⡴⠂⠀⣠⡶⠋⠁⠀⠈⠉⠀⠀⠀⠀⣠⠶│
│⠀⠀⠀⠀⠀⠀⠀⠀⣠⠾⠋⢀⣤⠞⠋⢀⠀⠀⠀⠀⠀⠀⢀⣴⠞⠉⣀⡰⠋⠁⣠⡾⠋⢀⡀⠀⠀⠀⠀⠀⠀⣠⠾⠋⢀⣴⠟⠁⢀⠀⠀⠀⠀⠄⠀⢀⣴⠟⠁⣀│
│⠀⠀⠀⠀⠀⢀⣴⠞⠁⣠⡴⠟⠁⣠⡶⠋⠀⠀⠀⠀⣠⡶⠛⠁⣠⡾⠋⢀⣴⠟⠁⣀⡴⠋⠀⠀⠀⠀⣀⣴⠟⠁⣠⡶⠋⢁⣠⠾⠋⢀⠀⠀⠀⣠⡾⠋⢀⣠⠾⠋│
│⡤⠀⠀⣠⡾⠛⢁⣤⡾⠋⢀⣤⠞⠁⠀⠀⠈⢀⣴⡾⠋⢀⣴⠟⠉⠀⠀⠈⠁⣠⠾⠋⠀⠀⠀⢀⣠⡾⠋⢁⣴⡾⠋⢀⡤⠈⠁⠀⠀⠉⢀⣴⠿⠋⣀⣴⠟⠁⠀⠀│
│⣁⣴⡿⠋⣠⣴⠟⠋⣠⡴⠛⠁⠀⠀⠀⠀⣼⠟⢉⣠⡾⠟⠁⣠⠄⠀⠀⢰⠟⠁⠀⠀⠀⠀⢰⡿⠋⣠⣶⠟⠁⣠⡶⠋⠀⠀⠀⠀⠀⣾⠟⢁⣠⡾⠋⠁⠀⠀⠀⠀│
│⡿⢋⣴⣾⠟⢁⣴⡾⠋⠀⠀⠀⠀⠀⠀⠀⢠⣴⡿⠋⣀⣴⠟⠁⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢀⣴⣿⠟⣁⣴⡿⠋⠀⠀⠀⠀⠀⠀⠀⣠⣴⡿⠋⣀⣴⠆⠀⠀⠀⠀│
│⣾⣿⢟⣡⣾⡿⠋⣀⡀⠀⠀⠀⠀⠀⠀⠀⠻⢋⣤⣾⠟⠁⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠸⢋⣥⣾⠟⠉⣀⠀⠀⠀⠀⠀⠀⠀⠀⠟⢋⣤⡾⠋⠀⠀⠀⠀⠀⠀│
└────────────────────────────────────────────────────────────┘
```

As you can see, the Unicode based plotting doesn't display nicely everywhere. It does display perfectly in e.g. VSCode, which is where it is the default printing type. Here is how it looks like in a dark background:

## Advanced Recurrence Plots

A text-based plot is cool, fast and simple. But often one needs the full resolution offered by the data of a recurrence matrix.

There are two more ways to plot a recurrence matrix using `RecurrenceAnalysis`

:

`RecurrenceAnalysis.coordinates`

— Function`coordinates(R) -> xs, ys`

Return the coordinates of the recurrence points of `R`

(in indices).

`RecurrenceAnalysis.grayscale`

— Function`grayscale(R [, bwcode]; width::Int, height::Int, exactsize=false)`

Transform the recurrence matrix `R`

into a full matrix suitable for plotting as a grayscale image. By default it returns a matrix with the same size as `R`

, but switched axes, containing "black" values in the cells that represent recurrent points, and "white" values in the empty cells and interpolating in-between for cases with both recurrent and empty cells, see below.

The numeric codes for black and white are given in a 2-element tuple as a second optional argument. Its default value is `(0.0, 1.0)`

, i.e. black is coded as `0.0`

(no brightness) and white as `1.0`

(full brightness). The type of the elements in the tuple defines the type of the returned matrix. This must be taken into account if, for instance, the image is coded as a matrix of integers corresponding to a grayscale; in such case the black and white codes must be given as numbers of the required integer type.

The keyword arguments `width`

and `height`

can be given to define a custom size of the image. If only one dimension is given, the other is automatically calculated. If both dimensions are given, by default they are adjusted to keep an aspect proportional to the original matrix, such that the returned matrix fits into a matrix of the given dimensions. This automatic adjustment can be disabled by passing the keyword argument `exactsize=true`

.

If the image has different dimensions than `R`

, the cells of `R`

are distributed in a grid with the size of the image, and a gray level between white and black is calculated for each element of the grid, proportional to the number of recurrent points contained in it. The levels of gray are coded as numbers of the same type as the black and white codes.

It is advised to use `width, height`

arguments for large matrices otherwise plots using functions like e.g. `imshow`

could be misleading.

For example, here is the representation of the above `R`

from the Roessler system using both plotting approaches:

```
using CairoMakie
fig = Figure(resolution = (1000,500))
ax = Axis(fig[1,1])
xs, ys = coordinates(R)
scatter!(ax, xs, ys; color = :black, markersize = 1)
ax.limits = ((1, size(R, 1)), (1, size(R, 2)));
ax.aspect = 1
ax2 = Axis(fig[1,2]; aspect = 1)
Rg = grayscale(R)
heatmap!(ax2, Rg; colormap = :grays)
fig
```

and here is exactly the same process, but using the embedded trajectory instead

```
y = tr[:, 2]
τ = estimate_delay(y, "mi_min")
m = embed(y, 3, τ)
R = RecurrenceMatrix(m, 5.0; metric = "euclidean")
xs, ys = coordinates(R)
fig, ax = scatter(xs, ys; markersize = 1)
ax.aspect = 1
fig
```

which justifies why recurrence plots are so fitting to be used in embedded timeseries.

It is easy when using `grayscale`

to not change the width/height parameters. The width and height are important when in `grayscale`

when the matrix size exceeds the display size! Most plotting libraries may resample arbitrarily or simply limit the displayed pixels, so one needs to be extra careful.

Besides graphical problems there are also other potential pitfalls dealing with the conceptual understanding and use of recurrence plots. All of these are summarized in the following paper which we suggest users to take a look at:

N. Marwan, *How to avoid potential pitfalls in recurrence plot based data analysis*, Int. J. of Bifurcations and Chaos (arXiv).

## Skeletonized Recurrence Plots

The finite size of a recurrence plot can cause border effects in the recurrence quantification-measures `rqa`

. Also the sampling rate of the data and the chosen recurrence threshold selection method (`fixed`

, `fixedrate`

, `FAN`

) plays a crucial role. They can cause the thickening of diagonal lines in the recurrence matrix. Both problems lead to biased line-based RQA-quantifiers and is discussed in:

K.H. Kraemer & N. Marwan, *Border effect corrections for diagonal line based recurrence quantification analysis measures*, Phys. Lett. A 2019.

`RecurrenceAnalysis.skeletonize`

— Function`skeletonize(R) → R_skel`

Skeletonize the `RecurrenceMatrix R`

by using the algorithm proposed by Kraemer & Marwan ^{[Kraemer2019]}. This function returns `R_skel`

, a recurrence matrix, which only consists of diagonal lines of "thickness" one.

Consider, e.g. a skeletonized version of a simple sinusoidal:

```
using DynamicalSystems, CairoMakie
data = sin.(2*π .* (0:400)./ 60)
Y = embed(data, 3, 15)
R = RecurrenceMatrix(Y, 0.25; fixedrate=true)
R_skel = skeletonize(R)
fig = Figure(resolution = (1000,600))
ax = Axis(fig[1,1]; title = "RP of monochromatic signal")
heatmap!(ax, grayscale(R); colormap = :grays)
ax = Axis(fig[1,2]; title = "skeletonized RP")
heatmap!(ax, grayscale(R_skel); colormap = :grays)
fig
```

This way spurious diagonal lines get removed from the recurrence matrix, which would otherwise effect the quantification based on these lines.

## Example

In the following we will plot recurrence plots of the Lorenz system for a periodic and chaotic regime (using scatter plot).

```
using DynamicalSystems, CairoMakie
lor = Systems.lorenz()
fig = Figure(resolution = (1000,600))
for (i, ρ) in enumerate((69.75, 28.0))
set_parameter!(lor, 2, ρ)
t, Δt = 20.0, 0.01
tr = trajectory(lor, t; Δt, Ttr = 2000.0)
tvec = 0:Δt:t
ax = Axis(fig[1, i]; title = "ρ = $ρ, " * (i != 1 ? "not periodic" : "periodic"))
lines!(ax, tr[:, 1], tr[:, 3]; color = Cycled(i), label = "X vs Z")
axislegend(ax)
ε = i == 1 ? 5.0 : 3.0
R = RecurrenceMatrix(tr, ε)
ax = Axis(fig[2, i])
ax.xlabel = "t"
i == 1 && (ax.ylabel = "t")
x, y = coordinates(R)
scatter!(ax, tvec[x], tvec[y]; markersize = 1, color = Cycled(i))
ax.limits = ((0, t), (0, t))
ax.aspect = 1
end
fig
```

On the left we see long (infinite) diagonals repeated over and over for different times. This is the case for periodic systems as they visit exactly the same area on the phase space again and again. The distance between the offset diagonals also coincides with the periodicity of the system, which is around `t ≈ 4`

.

On the right we see a structure typical of chaotic motion on a strange attractor such as the one of the Lorenz system: the orbit visits neighborhoods of previous points but then quickly diverges again. This results in many small diagonal lines.

## Distances

The distance function used in `RecurrenceMatrix`

and co. can be specified either as any `Metric`

instance from `Distances`

. In addition, the following function returns a matrix with the cross-distances across all points in one or two trajectories:

`RecurrenceAnalysis.distancematrix`

— Function`distancematrix(x [, y = x], metric = "euclidean")`

Create a matrix with the distances between each pair of points of the time series `x`

and `y`

using `metric`

.

The time series `x`

and `y`

can be `AbstractDataset`

s or vectors or matrices with data points in rows. The data point dimensions (or number of columns) must be the same for `x`

and `y`

. The returned value is a `n×m`

matrix, with `n`

being the length (or number of rows) of `x`

, and `m`

the length of `y`

.

The metric can be identified by a string, or any of the `Metric`

s defined in the `Distances`

package. The list of strings available to define the metric are:

`"max"`

or`"inf"`

for the maximum or L∞ norm (`Chebyshev()`

in the`Distances`

package).`"euclidean"`

for the L2 or Euclidean norm, used by default (`Euclidean()`

in`Distances`

).`"manhattan"`

,`"cityblock"`

,`"taxicab"`

or`"min"`

for the Manhattan or L1 norm (`Cityblock()`

in`Distances`

).

- Kraemer2019Kraemer, K.H., Marwan, N. (2019). Border effect corrections for diagonal line based recurrence quantification analysis measures. Physics Letters A 383(34).