# Measuring chaos in the Lorenz and Rössler models:
Fidelity tests for reservoir computing.

James Scully  Neuroscience Institute, Georgia State University, 100 Piedmont
Ave., Atlanta, GA 30303, USA.  Alexander Neiman  Department of Physics and
Astronomy, Ohio University, Athens, OH 45701, USA.  Andrey Shilnikov
Neuroscience Institute and Department of Mathematics & Statistics, Georgia
State University, 100 Piedmont Ave., Atlanta, GA 30303, USA.

###### Abstract

This study is focused on the qualitative and quantitative characterization of
chaotic systems with the use of symbolic description. We consider two famous
systems: Lorenz and Rössler models with their iconic attractors, and
demonstrate that with adequately chosen symbolic partition three measures of
complexity, such as the Shannon source entropy, the Lempel-Ziv complexity and
the Markov transition matrix, work remarkably well for characterizing the
degree of chaoticity, and precise detecting stability windows in the parameter
space.
The second message of this study is to showcase the utility of symbolic
dynamics with the introduction of a fidelity test for reservoir computing for
simulating the properties of the chaos in both models’ replicas. The results
of these measures are validated by the comparison approach based on one-
dimensional return maps and the complexity measures.

††preprint: AIP/123-QED

We employ the methods of qualitative theory and symbolic dynamics to measure
chaos and detect stability islands in one-parametric sweeps in the Lorenz and
Rössler models, as well as to compare chaotic properties in their reservoir
computed surrogates. We seek to test that reservoir computing algorithms are
able to pass all tests prepared for them: quantitive ones relaying on all
three measures extracted from binary sequences such as block entropy, Lempel-
Ziv complexity and Markov matrix structures, as well as qualitative ones based
on return maps. We hope that our algorithms and findings will be helpful for a
broad interdisciplinary audience including specialists and beginners in
dynamical systems and machine learning.

We dedicate this article to the memory of our colleague, teacher, and dear
friend – Professor Vadim Anischenko who made a fundamental contribution to the
field of nonlinear dynamics. Among so many aspects of nonlinear systems,
statistical properties of complex dynamics and methods of their
characterization were of his continuous interest.

##  I Introduction

Finding effective characterization of complex time series is a pivotal task
for the understanding of their underlying dynamics [1, 2, 3]. The Lorenz and
Rössler models are the classic examples of two types of deterministic chaos
observable in various low-dimensional dynamical systems [4], respectively,
with Lorenz-like attractors and spiral ones due to the Shilnikov saddle-focus.
A qualitatively different mechanism of formation and structure of chaos are
reflected in distinct statistical properties of these systems [5, 6, 7, 8]. As
such, the Lorenz and Rössler models serve as test-benches for testing and
development of new tools in the field of nonlinear science.

The first goal of this paper is to showcase how one may measure the degrees of
chaotic, homoclinic dynamics in such systems based upon the symbolic
description and how well the proposed approach complimentary agrees with the
conventional one employing the Lyapunov exponents. We also argue that the
symbolic approaches work exceptionally well to detect the stability windows in
parametric sweeps of such systems. The second goal of our paper is to
experiment to what degree machine learning tools such as reservoir computing
may well learn to reflect qualitatively and quantitatively on the dynamical
and probabilistic properties of original systems and surrogated ones.

While neither model needs to be introduced to the nonlinear community,
nevertheless let us first of all describe some of the key dynamical feathers
of both systems, as well as how their chaotic dynamics can be translated into
the symbolic description to generate long binary sequences to be further
quantified and analyzed using a simple technique of partitioning the phase
space or phase variables of these classic models.

The paper is organized as follows: first we discuss the Lorenz model, followed
by the Rössler model, and introduce the suitable symbolic description for
both. Next, we introduce the complexity measures using on the binary
framework; they include block-entropies and source entropy of symbolic
sequences, Lempel-Ziv (LZ) complexity. We argue that the entries of the Markov
transition matrix can effectively indicate structurally stable dynamics in
such pseudo-hyperbolic system, and therefore to detect the stability windows
in the parameter space. Next we compare these complexity measures with the
largest Lyapunov exponent as effective computational indicators of chaos and
periodic dynamics. Finally, we analyze and quantify the closeness of the
original chaotic dynamics occurring in the Lorenz and Rössler model and their
clones generated by various recurrent neural networks based on reservoir
computing principles.

##  II Models, symbolic partitions and binary framework

###  II.1 Lorenz model

Figure 1: (A) The Lorenz attractor at the classic value $r=29$. The
superimposed red dots defined by the $x$-variable critical events are well-
aligned on some straight-line intervals transverse to the wings of the Lorenz
butterfly in the phase space. (B) The $x$-variable plotted against time. Local
maxima and minima marked with red dots are detected to convert the
$x$-dynamics into binary sequences using the simple rule: $\\{\dot{x}=0\
|\,x>0\\}\to``0"$ and $\\{\dot{x}=0\ |\,x<0\\}\to``1"$ in this case. Figure 2:
(A) Chaotic transient converging to a stable attractor (green), encoded as
[$\overline{00001111}$], in the 3D phase space of the Lorenz model near a
stability window at $r\simeq$69.67. (B) The $x$-variable time traces passes
through a periodic pattern with the Markov transition probabilities
$p_{00}=3/4$ and $p_{01}=1/4$.

The Lorenz equation or model is given by:

$\dot{x}=-\sigma(x-y),\quad\dot{y}=r\,x-y+x\,z,\quad\dot{z}=bx+x\,y,$ (1)

with $x,y$, $z$ being the phase variables, and $\sigma,\ r,\ b>0$ being
bifurcation parameters; we will keep $\sigma=10$ and $b=8/3$ fixed through
this study. For $r>28$ the model starts exhibiting chaotic behavior associated
with an iconic butterfly-shaped strange attractor depicted in Figs. 1-3 below.
As this model is $\mathbb{Z}_{2}$-symmetric, i.e., it supports the group
symmetry $(x,y,z)\leftrightarrow(-x,-y,z)$. This is manifested in the shape
the Lorenz attractor shown in the projection in Fig. 1A, which is filled in
with flip-flopping patterns of a single solution of Eqs. (1). This is well
seen from Fig. 1B representing a typical evolution of the $x$-variable in
time. One can see from this and similar traces shown in Figs. 2B and 3B that
switching $x$-patterns changes with variations of the bifurcation
$r$-parameter. Specifically, one can see that the $x$-patterns in the last two
figures become periodic after some chaotic transients. These correspond to the
stable periodic orbits, shown in green, embedded in the chaotic attractors in
the phase space as shown in Figs. 2B and 3B. On the other hand, for other
$r$-parameters values, the Lorenz attractor remains chaotic as shown in Figs.
8A-9A and can be seen from time-progressions of the $z$-variable in Figs.
8B-9B. This property makes the Lorenz attractor pseudo-hyperbolic or a quasi-
attractor [9, 10, 11, 12, 13, 14]. Without going into details, the first means
that it constantly changes due to homoclinic bifurcations of its key
contributor – the saddle at the origin with two outgoing separatrices that
fill in the butterfly of the Lorenz attractor in the given phase space
projections. In contrast, a pseudo-hyperbolic attractor becomes a quasi-
attractor with homoclinic tangencies causing the emergence of stable periodic
orbits within it that remains such on $r$-parameter intervals, known as the
stability windows or islands.

In this study, we will demonstrate how time-progressions of both variables can
be used to compare and contrast the chaotic and periodic dynamics generated by
the Lorenz model and to detect stability islands in parameter space.

Let us first discuss the way symbolic, binary representations can be
introduced to describe the flip-flop dynamics of the Lorenz model. Figure 1
illustrates the concept: every turn of a phase trajectory around the right
"eye" in the wing of the butterfly, which is given by the condition
$x^{\prime}=0|~{}~{}x>0$, generates symbol "1" in the binary sequence.
Otherwise, a turn around the left eye, when $x^{\prime}=0|~{}~{}x<0$, adds "0"
symbol to the sequence. For example, the sequence $\\{10000..10111...\\}$ for
the $x$-variable progression shown in Fig.1B that corresponds to the right
separatrix of the saddle at the origin in the 3D phase space. Such a sequence
can be aperiodic/chaotic for the canonical parameter value $r=28$, used in all
textbook on nonlinear dynamics. For other values, the sequence can become
periodic with a repetitive block, e.g. [$\overline{00001111}$] of period 8,
after some short or long transient, as in the case depicted in Fig. 2. The
periodic orbit in Fig. 3 has a shorter periodic block [$\overline{000111}$] of
period 6, and so forth.

This outlines the method and ultimate goal of the symbolic approach: first,
one picks some typical, long steady-state trajectory of the model, and second,
examines the binary code/sequence extracted from the $x$-variable progression.
The question is how to determine efficiently whether the sequence is periodic
or aperiodic, i.e., chaotic and what is the degree of chaoticity/complexity?

Figure 3: (A) Convergence to the stable periodic orbit (green), encoded as
[$\overline{000111}$], after a long chaotic transient in the 3D phase space of
the Lorenz model at $r=92.5$. Superimposed red dots defined as critical events
$\\{x^{\prime}=0|\,x>0$ and $x<0\\}$ fill out two hooks on the bending wings
of the butterfly in the phase space. (B) The $x$-variable plotted against time
reveals the attracting periodic pattern with the Markov transition
probabilities $p_{11}=2/3$ and $p_{10}=1/3$.

###  II.2 Rössler model

Figure 4: (A) Long chaotic transient (grey) towards a stable periodic orbit
(green) in the 3D phase space of the Rössler model at $a=0.341$ and $c=4.8$.
Black dots indicate the location of $z^{\prime}=0$ events to generate binary
sequences $\\{...00011.,.\\}$ depending on where the critical events occur
below or above some $z$-threshold. (B) Spiking $z$-variable plotted against
time becomes regularized to produce a periodic pattern of low complexity.

The Rössler model [15, 16] is another classical example of deterministic chaos
occurring in many low-order systems. We use the following representation of
the Rössler model:

$\dot{x}=-y-z,\quad\dot{y}=x+ay,\quad\dot{z}=bx+z(x-c),$ (2)

with $x,y$, $z$ being the phase variables, and $a>0$ being bifurcation
parameters; here we keep $c=4.8$ and $b=0.3$ fixed. The convenience of the
representation (2) is that one equilibrium (EQ) state $O_{1}$ is always
located at the origin $(0,0,0)$, while the coordinates of the second one
$O_{2}$ are given by $(c-ab,b-c/a,-(b-c/a))$.

The best known feature of the Rössler model is the onset of chaotic dynamics
due to a Shilnikov saddle-focus [17, 18, 19], see Fig. 4. It begins with a
super-critical bifurcation of the stable equilibrium at the origin, followed
by a period-doubling bifurcation cascade as $a$ increases.

Recall that a 3D dissipative system with the Shilnikov saddle-focus cannot
produce a genuinely chaotic attractor but produces a quasi-attractor [10, 14]
instead. Homoclinic tangencies (literately stirred by two saddle-foci in the
Rössler model) inside such a quasi-attractor cause the emergence of stable
periodic orbits in the phase-space through saddle-node bifurcations followed
by period-doubling ones. Note that the chaotic attractor in the Lorenz model
is categorized as a pseudo-hyperbolic one according to [20, 21, 22, 23]; at
larger $r$-values it becomes a quasi-attractor [13] because of the presence of
various stability islands seen in Fig. 6. The same is true for the Rössler
model as illustrated by Fig. 7.

This observation lets us introduce the partition using critical events when
the $z$-variable reaches its maximal values on the attractor. The binary
sequence $\\{k_{n}\\}$ representing a trajectory is computed as follows:

$k_{n}=\begin{dcases*}1,&if~{}~{} $z_{max}$~{}~{}~{}~{}~{} > $z_{\rm th}$,\\\
0,&if~{}~{} $z_{min/max}$ $\leq$ $z_{\rm th}$,\end{dcases*},$ (3)

where the $z$-threshold can be set relative to the location of the secondary
equilibrium state $O_{2}$: $z_{\rm th}=0.1\,(c-ab)/a$ as was done in Ref.[19]
for example, or even set it fixed $z_{\rm th}=0.03$ or $z_{\rm th}=1.0$ as
done in this study. The choice of partition is motivated by its simplicity and
may differ if it satisfies the purpose, namely, to distinguish local dynamics
concentrated around the Shilnikov saddle-focus at the origin from large global
pathways associated with possible homoclinic excursions in such spiral
attractors, see Fig. 4 With this simple algorithm based on a threshold level
$z_{\rm th}$, we can convert the maximal values of the $z$-variable into the
binary framework: using $0$ when a solution of the model turns around the
saddle-focus, and $1$ when it transitions towards the other equilibrium state
$O_{2}$ and back to the origin. Figure 4A illustrates the concept in the 3D
phase space where the black dots indicate the critical events on the
attractors, while Fig. 4B represents the time-progression of the $z$-phase
variable.

##  III Methods: Complexity measures

There are several approaches available to assess such a degree. Perhaps the
easiest implementation is to convert a time-progression into a binary
sequence, apply a compression application, such as gzip, and then to compare
the lengths of the original, $L_{\text{orig}}$, and compressed,
$L_{\text{comp}}$, files. The compressibility measure can be introduced as
$R=1-L_{\text{comp}}/L_{\text{orig}}$ and it is related to the redundancy of
information contained in the sequence. As two benchmarks one can choose: (i)
non-redundant random Bernoulli sequence with the least compressibility, and
(ii) a redundant periodic sequence with the maximal compressibility. A
sequence generated by the deterministic chaotic system lies in between of
these two benchmark limits.

In what follows we will employ four different measures to examine the degree
of chaos in the Lorenz and Rössler models and to determine the stability
windows as the parameter is swept within some ranges of interest for these
systems. We state from the very beginning, that period-doubling bifurcations
are beyond the scope of this examination, and as such our partition designs
are not meant to detect such transitions. This can be obviously refined and
resolved with additional constrains on the chosen partitions.

The proposed chaos measures are the source entropy (SE), the Lempel-Ziv
complexity (LZ), the Markov transition matrix, and the largest Lyapunov
exponent (LE). As the Reader will see, all of these measures work quite well
to detect the stability islands in the chaos sea. The key to understanding why
this is the case is rooted in the fact that the Lorenz attractor is pseudo-
hyperbolic and then all four measures change abruptly with parameter
variations, except for stability intervals where they remain constant.
Specifically, SE, LE, and LZ vanish after converging to exponentially stable,
and therefore structurally stable, periodic orbits.

###  III.1 Block-entropies and the source entropy of symbolic sequences

A generic measure of complexity of a symbolic sequence is provided by the
entropy of the source [24]. Given a binary sequence $S$, the Shannon’s entropy
of words of $m$-symbols long, the so-called block entropy, is defined as [25,
26]

$H_{m}=-\sum_{\\{s_{m}\\}}P(s_{m})\log P(s_{m}),$ (4)

where $P(s_{m})$ is the probability of occurrence of a word of length $m$
within the sequence $\S$, and summation is carried over all words of length
$m$ occurring with nonzero probability. The $m$-block entropy, $H_{m},$ is
interpreted as average information contained in a word of length $m$. The
conditional entropies are defined as in reference[26],

$h_{m}=H_{m+1}-H_{m},\quad h_{0}:=H_{1},$ (5)

and provide the average information required for prediction of $(m+1)$ symbol,
given that the preceding $m$ symbols are known. The limit of $m\to\infty$
gives the quantity of interest, entropy of the source,

$h=\lim_{m\to\infty}h_{m}=\lim_{m\to\infty}\cfrac{H_{m}}{m}.$ (6)

For a dynamical system, the source entropy depends on a particular partition
of the phase space, that is, on a rule that maps the evolution of the
dynamical system to a discrete symbolic sequence. The Metric or Kolmogorov-
Sinai entropy is an upper bound of the source entropy over all finite
partitions; the entropies match for the so-called generating partition [27,
28].

For a symbolic sequence generated by a periodic source, such as a stable
periodic orbit of the Lorenz model, the source entropy equates to 0, while the
conditional entropy, $h_{m}$, drops to 0 when the word length reaches the
period, $m=p$, and so $h_{p}=h_{p+1}=...=0$. That is, it is enough to observe
a periodic sequence for just a period to predict with certainty the next
symbol. Equivalently, no new information is gained after one period of
periodic sequence is observed.

For a random or chaotic sequence the source entropy is positive [25, 29],
reflecting the sole fact of uncertainty in prediction of next symbol in
chaotic sequence even if the entire prehistory is known. The decay of the
conditional entropy, $h_{m}$, to its asymptotic value $h$ is a generic measure
of correlations in the sequence [25, 30]. In particular, for a Markov sequence
of memory $p$ the conditional entropy converges to the source entropy after
exactly $p$ steps, i.e. $h_{p}=h_{p+1}=...=h$.

The number of words grows quickly with the word length. For example, in the
case of Bernoulli binary sequence the number of possible words of length $m$
is $M_{m}=2^{m}$. In a numerical experiment the length of a sequence generated
by the dynamical system is always finite, which creates a well-known problem
in estimation of long-words probabilities [31, 30]. That is, some words cannot
be observed or encountered only few times, simply because the sequence is not
long enough. In result, the block entropies, $H_{m}$, are systematically
underestimated and finite-size correction must be applied. As we are
interested in an indicator of chaotic or periodic dynamics, in the following
we limit the word length to $m=6$ and use $h_{6}$ for estimation of the source
entropy. We collect long sequences, $N\gg 2^{6}$, and use the finite size
correction [30],

$H^{\text{observed}}_{m}\approx H_{m}-\frac{M_{m}-1}{2N},$ (7)

where $H^{\text{observed}}_{m}$ is the observed $m$-block entropy calculated
from the observed sequence of length $N$, $H_{m}$ is the true entropy, $M_{m}$
is the number of distinct $m$-words, $M_{m}\leq 2^{m}$.

###  III.2 The Lempel-Ziv complexity

Given a binary sequence, the Lempel-Ziv complexity (LZ) is related to the
number of substrings of increasing lengths that the given sequence is made of.
For example, for the sequence, $\\{1|0|10|010|0101|11|110\\}$, scanned from
left to right, that number $s$ is 8. Then $LZ$-complexity can be defined as

$LZ=s\log(N)/N,$ (8)

where $N$ is the length of the sequence. Indeed, the LZ algorithm can be used
as an effective estimator of the source entropy [32]. Thus, we expect the
source entropy estimate, $h_{6}$, and the LZ to be strongly correlated, as the
parameter of a dynamical system varies.

###  III.3 The Markov transition matrix

The simplest measure is introduced when it is based on elements of the
transition matrix of the single-step Markov process underlying the symbolic
sequence,

$\textbf{M}=\left[\begin{matrix}p_{11}&p_{01}\\\
p_{10}&p_{00}\end{matrix}\right],$

where $p_{11}$ and $p_{10}$ are probabilities of flop and flip, respectively,
i.e., of transitions $1\to 1$ and $1\to 0$. Because the Lorenz attractor is
symmetric, $p_{11}\simeq p_{00}$, which implies that $p_{10}\simeq p_{01}$ as
well. The situation is different when the attractor is a stable asymmetric
periodic orbit, or either one in a pair of asymmetric chaotic attractors that
emerged through a period-doubling cascade. The transition probabilities
$p_{11}$ and $p_{01}$ can be compared with $0.5$, e.g. with Bernoulli trial of
an unbiased coin. For example, the Markov matrices for the symmetric stable
periodic orbits shown in Figs. 2 and 3 are, resp., the following

$\left[\begin{matrix}4/5&1/5\\\
1/5&4/5\end{matrix}\right]\quad\mbox{and}\quad\left[\begin{matrix}3/4&1/4\\\
1/4&3/4\end{matrix}\right],$

and will remain such within the corresponding stability windows in the
parameter space.

###  III.4 The largest Lyapunov exponent

Figure 5: Conditional entropy versus the words length for the Lorenz model.
(A) Conditional entropy for chaotic sequences with $N=16988$ for $r=30$, and
$N=28475$ for $r=75$. (B) Conditional entropy for periodic sequences with
length $N=25672$ for $r=59.25$, and $N=32058$ for $r=92.5$. Open circles
refers to the conditional entropy $h_{m}$ estimated from a sequence generated
by the genuine Lorenz model. Filled black circles shows $h_{m}$ estimated from
sequences generated by the trained reservoir computer whith the same lengths
of binary sequences.

A Lyapunov exponent is meant to indicate how quickly nearby trajectories may
converge/diverge in the phase space. The sum of positive Lyapunov exponents is
related to the Kolmogorov-Sinai (KS) entropy via Ruelle’s inequality,
$\text{KS}\leq\sum_{\Lambda_{i}>0}\Lambda_{i}$ [33]. Since the Lorenz system
is a strongly dissipative system that may have a single positive Lyapunov
exponent on the chaotic attractor, then the largest Lyapunov exponent (LE) is
directly related to the KS entropy, $\text{KS}\leq\Lambda$. Thus, the measures
introduced for symbolic dynamics can be compared against the LE,

$\Lambda=\lim_{t\to\infty}\frac{1}{t}\log\frac{|\delta\textbf{x}(t)|}{|\delta\textbf{x}_{0}|},$
(9)

that measures average rate of convergence or divergence between two
trajectories. So, $\Lambda<0$ means a trajectory converges to a stable
equilibrium state; $\Lambda=0$ corresponds to the case of a stable periodic
orbit along which the distance between two solutions does not change over its
period. The case $\Lambda>0$ indicate that solutions of the system under
consideration run on some chaotic attractor.

For a fair comparison with symbolic-sequence measures, such as the source
entropy, the LE should be normalized using a characteristic time of the
system, $\tau$, so that the quantity $\lambda=\Lambda\tau$ becomes
dimensionless. For the Lorenz attractor case, as the characteristic time we
use the mean dwelling time intervals between the events $\dot{x}=0$, such as
ones marked as the red dots on the butterfly wings in Fig. 1. One can observe
from Eqs. (1) that increasing the parameter $r$ speeds up the time derivative
$\dot{y}$ and hence $\dot{x}$, which results in shortening of the dwelling
time intervals between transition events which is compensated by the growing
size of the Lorenz attractor, compare the coordinate axis in Figs. 1 and 9,
for example. Similarly, for the Rössler model below we used the dwelling times
between similar critical events.

Figure 6: Complexity measures vs the parameter $r$ of the Lorenz model. For
each parameter $r$ value, symbolic sequences of the length $10^{4}$ were
collected. A: The largest Lyapunov exponent (LE, blue), the source entropy
estimate, $h_{6}$ (Entropy, orange), and Lempel-Ziv complexity (LZ, green) vs
the parameter $r$. The largest Lyapunov exponent was normalized to the average
inter-symbol interval, $\lambda T$. B: Markov probabilities $p_{11}$ and
$p_{01}$ show sharp peeks and plateaus that are also indicative, resp., of
narrow and wide stability windows with periodic orbits within.

##  IV Data and Results

Figure 5 presents the conditional entropies for chaotic and periodic regimes
of the Lorenz model. The first examples in Fig. 5A illustrate the practicality
of the conditional entropy in distinguishing two chaotic regimes: larger
$h_{m}$ value for the aperiodic binary sequence generated at $r=30$ compared
to that at $r=75$ is a direct indicator of the higher degree of uncertainty of
symbols predictions in the former. This implies that one can treat chaos at
nearly canonical parameter value is superior than that at $r=75$, see Fig 9A.
One indirect justification for that is the Lorenz model at $r\geq 31$ no
longer exhibit genuinely chaotic attractors but quasi-attractors with
homoclinic tangencies causing the emergence of stable periodic orbits such as
one occurring at $r=68.75$ and shown in Fig. 2. This should imply that binary
sequences generated by chaotic solutions of the Lorenz model at this parameter
range include multiple “laminar” or periodic substrings leading to lower
values of conditional entropy. We will proceed with more direct arguments in
the next section.

As for stable periodic orbits generating stable periodic orbits after some
chaotic transients, the conditional entropy quickly vanishes when the word
length reaches the period, see Fig. 5B. The distinction between two periodic
regimes is clearly captured by the entropy. First, the obvious observation is
that the stable periodic orbit observed at $r=59.25$ has a longer period than
the orbit found at $r=92.5$. A less obvious observation is related to the
value of $H_{1}$ (or $h_{0}$), that is the entropy of words of length 1. For
$r=92.5$ $H_{1}$ is given by $\log 2$, reflecting that probabilities of "0"
and "1" are the same. This indeed reflects the symmetry of underlying stable
periodic orbit. On the contrary, for $r=59.25$, the entropy $H_{1}$ is
significantly smaller than $\log 2$, reflecting the asymmetry of the
underlying periodic orbit.

Figure 7: Complexity measures vs the parameter $a$ of the Rössler model. For
each parameter $a$ value, symbolic sequences of the length $10^{4}$ were
collected. A: The largest Lyapunov exponent (LE, blue), the source entropy
estimate, $h_{6}$ (Entropy, orange), and Lempel-Ziv complexity (LZ, green) vs
the parameter $a$. The largest Lyapunov exponent was normalized to the average
inter-symbol interval, $\lambda T$. B: Markov probabilities $p_{00}$ and
$p_{01}$. Sharp dropping peeks and plateaus are indicative of narrow and wide
stability windows with periodic orbits within.

Figures 6 and 7 with several charts are the graphical culmination of our
simulations in the first part of this study dealing with complexity measure of
dynamics demonstrated by the Lorenz and Rössler models. As we said above that
leave the chaos due to period-doubling bifurcations apart from our
consideration, as follows directly from the choice of partitions in either
case.

So, Figure 6 represents a $r$-parametric sweep of dynamics of the Lorenz model
starting with the classical parameter value $r=28$: its panel B depicts
empirically the way the complexity measures vary and all correlate with
$r$-parameter increase. The low panel B in Fig. 6 shows variations of two
(asymmetric) entries, $p_{11}$ and $p_{01}$, within the same $r$-range. Note
that here $p_{11}$ and $p_{01}$ taken from different columns of the transition
matrix and hence do not add up to $1$ always, are used to detect stability
windows with paired mirror-asymmetric periodic orbits within.

Let us first observe the sharp peeks and plaques in these sweeps where all
three measures, LE, SE, and LZ drop down to zero. One can conclude that such
$r$-intervals correspond to narrow or wider stability windows within which the
Lorenz model will eventually demonstrate periodic dynamics. The corresponding
periodic orbit can be symmetric like those shown in Figs. 2 and 3. By
analyzing the values of $p_{11}$ and $p_{01}$-probabilities from Fig. 6B, one
can deduct, without visual inspection of Fig. 3 that the plateau around
$r=92.5$ correspond to the symmetric stable orbit, while sharp peeks near
$r=100$ are indicative of the stability windows with a pair of co-existing
stable orbits within, to either one the transient converges eventually.

One can see from Fig. 6A that at the initial stage prior to two first
stability windows around $r=59$, the values of the SE and LE agree quite well
with each other, while after the LE-curve starts diverging from the LZ and SE
curves. As we pointed out above, increasing $r$-parameter increases on average
the time-derivatives of the $x$ and $y$-phase variables, accelerating the time
course in the Lorenz model. While we used the average time interval between
the critical events at assigning binary symbols as the normalizing factor for
comparison of LE and SE, the distinction is still noticeably growing with
$r$-parameter increase. A possible reason can be that at larger $r$-values the
dynamics of the Lorenz model is no longer flip-flop due to homoclinic
bifurcations of the saddle at the origin but becomes mostly determined by
period doubling bifurcations that our binary framework is not designed to
account. In other words, the used binary partition is far from a generating
partition. Such regimes will have inclusions of longer, self-similarly
regularized blocks like $....111000111000...$ that result in lower entropy-
and LZ-based measures, whereas the average rate of divergence of two nearby
trajectories in the phase space remains basically intact.

Figure 8: (A) The Lorenz attractor at $r=28$. Superimposed black dots are the
$z$-critical events well-aligned on two line segments transverse to the
butterfly wings in the phase space. (B) The $z$-variable plotted against time
is superimposed with local maxima given by
$\\{z^{\prime}=0|\,z^{\prime\prime}<0\\}$ to generate the Poincaré return maps
shown in (C). (C) 1D maps $T:z_{n}\to z_{n+1}$ generated by sequential pairs
of local maxima of the $z$-variable in the Lorenz model (black dots), while
the map filled out with red dots are generated from the ML-emulated trace
(red) shown in (D). The Lorenz $z$-map has a characteristic cusp-shape (due to
high outbursts of the separatrix of the saddle, see (A)) and an emerged fold
at $r=75$ giving rise to stable orbits due to tangent (SN) bifurcations. The
emulated map adequately captures some characteristics of the original on
average.

Figure 6A suggests non-ambiguously that chaos in the Lorenz model at low
$r$-values is more unpredictable and homogenous due to shorter flip-flop
dynamics caused by homoclinic bifurcations of the saddles than at higher
parameter values where the dynamics is mostly dominated by longer recurrent
patterns due to period-doubling bifurcations. Figure 6B supports this
assertion that the complexity is maximized till the entries of the transition
matrix remains close to 1/2 on the left in the sweep before the very first
stability islands.

One-parametric sweeps of the largest Lyapunov exponent, the LZ-complexity, and
the source entropy for the Rössler model are presented in Fig. 7A, while Fig.
7B depicts the way the Markov transition probabilities change accordingly with
variations of the $a$-parameter. One can see from its initial phase $a<0.33$
that complexity measures stay close to zero, while the LE is indicative of the
origin of chaotic dynamics alternating with periodic dynamics in the phase
space of the Rössler model. The examination of Fig. 7A reveals that the
probability $p_{00}=1$ till $a<0.33$, which indicates by the phase space
partition that the dynamics of the Rössler model remains “flat” near the
origin (see Fig. 4) that became a saddle-focus after a super-critical
Andronov-Hopf bifurcation followed by period-doubling cascades through which
the dynamics becomes weekly chaotic as indicated by the positive largest
Lyapunov exponents.

For $a>0.33$, the Shilnikov saddle-focus starts contributing to more
developing chaotic dynamics when its 2D stable manifold bends so that
trajectories come closer to the proximity of the origin. With large $a$-values
the chaotic attractor appears to look more like iconic spiral chaos in the
Rössler model. The growing values of the complexity measures along with the
$p_{01}$-probability in the transition matrix approaching 1/2 from below are
all synchronous indicators that the degree of chaos and unpredictability on
the Rössler model increases as it transitions further from period-doubling
type towards spiral and multi-funnel chaos with countably many homoclinic
orbits originated by and stirred by both Shilnikov saddle-foci in the model,
see Ref.[19] for more details. As was said earlier, in addition to chaos, the
Shilnikov homoclinic bifurcations lead to the emergence of homoclinic
tangencies, and hence the abundance of saddle-node bifurcations that determine
the borders of multiple stability windows in the parameter space that are
populated by sequential period-doubling cascades within. This is well-seen in
the sweep in Fig. 7. The reader is welcome to consult with Refs. [34, 35, 36,
37] that elucidate computationally the contribution of the Shilnikov saddle-
focus bifurcation in the formation of the so-called T-points for heteroclinic
connections between saddles and saddle-foci in various Lorenz-like systems and
the role of inclination-flip bifurcations in such unfolding that amplify
homoclinic tangencies in producing Shilnikov flames – stability islands nearby
such T-points.

To conclude this section we re-iterate that the conditional entropy captures
well the parameter dependence of the dynamics, showing qualitatively the same
dependence as the largest Lyapunov exponent, and LZ-complexity as de-facto
depicted by the parametric sweeps of the Lorenz and Rössler model. Let us
remind that for a fair comparison in the largest Lyapunov exponent was
normalized by the average time intervals between the critical events used in
the binary framework design. Finally, the transition probabilities cab be
effectively used as indicators to detect stability windows in parametric
sweeps as well.

##  V Reservoir-computing

Figure 9: (A) Chaos in the Lorenz model is shown at $r=75$. Superimposed black
dots are $z$-variable critical events form two 1D hooks on the bending
butterfly wings in the phase space. (B) The $z$-variable plotted against time
is superimposed with local maxima given by
$\\{z^{\prime}=0|\,z^{\prime\prime}<0\\}$ to generate the Poincaré return maps
shown in (C). (C) 1D maps $T:z_{n}\to z_{n+1}$ generated by sequential pairs
of local maxima of the $z$-variable in the Lorenz model (black dots) and one
(red dots) from the ML-emulated trace (red) shown in (D). The Lorenz $z$-map
(no longer 1D due to loss of strong contraction) has a characteristic cusp-
shape (due to high outburst of the separatrix of the saddle (see (A)) and an
emerging fold at $r=75$ to give rise to stable orbits due to tangent (saddle-
node) bifurcations. The emulated map captures adequately some characteristics
of the original on average.

Symbolic dynamics can be used to study the properties of systems
identification techniques for chaotic systems. Echo State Networks (ESNs), a
form of reservoir computing (RC), are one such method for systems
identification that has recently become popular in the nonlinear dynamics
community [38, 39]. In the following section we apply a fidelity test to echo
state networks, trained for time series prediction. The fidelity test includes
two components: a quantitative test through complexity measures, specifically
conditional entropies of various block sizes LZ complexity, and a qualitative
test through canonical 1D Poincaré return maps : $T$: $z_{max}^{n}\to
z_{max}^{n+1}$.

In an ESN, the input is embedded into the reservoir’s vector space by a random
linear transformation. The resulting input and reservoir vectors are summed
and passed through some generic point-wise nonlinearity, hyperbolic tangent in
this case. This is derived from the discretization of a nonlinear leaky
integrator, yielding the update relation for the reservoir [38]:

$R_{n+1}=(1-\alpha)R_{n}+\alpha\tanh(W_{R}R_{n}+W_{in}y_{n}),$ (10)

where $R_{n}$ is the reservoir and $y_{n}$ is either the teacher input, or the
last prediction: $W_{out}R_{n}$.

The input reservoir and input matrices are $W_{R}$ and $W_{in}$, respectively.
The final prediction is obtained from $R_{n}$ as

$x_{n}=W_{out}R_{n},$ (11)

where $W_{out}$ is the matrix representation of the best fitting linear
transformation of the reservoir sequence, found by ridge regression. Detailed
methods for both Lorenz and Rössler models can be found in the following
sections. The results are summarized for each and followed by a brief
discussion below.

###  V.1 Lorenz model’s replica

The results of the quantitative test are summarized in Fig. 5. Shown in
circles are the values of the conditional entropy $h_{m}$ estimated from
binary sequences generated by the Lorenz model, while the black dots
represents the $h_{m}$-values estimated from binary sequences generated by the
trained reservoir computer.

The results of the qualitative test are illustrated in the sequence of two
Figs. 8, and 9. Each figure includes four panels: two panels on the left and
right depict the phase space projections of the original Lorenz attractor
(grey) and its RC reconstruction (purple). The black/red dots in these panels
indicate the critical events when the $z$-variables become maximized. The
bottom panel show the time progression of the original $z$-variable overlapped
with that of the surrogate, as well as the $z_{max}$ map in the middle panel.
This is a standard approach that is also applicable to RC surrogate systems
[40, 41].

Let us begin with the return maps shown Fig. 8B. One can see that both maps,
ordinal and surrogate, produce the expected cusp-shape at $r=28$ with the
slope being steeper then $\pm 1$ by comparing the graph with the bisectrix or
the $45^{\circ}$-line presented in this figure as well, i.e., the map is an
expansion. One can see that the cusp graph populated by the $z{max}$-values
(black dots) extracted from the Lorenz model is exceptionally well fitted by
those (red dots) from the RC surrogate.

Such a cusp map constitutes a de-facto computational proof that the Lorenz
attractor has no stable orbits and is made only of countable many
unstable/saddle periodic orbits whose 2D stable and unstable manifolds cross
transversely in the 3D phase space of the model. This assertion follows also
indirectly from the visual inspection of the dots, given by the condition
$x^{\prime}_{max/min}=0$, that fill out two straight line segments on the
opposite sides of the symmetric wings of the Lorenz attractor shown in Fig. 8A
while its RC surrogate is shown in Fig. 8C.

Figure 10: True solution (A) and predicted solutions (C) for Rössler attractor
are shown with orange dots demarcating local maxima with $z>1$ and blue dots
demarcating local minima. Accompanying time series are given in (D) and (E)
respectively. (B) shows the characteristic $z_{max}$-return maps. Figure 11:
Conditional entropy versus the words length for the Rössler system and its
reservoir computer clone. The parameters are the same as in the previous
figure. Open circles refers to the conditional entropy $h_{m}$ estimated from
a sequence generated by the genuine Rössler system. Filled black circles shows
$h_{m}$ estimated from sequences generated by the trained reservoir computer
with the same lengths of binary sequences.

This is no longer the case at greater $r$-values in the Lorenz model, as was
first shown computationally in Ref. [13], namely, stable and unstable
manifolds of saddle periodic orbits populating such quasi-attractors may no
longer cross transversely, thereby causing homoclinic tangencies that give
rise to the onset of stable periodic orbits (inside stability windows in the
parameter space) through saddle-node bifurcations.

Let us first discuss the change that the shape of the return map undergoes at
higher $r$-values, specifically the maps at $r=75$ are contrasted with
corresponding maps with $r=28$. One can observe first of all from Fig. 9B that
the map becomes extended with the bending section on the right. This fold is
where the map looses its expansion property due to the flat portion with zero
derivative. This is a direct cause for forthcoming saddle-node bifurcations
producing stable periodic orbits which can coexist with the chaotic pseudo-
hyperbolic subset (comprised of countably many periodic orbits) in bistability
or become global attractors in the phase space. The images of such [close to,
due to strong contraction in the transverse direction] 1D folded return map
can be also observed in the phase space projections of the chaotic set in the
Lorenz model Fig. 9B and its RC surrogate in Fig. 9C. These maps populated by
the black and red dots (given by $x^{\prime}_{max/min}=\ 0$) demonstrate the
distinctive hooks/folds where the wings of the attractors bend in the phase
space, unlike that case $r=28$ where they remain straight and flat.

The ground truth trajectory was integrated with Matlab’s ODE45 with fixed
timesteps at .005. The Lorenz ESN was based on Matlab code by M. Lukosevicius
[42]. Detailed parameters can be found on our github, cited below [43].

###  V.2 Rössler model’s replica

The qualitative results for the Rössler model’s fidelity test are displayed in
Fig. 10. Panels A and D correspond to the real model and panels C and E
correspond to the predicted trajectories. The $z_{max}$-return map fits
remarkably well, shown in panel B. The geometric distortion of the predicted
attractor is an artifact of smoothing which applied to ameliorate the
detection of "false" critical points when calculating the symbolic dynamics,
but the distortion was included for illustrative purposes.

Quantitative results can be seen in Fig 11. The conditional entropies show
tight agreement, even as block size increases, demonstrating the effectiveness
of the appraoch.

Figure 12: 3D projections of four surrogate attractors generated by the
Rössler random search, each displaying different qualitative properties. The
top legend shows the results of the short-run scoring system, and the
trajectory in phase space (background) illustrates each. The actual and
predicted conditional entropies are plotted, along with the $z_{map}$-return
maps superimposed on the right.

We note that because of slow-fast nature of the Rössler model and the lack of
symmetry, it was not easy to achieve such a good representation, compared to
with the Lorenz model. Several additional steps were necessary to obtain a
close replication of the Rössler model by ESN. In the particular example of
Fig.10, the true solution of the Rössler system was obtained for the
parameters $a=0.341$, $b=0.3$, $c=4.8$ and was integrated with an adaptive
3rd-order Bogacki-Shampine method [44, 45]. The solution was stored at time
steps of 0.02. Noise sampled from a normal distribution with mean 0 and
variance 0.3 was subsequently added. Selecting an adequate noise level proved
critical in obtaining consistent results, but setting it too high tended to
yield periodic orbits rather than chaotic attractors.

A random search over a broad space of hyper-parameters of the ESN was
conducted at a reservoir dimension of 80. Each ESN generated in the search was
trained over 30000 time steps after discarding a 2000 time step transient
using the Julia’s library ReservoirComputing.jl [46]. and the score was
computed over a test prediction of an additional 30000 time-steps according to
the following system:

  * •

Score = LZ$\times$Transient$\times$Entropy + Symmetry;

  * •

Lz is the absolute value of the difference between the LZ complexity of the
test and predicted sequences;

  * •

Transient is the squared error over the first 1000 steps, testing for initial
synchronization;

  * •

Entropy is the sum of the differences in block entropies over block sizes
ranging from 2 to 10;

  * •

The symmetry term is designed to disqualify symmetric systems which were
commonly generated. It simply adds $10^{9}$ if the trajectory dips below z =
-1.

Symbolic sequences were calculated as a sequence where local min are marked as
“0” and local maxima with $z>1$ are marked as “1”. Thus when the trajectory
completes a rotation around the origin without crossing the z=1 plane, a “00”
is recorded. Computing the local minima for the predicted trajectories was
complicated by the existence of “false” local minima, which were ameliorated
by three applications of a moving average filter with window 20. This seemed
to outperform other combinations of higher order Savitzky-Golay filters,
although no rigorous analysis was performed. Critical points were computed
using FindPeaks.jl with minimum prominence set at 0.1 with a minimum distance
of 150.

Each ESN generated for the search was trained over 30000 time steps after
discarding a 2000 time step transient using ReservoirComputing.jl and the
score was computed over a test prediction of an additional 30000 time-steps.
The top 50 scoring trajectories were inspected in phase space by eye alongside
a plot of their associated block entropies calculated over a longer trajectory
of 168,000 steps. The top candidate was chosen and the $z_{max}$ return map
and long run conditional block entropies were computed, with the final
entropies calculated over 9,968,000 time steps yielding sequences of lengths
69,994 and 66,322 characters for the true and simulated solutions
respectively.

##  VI Discussion

In conclusion we would like to reiterate that the symbolic approach used for
the characterization of complex chaotic dynamics, followed by the application
of the quantitative measure such as Shannon’s block entropies

The quick block-entropy method, or a longer Lempel-Ziv complexity approach or
even a simple algorithm based on a Markov transition matrix works
exceptionally well for the given purpose. One may argue that its only “flaw”
is the choice of the proper partition for the phase space to introduce the
symbolic description. While such a partition for the Lorenz model seems
evident when we study the homoclinic portion of its parameter space that
results in chaotic flip-flop patterns, nevertheless it becomes ineffective in
the period-doubling region of the parameter space for larger $r$-values. We
stress again that period-doubling bifurcations were out of the scope of our
study focused specifically on homoclinic ones in both system: the Lorenz and
Rössler models. Lastly, let us add that this approach should be viewed not as
a substitution but a novel and welcome addition to the computational
approached based on the Lyapunov characteristic exponents. Our findings are
meant to validate this assertion.

There are two pivotal impressions gained as a result of our experimentation
with ESNs. First, it must be acknowledged that while the echo state networks
presented here show excellent statistical agreement, that is not necessarily
an inherent property of the ESN framework. Curiously, the Lorenz model showed
immediate statistical agreement upon finding a parameter set. The Rössler
model, on the contrary, did not. Many Rössler ESN surrogates had a relatively
long initial synchronization, but performed poorly statistically. One
potential reason for that is that the Rössler model is a slow-fast systems,
where the $z$-variable happens to be the fast one with sharp low and high
peaks. Furthermore, the Lorenz model was much easier to train then the Rössler
one. Rössler surrogates often displayed inappropriate symmetry, collapsed into
periodic orbits of various lengths, or converged to visually similar but
qualitatively distinct systems. A few examples with accompanying statistical
portraits are displayed in the concluding Fig. 12 below.

We found ESN’s to have poor stability in general, and require a significant
amount of tweaking. We also found that the generally accepted guidelines for
hyper-parameter selection were not particularly useful. The spectral radius,
for example, did not have any obvious relation to performance. The best
surrogates were found widely distributed throughout the hyper-parameter search
space. That said, we make no claim to expertise in the domain, but do note it
requires patience and luck. Perhaps someone else with more skill can do
better. As such, a method for systematic reservoir construction with improved
stability and intuitive parameters would be a welcome contribution.

The second worthy note is that the fidelity test proved to be highly valuable
as a tool for automatically filtering ESN search results. Error alone is not a
sufficient test for chaotic systems due to sensitivity to initial conditions.
The incorporation of statistical quantification into the search criteria
proved to be immensely practical, albeit imperfect. The further development of
search techniques based on symbolic dynamics, perhaps combined with cross
validation or Monte-Carlo search, is a promising area for future research
which would yield immediate practical value for chaos prediction methods of
many varieties.

In any case, we can attest that reservoir computing passed all our fidelity
tests, a qualitative one based on the methods of dynamical systems and a
quantitative test based on statistical properties of chaotic solutions of the
original systems such as the Lorenz and Rössler models.

## Acknowledgements

A. Shilnikov acknowledges partial funding support from the Laboratory of
Dynamical Systems and Applications NRU HSE, of the Ministry of Science and
Higher Education of Russian Federation, grant No. 075-15-2019-1931.

## Data/Code Availability

All codes and data are freely available upon request.

## References

## References

  * Abarbanel  _et al._ [1993] H. D. Abarbanel, R. Brown, J. J. Sidorowich, and L. S. Tsimring, “The analysis of observed chaotic data in physical systems,” Reviews of modern physics 65, 1331 (1993).
  * Kantz and Schreiber [2004] H. Kantz and T. Schreiber,  _Nonlinear time series analysis_ , Vol. 7 (Cambridge university press, 2004).
  * Bradley and Kantz [2015] E. Bradley and H. Kantz, “Nonlinear time-series analysis revisited,” Chaos: An Interdisciplinary Journal of Nonlinear Science 25, 097610 (2015).
  * Anishchenko [1990] V. S. Anishchenko,  _Complex oscillations in simple systems_ (Nauka, Moscow, 1990).
  * Anishchenko  _et al._ [2003] V. Anishchenko, T. Vadivasova, G. Okrokvertskhov, and G. Strelkova, “Correlation analysis of dynamical chaos,” Physica A: Statistical Mechanics and its Applications 325, 199–212 (2003).
  * Anishchenko  _et al._ [2005] V. Anishchenko, T. Vadivasova, G. Okrokvertskhov, and G. Strelkova, “Statistical properties of dynamical chaos,” Physics-Uspekhi 48, 151 (2005).
  * Anishchenko  _et al._ [2001] V. S. Anishchenko, T. E. Vadivasova, A. S. Kopeikin, J. Kurths, and G. I. Strelkova, “Effect of noise on the relaxation to an invariant probability measure of nonhyperbolic chaotic attractors,” Physical review letters 87, 054101 (2001).
  * Anishchenko  _et al._ [2004] V. S. Anishchenko, T. E. Vadivasova, J. Kurths, G. A. Okrokvertskhov, and G. I. Strelkova, “Autocorrelation function and spectral linewidth of spiral chaos in a physical experiment,” Physical Review E 69, 036215 (2004).
  * Afraimovich, Bykov, and Shilnikov [1977] V. S. Afraimovich, V. V. Bykov, and L. P. Shilnikov, “The origin and structure of the Lorenz attractor,” Sov. Phys. Dokl. 22, 253–255 (1977).
  * Afraimovich and Shilnikov [1983] V. S. Afraimovich and L. P. Shilnikov, “Strange attractors and quasiattractors,” (1983).
  * Afraimovich, Bykov, and Shilnikov [1983] V. S. Afraimovich, V. V. Bykov, and L. P. Shilnikov, “On structurally unstable attracting limit sets of Lorenz attractor type,” Trans. Moscow Math. Soc. 44, 153–216 (1983).
  * Shilnikov [1986] L. P. Shilnikov, “Bifurcation theory and turbulence,” Methods of the Qualitative Theory of Differential Equations , 150–163 (1986).
  * Bykov and Shilnikov [1989] V. Bykov and A. Shilnikov, “On the boundaries of the domain of existence of the Lorenz attractor,” Methods of Qualitative Theory and Theory of Bifurcations. Gorky State University, Gorky , 151–159 (1989).
  * Gonchenko, Shil’nikov, and Turaev [1993] S. V. Gonchenko, L. P. Shil’nikov, and D. V. Turaev, “On models with non-rough poincaré homoclinic curves,” Physica D 1-4, 1–14 (1993).
  * Rössler [1976] O. E. Rössler, “An equation for continuous chaos,” Physics Letters A 57, 397–398 (1976).
  * Rössler [1979] O. E. Rössler, “Continuous chaos four prototype equations,” Annals of the New York Academy of Sciences 316, 376–392 (1979).
  * Shilnikov [1965] L. P. Shilnikov, “A case of the existence of a denumerable set of periodic motions,”  _Doklady Akademii Nauk_ , 160, 558–561 (1965).
  * Afraimovich  _et al._ [2014] V. S. Afraimovich, S. V. Gonchenko, L. M. Lerman, A. L. Shilnikov, and D. V. Turaev, “Scientific heritage of L.P. Shilnikov,” Regular and Chaotic Dynamics 19, 435–460 (2014).
  * Malykh  _et al._ [2020] S. Malykh, Y. Bakhanova, A. Kazakov, K. Pusuluri, and A. L. Shilnikov, “Homoclinic chaos in the rössler model,” Chaos: An Interdisciplinary Journal of Nonlinear Science 30, 113126 (2020).
  * Turaev and Shil’nikov [1998] D. V. Turaev and L. P. Shil’nikov, “An example of a wild strange attractor,” Sbornik: Mathematics 189, 291 (1998).
  * Turaev and Shilnikov [2008] D. Turaev and L. P. Shilnikov, “Pseudohyperbolicity and the problem on periodic perturbations of lorenz-type attractors,” in  _Doklady Mathematics_ , Vol. 77 (SP MAIK Nauka/Interperiodica, 2008) pp. 17–21.
  * Gonchenko, Shil’nikov, and Turaev [1997] S. Gonchenko, L. P. Shil’nikov, and D. Turaev, “Quasiattractors and homoclinic tangencies,” Computers & Mathematics with Applications 34, 195–227 (1997).
  * Gonchenko, Kazakov, and Turaev [2021] S. Gonchenko, A. Kazakov, and D. Turaev, “Wild pseudohyperbolic attractor in a four-dimensional lorenz system,” Nonlinearity 34, 2018 (2021).
  * Shannon [1948] C. E. Shannon, “A mathematical theory of communication,” The Bell system technical journal 27, 379–423 (1948).
  * Ebeling and Nicolis [1992] W. Ebeling and G. Nicolis, “Word frequency and entropy of symbolic sequences: a dynamical perspective,” Chaos, Solitons & Fractals 2, 635–650 (1992).
  * Ebeling [1997] W. Ebeling, “Prediction and entropy of nonlinear dynamical systems and symbolic sequences with lro,” Physica D: Nonlinear Phenomena 109, 42–52 (1997).
  * Eckmann and Ruelle [1985] J.-P. Eckmann and D. Ruelle, “Ergodic theory of chaos and strange attractors,” in  _The theory of chaotic attractors_ (Springer, 1985) pp. 273–312.
  * Sinai [2009] Y. Sinai, “Kolmogorov-Sinai entropy,” Scholarpedia 4, 2034 (2009), revision #91407.
  * Basios and Mac Kernan [2011] V. Basios and D. Mac Kernan, “Symbolic dynamics, coarse graining and the monitoring of complex systems,” International Journal of Bifurcation and Chaos 21, 3465–3475 (2011).
  * Herzel, Schmitt, and Ebeling [1994] H. Herzel, A. Schmitt, and W. Ebeling, “Finite sample effects in sequence analysis,” Chaos, Solitons & Fractals 4, 97–113 (1994).
  * Grassberger [1988] P. Grassberger, “Finite sample corrections to entropy and dimension estimates,” Physics Letters A 128, 369–373 (1988).
  * Ziv and Lempel [1978] J. Ziv and A. Lempel, “Compression of individual sequences via variable-rate coding,” IEEE transactions on Information Theory 24, 530–536 (1978).
  * Ruelle [1978] D. Ruelle, “An inequality for the entropy of differentiable maps,” Boletim da Sociedade Brasileira de Matemática-Bulletin/Brazilian Mathematical Society 9, 83–87 (1978).
  * Xing, Barrio, and Shilnikov [2014] T. Xing, R. Barrio, and A. Shilnikov, “Symbolic quest into homoclinic chaos,” International Journal of Bifurcation and Chaos 24, 1440004 (2014).
  * Pusuluri and Shilnikov [2018] K. Pusuluri and A. Shilnikov, “Homoclinic chaos and its organization in a nonlinear optics model,” Physical Review E 98, 040202 (2018).
  * Pusuluri, Meijer, and Shilnikov [2020] K. Pusuluri, H. G. E. Meijer, and A. L. Shilnikov, “Homoclinic puzzles and chaos in a nonlinear laser model,” J. Communications in Nonlinear Science and Numerical Simulations (2020).
  * Xing, Pusuluri, and Shilnikov [2021] T. Xing, K. Pusuluri, and A. L. Shilnikov, “Ordered intricacy of homoclinics of the L. Shilnikov saddle-focus in symmetric systems,” J. Chaos (2021).
  * Jaeger [2001] H. Jaeger, “The" echo state" approach to analysing and training recurrent neural networks-with an erratum note’,” Bonn, Germany: German National Research Center for Information Technology GMD Technical Report 148 (2001).
  * Schrauwen, Verstraeten, and Campenhout [2007] B. Schrauwen, D. Verstraeten, and J. Campenhout, “An overview of reservoir computing: Theory, applications and implementations,” (2007) pp. 471–482.
  * Pathak  _et al._ [2018] J. Pathak, B. Hunt, M. Girvan, Z. Lu, and E. Ott, “Model-free prediction of large spatiotemporally chaotic systems from data: A reservoir computing approach,” Physical review letters 120, 024102 (2018).
  * Pathak  _et al._ [2017] J. Pathak, Z. Lu, B. R. Hunt, M. Girvan, and E. Ott, “Using machine learning to replicate chaotic attractors and calculate lyapunov exponents from data,” Chaos: An Interdisciplinary Journal of Nonlinear Science 27, 121102 (2017).
  * Lukoševičius [2012] M. Lukoševičius, “A practical guide to applying echo state networks,” in  _Neural networks: Tricks of the trade_ (Springer, 2012) pp. 659–686.
  * [43] J. Scully, “https://github.com/jamesjscully/rossleresn,” .
  * Rackauckas and Nie [2017] C. Rackauckas and Q. Nie, “Differentialequations.jl–a performant and feature-rich ecosystem for solving differential equations in julia,” Journal of Open Research Software 5 (2017).
  * Bezanson  _et al._ [2017] J. Bezanson, A. Edelman, S. Karpinski, and V. B. Shah, “Julia: A fresh approach to numerical computing,” SIAM Review 59, 65–98 (2017).
  * [46] F. Martinuzzi, “Reservoircomputing.jl, ” https://github.com/SciML/ReservoirComputing.jl .