ordpy: A Python package for data analysis with permutation entropy and ordinal network methods
oordpy: A Python package for data analysis with permutation entropy and ordinalnetwork methods
Arthur A. B. Pessa ∗ and Haroldo V. Ribeiro † Departamento de F´ısica, Universidade Estadual de Maring´a – Maring´a, PR 87020-900, Brazil (Dated: February 17, 2021)Since Bandt and Pompe’s seminal work, permutation entropy has been used in several applicationsand is now an essential tool for time series analysis. Beyond becoming a popular and successfultechnique, permutation entropy inspired a framework for mapping time series into symbolic se-quences that triggered the development of many other tools, including an approach for creatingnetworks from time series known as ordinal networks. Despite the increasing popularity, the com-putational development of these methods is fragmented, and there were still no efforts focusingon creating a unified software package. Here we present ordpy , a simple and open-source Pythonmodule that implements permutation entropy and several of the principal methods related to Bandtand Pompe’s framework to analyze time series and two-dimensional data. In particular, ordpy implements permutation entropy, complexity-entropy plane, complexity-entropy curves, missing or-dinal patterns, ordinal networks, and missing ordinal transitions for one-dimensional (time series)and two-dimensional (images) data as well as their multiscale generalizations. We review sometheoretical aspects of these tools and illustrate the use of ordpy by replicating several literatureresults.
I. INTRODUCTION
Stemming from a combination of ideas from nonlin-ear time series analysis [1, 2] and information theory [3],permutation entropy was first introduced in 2002 byBandt and Pompe [4] as a simple, robust, and compu-tationally efficient complexity measure for time series.This complexity measure is defined as the Shannon en-tropy of a probability distribution associated with ordi-nal patterns estimated from partitions of a time series– a procedure known as the Bandt-Pompe symboliza-tion approach. Permutation entropy and its underlyingsymbolization approach have become increasingly pop-ular among researchers working with time series analy-sis, leading to successful applications in fields as diverseas biomedical sciences [5], econophysics [6], physical sci-ences [7], and engineering [8]. The uses of permutationentropy also span a large variety of goals such as mon-itoring the dynamical regime of a system [8], detectinganomalies in time series [7], characterizing time seriesdata [5], testing for serial independence [9], and are fur-ther documented in review articles by Zanin et al. [10],Reidl et al. [11], Amig´o et al. [12], and Keller et al. [13].Permutation entropy’s success is not limited to itspractical usage as this approach has inspired numeroustime series analysis tools. Some of these related meth-ods consider different quantifiers for the ordinal proba-bility distribution [14–23], generalize the Bandt-Pompesymbolization algorithm to evaluate ordinal structureson multiple temporal scales [24–27], include signal am-plitude information [28–31], and account for equal val-ues in time series [32, 33]. Other works have generalized ∗ arthur [email protected] † hvr@dfi.uem.br permutation entropy and its ordinal approach to two-dimensional data such as images [34, 35]. The ordinalpatterns underlying permutation entropy have also beenused for mapping time series and images into networksknown as ordinal networks [36–42].The original version of permutation entropy and itsvarious generalizations represent an essential and appeal-ing framework for data analysis, especially when consid-ering the increasing availability of large data sets [43]and the steady demand for reliable and computationallyefficient methods for extracting meaningful informationfrom these data sets [44, 45]. However, most methodsemerging from Bandt and Pompe’s seminal work lackfreely available computational implementation, and theexceptions are limited to a single or very few approaches.This poses significant limitations to further spreading theuse of these methods, especially to researchers in fieldswith less tradition in developing computational tools.Here we help to fill this gap by presenting ordpy – asimple and open-source Python module implementingseveral of the principal methods related to Bandt andPompe’s framework. This module has been designed tobe easily set up and installed as its only dependency is numpy [46], a fundamental Python library implementingarray objects and fast math functions that operate onthese objects. Beyond our preferences, the Python pro-gramming language has been chosen for its widespreaduse in scientific computing [46] and extensive communitysupport [47].We present ordpy ’s functions and illustrate their us-age along with a review of the pertinent theoretical de-velopments of permutation entropy and its related ordi-nal methods. Our work alternates between the mathe-matical description of the different techniques and thepresentation of functions and code snippets that imple-ment these data analysis tools. We further use ordpy to replicate several literature results. The source-code of a r X i v : . [ phy s i c s . d a t a - a n ] F e b ordpy is freely available on its git repository ( github.com/arthurpessa/ordpy ) together with the documenta-tion of all ordpy ’s functions ( arthurpessa.github.io/ordpy ). We can install ordpy using the Python PackageIndex (PyPI) via: $ pip install ordpy We further provide the code and data for replicatingall analyses presented in this article as a Jupyter note-book [48, 49] on ordpy ’s website.
II. AN OVERVIEW OF ORDINALDISTRIBUTIONS, PERMUTATION ENTROPY,AND COMPLEXITY-ENTROPY PLANE
We start by presenting a short review of Bandt andPompe’s seminal permutation entropy [4]. As we havealready mentioned, permutation entropy is the Shannonentropy of a probability distribution related to ordinal(or permutation) patterns estimated using sliding parti-tions over a time series. This probability distribution isthe so-called ordinal distribution or distribution of ordi-nal patterns, and the symbolization process used to esti-mate this distribution is the Bandt-Pompe approach. Todescribe this process, let us consider an arbitrary timeseries { x t } t = ,...,N x . First, we divide this time series into n x = N x − ( d x − ) τ x overlapping partitions comprised of d x > τ x ≥ d x and τ x , each data partition can berepresented by w p = ( x p , x p + τ x , x p + τ x , . . . , x p +( d x − ) τ x , x p +( d x − ) τ x ) , (1)where p = , . . . , n x is the partition index. The parame-ters d x and τ x are the only two parameters of the Bandt-Pompe method: d x is the embedding dimension [4] and τ x is the embedding delay [25]. It is worth remarking thatBandt and Pompe’s original proposal was restricted to τ x = et al. [50] and Zunino et al. [25].As we shall see, the choices of d x and τ x are important,and there is research exclusively devoted to determiningoptimal values for these parameters [11, 51].Next, for each partition w p , we evaluate the per-mutation π p = ( r , r , . . . , r d x − ) of the index numbers ( , , . . . , d x − ) that sorts the elements of w p in ascend-ing order, that is, the permutation of the index numbersdefined by the inequality x p + r ≤ x p + r ≤ ⋅ ⋅ ⋅ ≤ x p + r dx − .In case of equal values, we maintain the occurrence or-der of the partition elements, that is, if x p + r k − = x p + r k then r k − < r k for k = , . . . , d x − x t = ( , , , , , ) and set d x = τ x =
1. The firstpartition is w = ( , , , ) , and sorting its elements wefind 2 ≤ < < x + ≤ x + < x + < x + . Thus, thepermutation symbol (ordinal pattern) associated with w is π = ( , , , ) .After estimating the permutation symbols associatedwith all data partitions, we obtain a symbolic sequence { π p } p = ,...,n x . The ordpy ’s function ordinal sequence estimates this quantity as illustrated in the followingcode: >>> from ordpy import ordinal_sequence>>> x = [5, 3, 2, 2, 7, 9]>>> ordinal_sequence(x, dx=4, taux=1)array([[2, 3, 1, 0],[1, 2, 0, 3],[0, 1, 2, 3]]) The ordinal probability distribution P = { ρ i ( Π i )} i = ,...,n π is simply the relative frequency of all possible permuta-tions within the symbolic sequence, that is, ρ i ( Π i ) = number of partitions of type Π i in { π p } n x , (2)where Π i represents each one of the n π = d x ! differ-ent ordinal patterns. The following code shows howto obtain the ordinal distribution with ordpy ’s function ordinal distribution : >>> from ordpy import ordinal_distribution>>> x = [5, 3, 2, 2, 7, 9]>>> pis, rho = ordinal_distribution(x, dx=3)>>> pisarray([[0, 1, 2],[1, 2, 0],[2, 1, 0]])>>> rhoarray([0.5 , 0.25, 0.25]) The two arrays returned by ordinal distribution are the ordinal patterns and their correspond-ing relative frequencies, respectively. By default, ordinal distribution does not return non-occurringpermutations (that is, those with Π i = return missing modifies this behavior as in: >>> from ordpy import ordinal_distribution>>> x = [5, 3, 2, 2, 7, 9]>>> pis, rho = ordinal_distribution(x, dx=3,... return_missing=True)>>> pisarray([[0, 1, 2],[1, 2, 0],[2, 1, 0],[0, 2, 1],[1, 0, 2],[2, 0, 1]])>>> rhoarray([0.5 , 0.25, 0.25, 0. , 0. , 0. ]) Missing permutation symbols are always the latest ele-ments of the returned array.Having the ordinal probability distribution P , we cancalculate its Shannon entropy [3] and define the permu-tation entropy as S ( P ) = − n π ∑ i = ρ i ( Π i ) log ρ i ( Π i ) , (3)where log ( . . . ) stands for the base-2 logarithm. Permu-tation entropy quantifies the randomness in the order-ing dynamics of a time series such that S ≈ log n π in-dicates a random behavior, while S ≈ S is S max = log n π , we can further define the normalized per-mutation entropy as H ( P ) = S ( P ) log n π , (4)where the values of H are restricted to the interval [ , ] .The ordpy ’s function permutation entropy estimatesthe values of S and H directly from a time series as il-lustrated in: >>> from ordpy import permutation_entropy>>> x = [5, 3, 2, 2, 7, 9]>>> permutation_entropy(x)0.5802792108518123>>> permutation_entropy(x, normalized=False)1.0397207708399179 The permutation entropy function uses the base-2 log-arithm function by default; however, the parameter base can modify this behavior.The embedding dimension d x defines the number ofpossible permutations ( n π = d x ! ) , and following Bandtand Pompe’s recommendation [4], it is common to choosethe values of d x to satisfy the condition d x ! ≪ N x to havea reliable estimate of the ordinal probability distribution.Another less common choice is to use a value of d x suchthat 5 d x ! ≤ N x [53]. More recently, however, Cuesta-Frau et al. [51] have shown that these requirements on d x canbe considerably loosened in several situations related toclassification tasks.The permutation entropy framework was extended totwo-dimensional data by Ribeiro et al. [34] and Zuninoand Ribeiro [35]. To present this generalization, letus consider an arbitrary two-dimensional data array { y ut } u = ,...,N y t = ,...,N x whose elements may represent pixels of animage. We further define the embedding dimensions d x and d y along the horizontal and vertical directions (re-spectively), and the corresponding embedding delays τ x and τ y . Similarly to the one-dimensional case, we slicethe data array in partitions of size d x × d y defined by w qp = ⎛⎜⎜⎜⎜⎜⎝ y qp y q + τ y p . . . y q +( d y − ) τ y p y qp + τ x y q + τ y p + τ x . . . y q +( d y − ) τ y p + τ x ⋮ ⋮ ⋱ ⋮ y qp +( d x − ) τ x y q + τ y p +( d x − ) τ x . . . y q +( d y − ) τ y p +( d x − ) τ x ⎞⎟⎟⎟⎟⎟⎠ , (5) where the indexes p = , . . . , n x and q = , . . . , n y , with n x = N x −( d x − ) τ x and n y = N y −( d y − ) τ y , cover all n x n y data partitions. To associate a permutation symbol witheach two-dimensional partition, we flatten the partitions w qp line by line, that is, w qp = ( y qp , y q + τ y p , . . . , y q +( d y − ) τ y p ,y qp + τ x , y q + τ y p + τ x , . . . , y q +( d y − ) τ y p + τ x , . . . ,y qp +( d x − ) τ x , y q + τ y p +( d x − ) τ x , . . . , y q +( d y − ) τ y p +( d x − ) τ x ) . (6)As this procedure does not depend on a particular par-tition, we can simplify the notation by representing w qp as w qp = ( ˜ y , ˜ y , . . . , ˜ y d x d y − , ˜ y d x d y − ) , (7)where ˜ y = y qp , ˜ y = y q + τ y p , and so on. Then, we evaluatethe permutation symbol associated with each data par-tition as in the one-dimensional case to define the sym-bolic array { π qp } q = ,...,n y p = ,...,n x related to the data set. Fromthis array, we estimate the relative frequency for all n π = ( d x d y ) ! possible permutations Π i via ρ i ( Π i ) = number of partitions of type Π i in { π qp } n x n y , (8)where i = , . . . , n π , and so the ordinal probability distri-bution is P = { ρ i ( Π i )} i = ,...,n π . It is worth noticing thatthe ordering procedure defining the permutation symbolsis no longer unique as in the one-dimensional case. Forinstance, we would find a different symbolic array by flat-tening the partitions w qp column by column. However,different ordering procedures do not modify the set ofelements comprising the ordinal probability distribution(only their order is changed) [34].As in the one-dimensional case, the two-dimensionalpermutation entropy is simply the Shannon entropy ofthe ordinal distribution P = { ρ i ( Π i )} i = ,...,n π , so we canestimate the two-dimensional permutation entropy andits normalized version using Eqs. 3 and 4, respectively.Only the total number of possible ordinal patterns ( n π =( d x d y ) ! in the two-dimensional case) is modified.Similarly to the one-dimensional case, the values of d x and d y are usually constrained by the condition ( d x d y ) ! ≪ N x N y in order to obtain a reliable esti-mate of the ordinal distribution P [34, 35]. Natu-rally, this two-dimensional formulation recovers the one-dimensional case ( N y = d y = τ y =
1. In ordpy , the functions ordinal sequence , ordinal distribution and permutation entropy au-tomatically implement this two-dimensional generaliza-tion when the input data is a two-dimensional array asin: >>> from ordpy import ordinal_sequence,... ordinal_distribution, permutation_entropy>>> y = [[5, 3, 2], [2, 7, 9]]>>> ordinal_sequence(y, dx=2, dy=2) array([[[2, 1, 0, 3],[1, 0, 2, 3]]])>>> ordinal_distribution(y, dx=2, dy=2)(array([[1, 0, 2, 3],[2, 1, 0, 3]]), array([0.5, 0.5]))>>> permutation_entropy(y, dx=2, dy=2)0.21810429198553155 In addition to permutation entropy, the complexity-entropy plane proposed by Rosso et al. [14] is an-other popular time series analysis tool directly relatedto Bandt and Pompe’s symbolization approach. Thismethod was initially introduced for distinguishing be-tween chaotic and stochastic time series but has beensuccessfully used as an effective discriminating tool inseveral other contexts [54–58]. The complexity-entropyplane combines the normalized permutation entropy H (Eq. 4) with an intensive statistical complexity measure C (also estimated from the ordinal distribution) to builda two-dimensional representation space with the valuesof C versus H . The statistical complexity C used byRosso et al. is inspired by the work of Lopez-Ruiz etal. [59] and is defined by the product of the normal-ized permutation and a normalized version of the Jensen-Shannon divergence [60] between the ordinal distribu-tion P = { ρ i ( Π i )} i = ,...,n π and the uniform distribution U = { / n π } i = ,...,n π (it is worth remembering that n π isthe number of possible ordinal patterns). Mathemati-cally, we can write this measure as C ( P ) = D ( P, U ) H ( P ) D max , (9)where D ( P, U ) = S [( P + U )/ ] − S ( P ) − S ( U ) (10)is the Jensen-Shannon divergence and D max = − ( n π ! + n π ! log ( n π ! + ) − ( n π ! ) + log n π ! ) is a normalization constant. This latter constant ex-presses the maximum possible value of D ( P, U ) occurringfor P = { δ ,i } i = ,...,n π [61, 62].Differently from permutation entropy, the statisticalcomplexity C is zero in both extremes of order (whenonly one permutation symbol occurs) and disorder (whenall permutations are equally likely to happen). The valueof C quantifies structural complexity and provides addi-tional information that is not carried by the value of H .Furthermore, C is a nontrivial function of H in the sensethat for a given value of H , there exists a range of possi-ble values for C [14, 61, 62]. This happens because H and D are expressed by different sums of ρ i ( Π i ) and there isthus no reason for assuming a univocal relationship be-tween H and C .To better illustrate this feature, let us assume (for sim-plicity) we replace the Jensen-Shannon divergence by theEuclidean distance between P and U (as in the semi-nal work of Lopez-Ruiz et al. [59]), that is, D ( P, U ) = ∑ n π i = ( ρ i ( Π i ) − / n π ) . In this case, the statistical com-plexity is C ( P ) ∝ − ( n π ∑ i = ρ i ( Π i ) log ρ i ( Π i )) ( n π ∑ i = ( ρ i ( Π i ) − / n π ) ) , and we can readily observe that different ordinal dis-tributions P = { ρ i ( Π i )} i = ,...,n π may lead to the samevalue of H but different values of C (or vice-versa).Let us further consider a particular ordinal distributionwith three possible permutation symbols (this would beequivalent to have d x ! = ( d x d y ) ! =
3, if possi-ble), that is, P = { a, b, − ( a + b )} , where a > b > ( a + b ) ≤ P ). For this case, we have S = − a log a − b log b − [ − ( a + b )] log [ − ( a + b )] and D = ( a − / ) + ( b − / ) + ([ − ( a + b )] − / ) . Thus, forinstance, if a = .
79 and b = .
18 or a = .
80 and b = . H = S / log 3 ≈ .
55, but differ-ent values for D (0 .
32 in the first case and 0 .
33 in thesecond) and, consequently, for C .In ordpy , the complexity entropy function simulta-neously estimates the values of H and C from time seriesas illustrated in: >>> from ordpy import complexity_entropy>>> complexity_entropy([4,7,9,10,6,11,3], dx=2)(0.9182958340544894, 0.06112816548804511) Furthermore, the complexity-entropy plane was general-ized for two-dimensional data [34, 35] (notice that theonly changes are related to the process of estimatingthe ordinal distribution) and the complexity entropy function also accepts two-dimensional arrays as input asshown in: >>> from ordpy import complexity_entropy>>> complexity_entropy([[1,2,1],[8,3,4],[6,7,5]],... dx=2, dy=2)(0.3271564379782973, 0.2701200547320647)
III. APPLICATIONS OF BANDT ANDPOMPE’S FRAMEWORK WITH
ORDPY
This section presents more engaging applications of ordpy ’s functions by replicating literature results. Westart by estimating the ordinal probability distributionsof two different time series of stochastic and chaotic na-ture, namely, a random walk with Gaussian steps and thelogistic map at fully developed chaos (see Appendix Afor definitions). We choose these two examples be-cause their ordinal distributions are exactly known forsome combinations of the embedding parameters [63, 64].More specifically, for d x = τ x =
1, the probabil-ity distributions associated with the permutation sym-bols {( , , ) , ( , , ) , ( , , ) , ( , , ) , ( , , ) , ( , , )} are P walk = { / , / , / , / , / , / } and P logistic ={ / , / , / , / , / , } for the random walk [64]and the logistic map [63], respectively.
012 021 102 120 201 2100 P r obab ili t y , i ( i ) (a) random walkSimulationTheory
012 021 102 120 201 210 P r obab ili t y , i ( i ) (b) logistic mapSimulationTheory FIG. 1. Probability distributions of ordinal patterns for stochastic and deterministic series. (a) Comparison between the empir-ical probability distribution of ordinal patterns obtained from a simulated Gaussian random walk with 10 steps and the exactdistribution P walk (dashed horizontal lines) for d x = τ x =
1. (b) Comparison between the empirical probability distributionof ordinal patterns obtained from 10 iterations of the logistic map at fully developed chaos and the exact distribution P logistic (dashed horizontal lines) for d x = τ x =
1. All results in this figure can be replicated by running a Jupyter notebookavailable at ordpy ’s webpage.
To numerically estimate these two ordinal distribu-tions, we generate a time series from a Gaussian randomwalk process and another time series from iterations ofthe fully chaotic logistic map. For both cases, we havesimulated one realization of each process with 10 obser-vations and used the ordinal distribution function.Figure 1 shows that the exact ordinal distributions arein excellent agreement with simulated results obtainedwith ordpy . It is intriguing to observe that the ordinalpattern ( , , ) (“descending permutation”) does not oc-cur in the logistic series (it has probability zero). Thisfact is best understood as a feature directly associatedwith the intrinsic determinism of the logistic map dy-namics [63, 65]. As we shall discuss in the next section,investigations about such “missing ordinal patterns” arealso useful for characterizing time series dynamics.To better illustrate the use of the permutation entropy function, we partially repro-duce Bandt and Pompe’s analysis of the logistic map(Fig. 2 of Ref. [4]). We generate time series consisting of10 iterations of the logistic map for each value of param-eter r ∈ { . , . , . , . . . , . } (see Appendix Afor definitions). Next, we calculate the permutationentropy S for each of these 5001 time series using permutation entropy with embedding parameters d x = τ x =
1. We further divide the permutationentropy by 5 to obtain the permutation entropy persymbol of order 6, that is, h = S / h as a function of the parameter r . We note that the permutation entropy per symbolhas an overall increasing trend with the parameter r marked by abrupt drops in intervals of r related toperiodic behaviors. As noticed by Bandt and Pompe,the behavior of the permutation entropy is similar to the one observed for the Lyapunov exponent [4].In another example with permutation entropy , wereplicate a numerical experiment of Cao et al. [50] (seetheir Fig. 1) that searches for dynamical changes in thetransient logistic map time series (see Appendix A fordefinitions). This problem illustrates the role of the em-bedding delay τ x . As in the original article, we iterate thetransient logistic map starting with the initial condition x = .
65 and incrementing the logistic parameter r from2 . − . This process generates atime series with 120001 observations as shown in Fig. 2c.Using this time series, we calculate the normalized per-mutation entropy within a sliding window with 1024 ob-servations for the embedding dimension d x = τ x = τ x = et al. [50], we denote the permutation en-tropy values by H [ r ( t )] , where r ( t ) represents the lo-gistic parameter at the end of the sliding window. Fig-ure 2d shows the values of H [ r ( t )] , where abrupt changesclearly associate with dynamical changes observed in thetime series (Fig. 2c). Despite the overall similarities, wenote that the embedding delay τ x = τ x = r ≈ .
56) is missed when τ x = τ x = et al. [34]. To illustrate the use of the permutation entropy function with two-dimensionaldata, we replicate a numerical experiment related to Isingsurfaces (see Appendix B for definitions) present in thatwork (Fig. 8 of Ref. [34]). These surfaces represent theaccumulated sum of spin variables of the canonical two-dimensional Ising model in a Monte Carlo simulation.Figure 2e shows four examples of these surfaces (square Parameter, r P opu l a t i on , x (a) Parameter, r ( t ) P opu l a t i on , x [ r ( t )] (c) (e) T r = 0.8 T r = 0.9 T r = 1.0 T r = 1.1 Parameter, r P e r m u t a t i on en t r op y pe r sy m bo l , h (b) d x = 6 Parameter, r ( t ) P e r m u t a t i on en t r op y , H [ r ( t )] r = 3.56 (d) d x = 5, x = 1 d x = 5, x = 2 Reduced temperature, T r P e r m u t a t i on en t r op y , H (f) d x = 3 d y = 2 d x = 2 d y = 3 FIG. 2. Permutation entropy of one- and two-dimensional data. (a) Bifurcation diagram of the logistic map for r between 3 . − . (b) Permutation entropy per symbol of order six ( h ) calculated from logistic time series with 10 observations (random initial conditions) and r ∈ { . , . , . , . . . , . } . The embedding parameters are d x = τ x = x = .
65, and by incrementing the logisticparameter r at each iteration from 2 . − . Despite appearing very similar to a bifurcation diagram, thisresult refers to a time series where each observation x [ r ( t )] corresponds to a value r ( t ) . (d) Dependence of the normalizedpermutation entropy estimated within a sliding window with 1024 observations of the original time series. Here r ( t ) representsthe logistic parameter at the end of each sliding window. The different curves show the results for d x = τ x = d x = τ x = r = .
56 indicates the period-8 to period-16 bifurcation. (e) Ising surfacesobtained after 10 Monte Carlo steps with reduced temperatures T r ∈ { . , . , . , . } . In these surfaces, dark gray shadesindicate high lattice sites while light gray regions indicate the opposite. (f) Normalized permutation entropy as a function ofthe reduced temperature T r ∈ { . , . , . . . , . } for Ising surfaces of size 250 ×
250 obtained after 10 Monte Carlo steps. Thedifferent curves show the results for embedding parameters d x = d y = d x = d y = τ x = τ y =
1. All results in this figure can be reproduced by running a Jupyter notebook available at ordpy ’s webpage. lattices of size 250 × Monte Carlosteps for different reduced temperatures T r . We noticenon-trivial patterns emerging when the reduced tempera-ture is equal to the critical temperature ( T r = ) of phasetransition for the Ising model [66]. Following the originalarticle, we generate Ising surfaces (size 250 × T r ∈ { . , . , . . . , . } and calculatetheir normalized permutation entropy with d x = d y =
2, and d x = d y =
3, both for τ x = τ y =
1. Inagreement with Ribeiro et al. [34], Fig. 2f shows that thepermutation entropy precisely identifies the phase transi-tion of the Ising model (the sudden decrease around thecritical temperature) and that these Ising surfaces aresymmetric under reversal of the embedding dimensions.The complexity entropy function simultaneously cal-culates the permutation entropy and the statistical com-plexity from time series and image data. To illustrateits usage, we partially reproduce the results of Rosso etal. [14] (Fig. 1 in that work) on distinguishing chaoticfrom stochastic time series. By following their article, we iterate four discrete maps to generate chaotic series.Specifically, we obtain chaotic time series from skew tentmap (parameter w = . x -component,parameters a = . b = . ( r = ) ,and Schuster map (parameter z ∈ { / , , / } ) – see Ap-pendix A for definitions. We further generate stochasticseries from three stochastic processes: noises with 1 / f − k power spectrum (for k ∈ { . , . , . . . , . } ), fractionalBrownian motion (Hurst exponent h ∈ { . , . , . . . , . } ),and fractional Gaussian noise (also h ∈ { . , . , . . . , . } )– see Appendix B for definitions. For each of these mapsand stochastic processes, we generate ten time series with2 observations and random initial conditions. Next,we use complexity entropy with embedding parameters d x = τ x = et al. [14], Fig. 3ashows that chaotic series usually have high complex-ity and low entropy values. Stochastic time series, in Permutation entropy, H S t a t i s t i c a l c o m p l e x i t y , C z =3/2 z =2 z =5/2 k =0 k =1 k =2 k =3 h =0.1 h =0.5 h =0.9 d x = 6 (a) Skew Tent MapHenon MapLogistic MapSchuster Map K-NoisefBmfGn
Permutation entropy, H S t a t i s t i c a l c o m p l e x i t y , C PollockReinhardt Amaral d x = 2 d y = 2 (b) FIG. 3. Complexity-entropy plane for one- and two-dimensional data. (a) Average values of the statistical complexity C versus the normalized permutation entropy H (over ten realizations) estimated from time series of chaotic maps and stochasticprocesses. The embedding parameters are d x = τ x =
1. The solid lines represent the maximum and minimal possiblevalues of complexity for a given entropy (for d x = τ x = d x = d y = τ x = τ y =
1. All data and code necessary to reproduce this figure are availablein a Jupyter notebook at ordpy ’s webpage. turn, display high entropy and intermediary complex-ity values. It is also interesting to note that stochas-tic time series approach the lower-right corner of thecomplexity-entropy plane ( H → C →
0) as theserial auto-correlation decreases [14]. These results alsoillustrate that some stochastic and chaotic series havevery similar entropy values but different statistical com-plexity (for instance, fractional Brownian motion with h = . z = / ordpy , the functions maximum complexity entropy and minimum complexity entropy generate these curves, asshown in the following code snippet: >>> from ordpy import maximum_complexity_entropy,... minimum_complexity_entropy>>> maximum_complexity_entropy(dx=4)array([[-0. , -0. ],[ 0.21810429, 0.19670592],[ 0.34568712, 0.28362016],...[ 0.98660828, 0.02388382]])>>> minimum_complexity_entropy(dx=4)array([[-0.00000000e+00, -0.00000000e+00],[ 2.67076969e-02, 2.55212327e-02],...[ 1.00000000e+00, -3.66606083e-16]]) The complexity entropy function also works withtwo-dimensional data, and to illustrate its usage, we fol-low Sigaki et al. [67] and use the complexity-entropyplane to investigate patterns in art paintings. Due tothe large-scale of the data analyzed by Sigaki et al. andto keep our examples self-contained, we do not reproducetheir original results but simply use their ideas to illus-trate how complexity and entropy extract useful informa- tion from images. To do so, we handpick three paintingsfrom wikiart.org (in the original article, the authors stud-ied 137,364 images obtained from the same webpage).These are a Color Field Painting artwork (Blue, 1953by Ad Reinhardt, image size 768 ×
435 [68]), a BrazilianModernist artwork (Abaporu, 1928 by Tarsila do Amaral,image size 1200 × × d x = d y = τ x = τ y = et al. [67], these re-sults show that paintings portraying objects with clearlydefined borders (such as the squares in Reinhardt’s art-work) tend to present large values of statistical complex-ity and low values of entropy. On the other extreme,paintings with smudged and diffuse contours (such asPollock’s drip paintings) have high entropy and low com-plexity values. Between these somewhat opposite behav-iors, we have a whole continuum of images, as exempli-fied here by the work of the Brazilian painter Tarsila doAmaral. As argued by Sigaki et al. [67], the complexity-entropy plane maps the local degree of order of artworksinto a scale of order-disorder and simplicity-complexitythat is similar to qualitative descriptions of artworks pro-posed by art historians such as W¨olfflin (the linear versuspainterly dichotomy) and Riegl (the haptic versus opticdichotomy).
60 2000 4000 6000
Time series length, N x N u m be r o f m i ss i ngpa tt e r n s , logistic maprandom walk (a) Noise amplitude, N u m be r o f m i ss i ngpa tt e r n s , noisy logistic map (b) d x = 5 d x = 6 d x = 5 d x = 6 d x = 5 d x = 6 FIG. 4. Missing ordinal patterns in time series. (a) Number of missing ordinal patterns ( η ) in random walk (blue) andlogistic map (red) time series as a function of sequence length ( N x ) for embedding parameters d x = d x =
6, bothwith τ x =
1. Results represent the average number of missing permutations over ten time series replicas for each N x ∈{ , , , . . . , } . (b) Dependence of the number of missing ordinal patterns on the noise intensity ( ξ ) for noisy logistictime series with 6000 observations. The noise added to the logistic series is uniformly distributed in the interval [− ξ, ξ ] with ξ ∈ { , . , . , . . . , . } . Results represent average values over ten time series replicas for each noise level. The embeddingdimensions are indicated within the plot and the embedding delay is τ x =
1. We use random initial conditions and set theparameter r = ordpy ’s webpage. IV. MISSING ORDINAL PATTERNS
As we have commented, the logistic map at fully de-veloped chaos does not exhibit the “descending permuta-tion” ( , , ) for d = etal. [63, 65] are seminal in this regard, and by followingtheir classification, we can divide these forbidden ordinalpatterns into two categories: true or false [65]. True for-bidden patterns (such as the ( , , ) in the logistic map)are a fingerprint of determinism in a time series dynam-ics and represent an intrinsic feature of the underlyingdynamical process [63]; that is, these patterns are not anartifact related to the finite length of empirical observa-tions. In turn, false forbidden patterns are related to thefinite length of time series [65] and can emerge even fromcompletely random processes.This distinction is not straightforward when dealingwith empirical data, but a typical analysis in this contextconsists in investigating the number of missing patterns( η ) as a function of the time series length ( N x ). Thebehavior of this curve is useful for discriminating timeseries. In ordpy , the missing patterns function identi-fies missing ordinal patterns and estimates their relativefrequency as in: >>> from ordpy import missing_patterns>>> missing_patterns([4,7,9,10,6,11,3,5,6,2,3,1],... dx=3)(array([[0, 2, 1], [2, 1, 0]]),0.3333333333333333) To better illustrate the use of this function, we inves-tigate missing ordinal patterns in time series obtainedfrom the logistic map at fully developed chaos ( r =
4) andGaussian random walks. In both cases, we use the em-bedding dimensions d x = d x = τ x =
1) andseries lengths N x ∈ { , , , . . . , } . Figure 4ashows the results. We observe that the number of miss-ing ordinal patterns approaches zero as the time serieslength of random walks increases. Conversely, the num-ber of missing permutations related to the logistic mapdisplays an initial decay with the time series length butit rapidly saturates in considerably large numbers, indi-cating that these missing patterns are intrinsically associ-ated with the underlying determinism of the process [65].In another application with the missing patterns function, we replicate a result of Amig´o et al. [65] (Fig. 4in their work) to further show that the number of miss-ing patterns is a good indicator of determinism in timeseries [53, 65]. By following the original work, we gener-ate time series from the logistic map at fully developedchaos (6000 iterations) and add to them uniformly dis-tributed noise in the interval [− ξ, ξ ] , where ξ is the noiseamplitude. Next, we estimate the average number ofmissing patterns (over ten time series replicas) for eachnoise level ξ ∈ { , . , . , . . . , . } , and embeddingdimensions d x = d x = τ x = ξ for both embedding dimensions.We observe that the number of missing patterns relatedto these deterministic series contaminated with noise ap-proaches zero as noise amplitude grows. However, sig-nificantly higher noise levels are necessary to remove allsigns of determinism expressed by the lack of permuta-tion patterns when d x = V. TSALLIS AND R´ENYI ENTROPY-BASEDQUANTIFIERS OF THE ORDINALDISTRIBUTION
In addition to Shannon’s entropy and the statisticalcomplexity, researchers have proposed to use other quan-tifiers of the ordinal probability distribution. As we haveexplicitly verified for the statistical complexity, these dif-ferent quantifiers are supposed to extract additional in-formation from a time series dynamics that is not cap-tured by permutation entropy and statistical complex-ity. In this context, a productive approach is to considerparametric generalizations of Shannon’s entropy, such asthose proposed by Tsallis [75] and R´enyi [76]. The workof Zunino et al. [15] was the first to consider the Tsal-lis entropy in place of Shannon’s entropy to define theTsallis permutation entropy as S β ( P ) = β − n π ∑ i = ( ρ i ( Π i ) − ρ i ( Π i ) β ) , (11)where β is a real parameter ( β → S max β = −( n π ) − β β − . Thus, the normalizedTsallis permutation entropy is H β ( P ) = ( β − ) S β ( P ) − ( n π ) − β . (12)Similarly, Liang et al. [19] have proposed the R´enyipermutation entropy S α ( P ) = − α ln ( n π ∑ i = ρ i ( Π i ) α ) , (13)where α > α → S max α = ln n π , as theusual Shannon entropy). Thus, the normalized R´enyipermutation entropy is H α ( P ) = S α ( P ) ln n π . (14)In both cases, the generalized entropic form is mono-parametric and has a term where the ordinal proba-bilities appear raised to the power of the entropic pa-rameter (that is, ρ i ( Π i ) β and ρ i ( Π i ) α ). These param-eters assign different weights to the underlying ordinalprobabilities, allowing us to access different dynamicalscales and produce a family of quantifiers for the ordi-nal distribution. In ordpy , the tsallis entropy and renyi entropy functions implement these two quanti-fiers as in: >>> from ordpy import tsallis_entropy,... renyi_entropy>>> tsallis_entropy([4,7,9,10,6,11,3], q=[1,2],... dx=2) In a similar direction, there are also the develop-ments of complexity-entropy curves proposed by Ribeiro et al. [22] and Jauregui et al. [23]. These works havefurther extended the complexity-entropy plane conceptby considering the Tsallis and R´enyi entropies combinedwith proper generalizations of statistical complexity [62].Thus, instead of having a single point in the complexity-entropy plane for a given time series, Ribeiro et al. [22]and Jauregui et al. [23] have created parametric curvesby varying the entropic parameter ( β or α ) and simulta-neously estimating the generalized entropy and the gen-eralized statistical complexity.To define the Tsallis complexity-entropy curves [22],we first extend the statistical complexity (Eq. 9) usingthe Tsallis entropy, that is, C β ( P ) = D β ( P, U ) H β ( P ) D max β , (15)where D β ( P, U ) = K β ( P ∣ P + U ) + K β ( U ∣ P + U ) (16)is the Jensen-Tsallis divergence [62] written in terms ofthe corresponding Kullback-Leibler divergence [62, 77] K β ( V ∣ R ) = β − n π ∑ i v βi [ r − βi − v − βi ] , (17)where V = { v i } i = ,...,n π and R = { r i } i = ,...,n π are two ar-bitrary distributions. In Eq. 15, D max β = − β n π − ( + n π ) − β − n π ( + / n π ) − β − n π + − β n π ( − β ) is a normalization constant representing the maxi-mum possible value of D β ( P, U ) that occurs for P ={ δ ,i } i = ,...,n π (as in the usual Jensen-Shannon diver-gence). By following Ribeiro et al. [22], we con-struct a parametric representation of the ordered pairs ( H β ( P ) , C β ( P )) for β >
0, obtaining the Tsalliscomplexity-entropy curves.Similarly, to define the R´enyi complexity-entropycurves [23], we generalize the statistical complexity inR´enyi’s formalism as C α ( P, U ) = D α ( P, U ) H α ( P ) D max α , (18)where D α ( P, U ) = K α ( P ∣ P + U ) + K α ( U ∣ P + U ) (19)0is the Jensen-R´enyi divergence [62] written in terms of K α ( V ∣ R ) = α − ( d ! ∑ i = v αi r − αi ) , (20)the corresponding Kullback-Leibler divergence forR´enyi’s entropy [62, 78]. The normalization constant D max α = ( α − ) ln [ ( n π + ) − α + n π − n π ( n π + n π ) − α ] corresponds to the maximum possible value of D α ( P, U ) occurring for P = { δ ,i } i = ,...,n π (as in the usualJensen-Shannon divergence). Again, we can con-struct a parametric representation of the ordered pairs ( H α ( P ) , C α ( P )) for α >
0, obtaining the R´enyicomplexity-entropy curves proposed by Jauregui etal. [23].In ordpy , the functions tsallis complexity entropy and renyi complexity entropy implement the Tsallisand R´enyi complexity-entropy curves as shown in the fol-lowing code snippet: >>> from ordpy import tsallis_complexity_entropy,... renyi_complexity_entropy>>> tsallis_complexity_entropy([4,7,9,10,6,11,3],... dx=2, q=[1,2])
To better illustrate the use of these ordpy ’s functions,we replicate some numerical experiments involving thelogistic map and random walks present in the originalworks of Ribeiro et al. [22] (Figs. 1 and 6 in that work)and Jauregui et al. [23] (Figs. 1 and 3 in that work).We start by generating time series from the logistic mapat fully developed chaos ( r =
4, random initial condition)and a Gaussian random walk. For the logistic map series,we discard the first 10 iterations to avoid transient ef-fects and iterate other 10 steps. The random walk seriesalso has 10 observations. By using these time series, wegenerate their corresponding Tsallis complexity-entropycurves for d x = τ x = log-spacedvalues of the entropic parameter β between 0 .
01 and 100for the logistic map, and between 0 .
001 and 100 for therandom walk.Figures 5a and 5b show the empirical complexity-entropy curves in comparison with their exactshape (dashed lines). These theoretical curvescan be estimated for these time series becausethe ordinal distributions of the logistic map( P logistic = { / , / , / , / , / , } ) and ran-dom walks ( P walk = { / , / , / , / , / , / } ) areexactly known for d x = et al. [22], random series tend to form closed complexity-entropy curves (Fig. 5a),while chaotic time series are usually represented by opencomplexity-entropy curves (Fig. 5b). These featuresemerge as a direct consequence of the existence or notof missing ordinal patterns captured by the limitingbehavior of H β as β → β → ∞ [22].By following a similar approach, we also estimate theR´enyi complexity-entropy curves for the two previoustime series for d x = τ x = et al. [23]have found that the initial curvature of R´enyi complexity-entropy curves ( dC α / dH α for small α ) can be used asan indicative of determinism in time series. Specifically,they found that positive curvatures are associated withtime series of stochastic nature, while negative ones arerelated to chaotic phenomena. This pattern also occursin the results of Figs. 5c and 5d. VI. ORDINAL NETWORKS
Among the more recent developments related to theBandt-Pompe framework, we have the so-called ordinalnetworks. First proposed by Small [36] for investigatingnonlinear dynamical systems, ordinal networks belong toa more general class of methods designed to map timeseries into networks, collectively known as time seriesnetworks [79]. Beyond counting ordinal patterns, thisapproach considers first-order transitions among ordinalsymbols within a symbolic sequence. In this network rep-resentation, the different ordinal patterns occurring in adata set are mapped into nodes of a complex network.The edges between nodes indicate that the associatedpermutations symbols are adjacent to each other in asymbolic sequence. Furthermore, edges can be directedaccording to the temporal succession of ordinal symbolsand weighted by the relative frequencies in which the cor-responding successions occur in a symbolic sequence [37].After applying the Bandt-Pompe method with embed-ding parameters d x and τ x to a time series { x t } t = ,...,N x and obtaining the symbolic sequence { π p } p = ,...,n x , wecan define the elements of the weighted adjacency ma-trix of the corresponding ordinal network as [37, 39] ρ i,j = total of transitions Π i → Π j in { π p } p = ,...,n x n x − , (21)where i, j = , , . . . , n π (with n π = d x !), Π i and Π j rep-resent all possible ordinal patterns, and the denominator n x − ordpy ,the ordinal network function returns the nodes, edges,and edge weights of an ordinal network mapped from atime series as in:1 Permutation entropy, H S t a t i s t i c a l c o m p l e x i t y , C logistic map d x = 3 (a) TheorySimulation
Permutation entropy, H S t a t i s t i c a l c o m p l e x i t y , C random walk d x = 3 (b) TheorySimulation
Permutation entropy, H S t a t i s t i c a l c o m p l e x i t y , C logistic map d x = 4 (c) Permutation entropy, H S t a t i s t i c a l c o m p l e x i t y , C random walk d x = 4 (d) FIG. 5. Tsallis and R´enyi complexity-entropy curves. Tsallis complexity-entropy curves for time series obtained from (a)the logistic map at fully developed chaos and (b) a Gaussian random walk, both with embedding parameters d x = τ x =
1. The solid lines represent the empirical results and the dashed lines indicate the exact form of these complexity-entropycurves. Panels (c) and (d) show the R´enyi complexity-entropy curves estimated from the same two time series with embeddingparameters d x = τ x =
1. In all panels, star markers indicate the beginning of the curves ( β ≈ α ≈ β and α ). Data and code necessary to reproduce these results areavailable in a Jupyter notebook at ordpy ’s webpage. >>> from ordpy import ordinal_network>>> ordinal_network([4,7,9,10,6,11,8,3,7], dx=2,... normalized=False)(array(['0|1', '1|0'], dtype=' 1) or vertically( π qp → π qp + for p = , . . . , n x − 1) adjacents in the sym-bolic sequence. The directed link between a pair of per-2mutation symbols (Π i and Π j ) is weighted by the totalnumber of occurrences of this particular transition in thesymbolic sequence. Thus, the weighted adjacency ma-trix representing the ordinal network mapped from two-dimensional data is [40] ρ i,j = total of transitions Π i → Π j in { π qp } q = ,...,n y p = ,...,n x n x n y − n x − n y , (22)where i, j = , . . . , n π (with n π = ( d x d y ) !) and the de-nominator represents the total number of horizontal andvertical transitions. The ordinal network function alsohandles two-dimensional data as in: >>> from ordpy import ordinal_network>>> ordinal_network([[1,2,1],[8,3,4],[6,7,5]],... dx=2, dy=2, normalized=False)(array(['0|1|3|2', '1|0|2|3', '1|2|3|0'],dtype=' An intriguing feature of ordinal networks is the ex-istence of intrinsic connectivity constraints [39, 40] in-herited from Bandt and Pompe’s symbolization method.These constraints are directly related to the fact that ad-jacent partitions share elements, such that ordering rela-tions in one partition are partially carried out to neigh-boring partitions. For one-dimensional data, these re-strictions imply that all nodes in an ordinal network havein-degree and out-degree limited to numbers between 0and d x ; consequently, the maximum number of edgesis d x × ( d x ! ) [39]. Ordinal networks mapped from one-dimensional data can only have self-loops in nodes asso-ciated with solely ascending or solely descending ordinalpatterns [39].The horizontal and vertical transitions related to net-works mapped from two-dimensional data impose similar but trickier connectivity constraints [40]. In this case, themaximum number of outgoing connections emerging fromhorizontal and vertical transitions are C ( d x d y , d y ) × d y !and C ( d x d y , d x )× d x !, respectively [40]. However, the setsof horizontal and vertical transitions are not disjoint, andtheir union defines all possible outgoing edges. Finding ageneral expression for the latter set operation is cumber-some because it depends on the ordinal pattern associ-ated with the node under analysis. Thus, while limited,the maximum number of edges varies among the ordi-nal patterns and needs to be numerically obtained [40].Furthermore, differently from the one-dimensional case,ordinal networks mapped from two-dimensional data candisplay self-loops in several nodes [40].A direct consequence of these intrinsic connectivityconstraints is that ordinal networks mapped from com-pletely random arrays (in one or two dimensions) are notrandom graphs [39, 40]. Even more counter-intuitive isthe existence of different edge weights in random ordi-nal networks, albeit all permutations are equiprobable inrandom arrays [39, 40]. This non-trivial property resultsfrom the fact that, among all possible amplitude relationsinvolved in an ordinal transition between a fixed per-mutation and all its possible neighboring permutations,some permutations appear more than once. For one-dimensional data, random ordinal networks only havetwo different edge weights: 1 /( d x + ) ! and 2 /( d x + ) ! (thedenominator represents the sum of weights) [39]. A ruleof thumb for determining the edges with double weightis to pick all transitions in which the symbol equal to“ d x − 1” in the next permutation fits the position of thesymbol “0” in the first permutation [39], for instance,the edge weight between permutations ( , , , ) and ( , , , ) has double weight. Ordinal networks mappedfrom two-dimensional random data have more than twodifferent edge weights, and there is no simple rule (atleast up to now) for obtaining these weights [40]. How-ever, these values can be numerically calculated by ex-plicitly considering each possible ordinal pattern [40].In ordpy , the random ordinal network function gen-erates the exact form of ordinal networks expectedfrom the mapping of one- and two-dimensional ran-dom data with arbitrary embedding dimensions ( d x and d y ). The following code illustrates the usage of random ordinal network : >>> from ordpy import random_ordinal_network>>> random_ordinal_network(dx=2)(array(['0|1', '1|0'], dtype=' Similarly, embedding delays larger than one modify howelements are shared among partitions and impose con-nectivity constraints to high-order transitions. The random ordinal network function is thus restricted tothe case τ x = τ y = s i = − ∑ j ∈O i ρ ′ i,j log ρ ′ i,j , (23)where the index i refers to a node related to a given per-mutation Π i , ρ ′ i,j = ρ i,j / ∑ k ∈O i ρ i,k represents the renor-malized probability of transitioning from node i to node j (permutations Π i and Π j ), and O i is the outgoing neigh-borhood of node i (set of all edges leaving node i ). Thisquantity measures the determinism of ordinal transitionsat the node level such that s i is maximum when all edgesleaving i have the same weight, while s i = i . Using the local node entropy,we can further define the global node entropy [38–40, 81] S GN = n π ∑ i = ρ i s i , (24)where ρ i is the probability of finding the permutationΠ i (Eqs. 2 and 8). Thus, the value S GN represents aweighted average of the local determinism over all nodesof an ordinal network (see also Unakafov and Keller [18]for the definition of conditional entropy of ordinal pat-terns). In image classification tasks [40], global node en-tropy network has proven to outperform different imagequantifiers derived from gray-level co-occurrence matri-ces (GLCMs) [82, 83], a traditional technique for textureanalysis.Contrarily to permutation entropy [4] and because ofthe intrinsic connectivity constraints of ordinal networks,the global node entropy is not maximized by random data [39, 40]. For one-dimensional data, the global nodeentropy calculated from a random ordinal network is [39] S randomGN = log ( d x + ) − ( log 4 )/( d x + ) . (25)While there is no equivalent expression for two-dimensional data, it is possible to numerically estimate S randomGN using random ordinal networks numerically gen-erated [40]. In both cases, the global node entropy canbe normalized by the value of S randomGN , that is, H GN = S GN / S randomGN .In ordpy , the global node entropy function esti-mates S GN directly from data arrays or using an ordinalnetwork as returned by ordinal network . The followingcode shows simple usages of global node entropy : >>> from ordpy import global_node_entropy>>> global_node_entropy(... [1,2,3,4,5,6,7,8,9], dx=2)0.0>>> global_node_entropy(... ordinal_network([1,2,3,4,5,6,7,8,9], dx=2))0.0>>> global_node_entropy(... np.random.uniform(size=100000), dx=3)1.4988332319747597>>> global_node_entropy(... random_ordinal_network(dx=3))1.5 VII. APPLICATIONS OF ORDINALNETWORKS WITH ORDPY To better illustrate the use of ordpy in the contextof ordinal networks, we review and replicate some litera-ture results. Before starting, we remark that ordpy doesnot have functions for network analysis or graph visu-alization. The ordinal network function processes dataand generates output data (nodes, edges and weight lists)that can feed graph libraries such as graph tool [84], networkx [85], and igraph [86]. Here, we have used networkx and igraph .We start by partially reproducing Small’s [36] pioneer-ing work in which “ordinal partition networks” first ap-peared (see Fig. 3 in that work). By following Small [36],we numerically solve the differential equations of theR¨ossler system (with parameters a = . b = c = 4, see Appendix A for definitions) and sample the x -coordinate to obtain a time series with 10 observations.Figure 6a illustrates the periodic behavior of this time se-ries. We then estimate the ordinal network from this dataset with embedding parameters d x = 16 and τ x = 1. It isworth remembering that Small’s original algorithm usesnon-overlapping partitions and the edges of the result-ing ordinal network are undirected and unweighted. Theparameter overlapping in ordinal network should beequal to False to properly use Small’s original algorithm.Figure 6b shows a visualization of this ordinal network,where the circular structure alludes to the periodicity ofthe original time series.4 Time index, t x - c oo r d i na t e R ö ss l e r , x t (a) Time index, t F r a c t i ona l B r o w n i an m o t i on , x t (c) h = 0.8 (e)(b) d x = 16 (d) d x = 3 (f) d x = d y = 2 FIG. 6. Ordinal networks mapped from one- and two-dimensional data. (a) Time series obtained from the x -coordinate of theR¨ossler system (with parameters a = . b = c = observations. This timeseries exhibits a periodic behavior after a short initial transient. (b) Visualization of the ordinal network mapped from the x -coordinate of the R¨ossler system with embedding parameters d x = 16 and τ x = 1. This ordinal network uses Small’s originalalgorithm [36] with non-overlapping data partitions and undirected and unweighted edges. (c) Time series obtained from arealization of the fractional Brownian motion with h = . observations are shown). (d) Ordinalnetwork representation of the fractional Brownian motion time series with d x = τ x = 1. Here we have used overlappingdata partitions and made edge thickness proportional to edge weight. (e) Example of a periodic ornament with size 250 × d x = d y = τ x = τ y = 1. In thisvisualization, edge thickness is made proportional to edge weight. Data and code necessary to reproduce these results areavailable in a Jupyter notebook at ordpy ’s webpage. In another simple example with ordinal networks, wepartially replicate Pessa and Ribeiro’s [39] results on frac-tional Brownian motion (see Fig. 6 in their work). To doso, we generate a time series from this stochastic processwith Hurst exponent h = . observations, as illustrated in Fig. 6c.Next, we map this time series into an ordinal networkwith embedding parameters d x = τ x = h = . ( , , ) and ( , , ) (that is, the upward and down-ward trends of this time series). Pessa and Ribeiro [39]have also shown that local properties of ordinal networks(for instance, average weighted shortest path) are quiteeffective for estimating the Hurst exponent of time series,having performance superior to widely used approachessuch as detrended fluctuation analysis (DFA) [87].We also consider ordinal networks mapped from two-dimensional data. We map a periodic ornament previ- ously explored in Ref. [40] (see Fig. 2 in that reference).Figure 6e shows the ornament of size 250 × 250 (see Ap-pendix B for more details), while Fig. 6f presents a visu-alization of the corresponding ordinal network with em-bedding parameters d x = d y = τ x = τ y = 1. Wehave made edge thickness proportional to the edge weight(Eq. 22). We observe that a few edges concentrate mostof the transition probability of the network and only asmall fraction of all possible nodes and edges (24 nodesand 416 edges) appear in this network.In addition to the previous more qualitative examples,we have also replicated some results related to the globalnode entropy of ordinal networks. For time series, we fol-low Pessa and Ribeiro [39] (see Fig. 5 in their work) andgenerate a periodic sawtooth-like signal (Fig. 7a) with10 observations and add to it uniform white noise inthe interval [− ξ, ξ ] , where ξ represents the noise ampli-tude. We generate these noisy sawtooth-like time seriesfor each ξ ∈ { , . , . , . . . , } and estimate the averagevalues of the normalized permutation entropy ( H ) andthe normalized global node entropy ( H GN ) over ten timeseries replicas with d x = τ x = Time index, t T i m e s e r i e s , x t (a) = 0 Noise amplitude, E n t r op y d x = 4 (b) Permutation entropy, H Global node entropy, H GN D15 D21 D68 D70D38 D55 D56 D65 (c) S H o r i z o n t a l G N S V e r t i c a l G N D15D21 D38 D55D56 D65D68D70 (d) FIG. 7. Global node entropy of one- and two-dimensional data. (a) The initial data points of a sawtooth-like time series definedas x t = { , / , / , , . . . } . (b) Normalized permutation entropy ( H ) and normalized global node entropy ( H GN ) as a functionof the amplitude of the uniform white noise ( ξ ) added to the periodic sawtooth-like signals. The different curves representaverage values of H and H GN over ten realizations for ξ = { , . , . , . . . , } . (c) Eight examples (out of 112) of the normalizedBrodatz textures. These are grayscale images (256 gray levels) with size 640 × 640 [88]. (d) Differences between the global nodeentropy estimated from the horizontal and vertical ordinal networks ( S HorizontalGN − S VerticalGN ) mapped from each Brodatz texture.We highlight eight textures (the same shown in panel c) with the largest differences. Data and code necessary to reproducethese results are available in a Jupyter notebook at ordpy ’s webpage. Figure 7b shows the average values of H and H GN as afunction of the noise amplitude ξ . We note that bothmeasures approach one with the increase of the noiseamplitude. However, permutation entropy saturates for ξ ≈ 1, while global node entropy requires significantlyhigher values of ξ . This result indicates that global nodeentropy is more robust to noise addition and has a higherdiscrimination power than permutation entropy [39].To demonstrate the use of global node entropy withtwo-dimensional data, we calculate the global node en-tropy for a set of 112 8-bit images of natural texturesknown as the normalized Brodatz textures [88, 89]. Fig-ure 7c shows examples of these images. By follow-ing Pessa and Ribeiro [40] (see Fig. 5 in their work),we estimate the global node entropy from the horizon-tal ( S HorizontalGN ) and vertical ( S VerticalGN ) ordinal networksmapped from the Brodatz textures with d x = d y = τ x = τ y = S HorizontalGN − S VerticalGN ) for each Bro-datz texture. We have also highlighted eight textures with extreme values for this difference. Most of theseimages are characterized by stripes or line segments pre-dominantly oriented in vertical or horizontal directionswhich, in turn, suggests that properties of vertical andhorizontal ordinal networks can detect simple image sym-metries.In a final application with ordinal networks, we ex-plore the concept of missing links or missing transitionsamong ordinal patterns [39]. Similarly to the missing or-dinal patterns described by Amig´o et al. [63, 65], ordinalnetworks can display true and false forbidden transitionsamong ordinal patterns. In this case, true missing linksare related to the intrinsic dynamics of the process un-der analysis, while false missing links are associated withthe finite size of empirical data sets. Because we knowthe exact form of random ordinal networks [39, 40] (herethese networks represent all possible connections) we canreadily find all missing links of an empirical ordinal net-work. In ordpy , the missing links function evaluatesall missing ordinal transitions directly from a data set orthe returned arrays of ordinal network as in:6 Time series length, N x F r a c t i on o f m i ss i ng li n ks , f white noise (a) d x = 3 d x = 4 d x = 5 d x = 6 Time series length, N x F r a c t i on o f m i ss i ng li n ks , f logistic map (b) d x = 3 d x = 4 d x = 5 d x = 6 FIG. 8. True and false missing links in ordinal networks. (a) Dependence of the fraction of missing links ( f ) estimated fromGaussian white noise time series as a function of the time series length ( N x ). The different curves represent average values overten realizations (for each series length) and embedding dimension d x ∈ { , , , } with τ x = 1. (b) Dependence of the fractionof missing links ( f ) estimated from fully chaotic logistic time series ( r = 4) as a function of the time series length ( N x ). In bothpanels, the different curves represent average values over ten realizations (for each series length) and embedding dimension d x ∈ { , , , } with τ x = 1. We also use 194 values for N x logarithmically spaced in the interval [ , ] in both panels. Dataand code necessary to reproduce these results are available in a Jupyter notebook at ordpy ’s webpage. >>> from ordpy import missing_links>>> missing_links([4,7,9,10,6,11,3], dx=2,... return_fraction=False)(array([['1|0', '1|0']], dtype=' 1. Figure 8a shows these fractions of missing linksas a function of the time series length. We observe thatthis quantity approaches zero as N x becomes sufficientlylarge. Furthermore, the smaller the embedding dimen-sion, the faster the missing links vanish. This pattern isa fingerprint of false missing links. We have also carriedout the same analysis with time series generated from lo-gistic map iterations at fully developed chaos. Figure 8bshows the corresponding results. Unlike white noise, thelogistic map produces ordinal networks with missing linksthat persist even in considerably long time series. Thisbehavior is typical of true missing links. VIII. CONCLUSIONS We have introduced ordpy – an open-source Pythonmodule for data analysis that implements several ordi-nal methods associated with the Bandt-Pompe frame-work. Specifically, ordpy has functions implementing the following methods: permutation entropy, complexity-entropy plane, missing ordinal patterns, Tsallis andR´enyi complexity-entropy curves, ordinal networks, andmissing ordinal transitions. All ordpy ’s functions au-tomatically deal with one-dimensional (time series) andtwo-dimensional (images) data. Furthermore, most ofthese functions are also ready for multiscale analysisvia the embedding delay parameters. Along with thedescription of ordpy functionalities, we have also pre-sented a literature review of several of the principal meth-ods related to Bandt and Pompe’s framework. This re-view further includes a reproduction of several litera-ture results with ordpy ’s functions. Beyond the sum-marized description of ordpy ’s functions presented here,we notice that a complete documentation is available at arthurpessa.github.io/ordpy . All data and code usedin this work are also freely available at ordpy ’s website.We believe ordpy will help to popularize ordinal meth-ods even further, particularly in research fields with morelimited tradition in scientific computing. In addition toa myriad of possible empirical applications, we also be-lieve ordpy can further promote the development of newmethods related to Bandt and Pompe’s framework. Weremark that some techniques available in ordpy have re-ceived little attention or have not even been formallyproposed. These possible developments already imple-mented in ordpy include the use of complexity-entropycurves for two-dimensional data, multiscale complexity-entropy curves, ordinal networks with different embed-ding delays (particularly for two-dimensional data), anal-ysis of missing patterns in two-dimensional data, andmissing ordinal transitions. We also plan to implementmore techniques based on the Bandt and Pompe’s frame-work and include them in future versions of ordpy .Finally, we hope our module helps making re-7search methods more accessible and reproducible [90,91] as well as other open-source software efforts suchas the tisean [92] (nonlinear time series analysis), pyunicorn [93] (time series networks and recurrenceanalysis), and powerlaw [94] (analysis of heavy-tailed dis-tributions) packages. ACKNOWLEDGMENTS This research was supported by Coordena¸c˜ao de Aper-feicoamento de Pessoal de N´ıvel Superior (CAPES)and Conselho Nacional de Desenvolvimento Cient´ıficoe Tecnol´ogico (CNPq – Grants 407690/2018-2 and303121/2018-1). Appendix A: Definitions of dynamical systems In this appendix, we present a brief definition of thedynamical systems used in this manuscript.1. The logistic map is defined by the following differenceequation [95]: x t + = rx t ( − x t ) , (A1)where r is a parameter. We have used r = x t + = r ( t ) x t ( − x t ) , (A2)where the parameter r ( t ) changes at each iteration. Theresults of Figs. 2e and 2f were obtained with r ( t = ) = . − up to r ( t ) = { x / ω for x ∈ [ , ω ]( − x )/( − ω ) for x ∈ [ ω, ] , (A3)where ω is a parameter. In the complexity-entropy planeapplication shown in Fig. 3a, we have used ω = . { x t + = − ax t + y t y t + = bx t , (A4)where a and ∣ b ∣ < a = . b = . x t + = ( x t + x zt ) mod 1 , (A5)where z is a parameter, and the modulo operation returnsthe fractional (decimal) part of a number. We have used z ∈ { , . , } , as shown in Fig. 3a. 6. The R¨ossler system is a continuous time dynamicalsystem defined as [99, 100] dxdt = − y − zdydt = x + aydzdt = b + z ( x − c ) , (A6)where a , b , and c are parameters. We have numericallysolved this differential equation system using the scipy Python module [101] with parameters a = . b = c = Appendix B: Definitions of stochastic processes In this appendix, we briefly describe the stochastic pro-cesses used in the manuscript.1. An Ising surface [102, 103] is a square lattice in whichthe height at each lattice site represents the accumulatedsum of spin variables of particles in a Monte Carlo simu-lation [66]. If we assume σ i ∈ {− , } represents the spinvariable at site i , we can write the Hamiltonian of thissystem as H = − ∑ ⟨ i,j ⟩ σ i σ j , (B1)where the summation is over all pairs of first neighborsin a square lattice. The height S i at site i of the corre-sponding Ising surface is then defined as S i = ∑ t σ i ( t ) , (B2)where σ i ( t ) is the spin value in step t of the MonteCarlo simulation. For each surface, we define the re-duced temperature T r as the ratio between the temper-ature T and the critical temperature T c of the Ising sys-tem ( T c = / ln ( + √ )) . Finally, we have used periodicboundary conditions in our numerical experiments.2. A fractional Brownian motion is a continuous, self-similar, and non-stationary stochastic process introducedby Mandelbrot and Van Ness [104]. The Hurst exponent h ∈ ( , ) controls the roughness observed in samples ofthis process, such that the smaller the values of h , therougher the time series. The case h = / h ∈ ( , ) controls the range of auto-correlationof the time series. For h > / 2, the process presents long-range persistent memory. For h < / 2, samples present8anti-persistent behavior. We have Gaussian white noiseif h = / 2. To generate samples of a fractional Gaussiannoise, we have also used the Hosking method [105, 106].More detailed information about simulations of fractionalGaussian noise and fractional Brownian motion can befound in Ref. [106]. The C source code used in this workis publicly available in Ref. [107].4. A 1 / f noise or a Flicker noise [108, 109] is a class ofstochastic processes presenting a power-law power spec-tral density [108–110], that is, P ∼ / f − k . The case k = k = / f − k noise for k ∈ { , . , . , . . . , . } with the algorithm proposed byTimmer and K¨onig [110] as implemented in Ref. [111].5. The periodic ornament used in this work can be gen- erated by first defining two square arrays X i,j = π ( j − )( n − ) and Y i,j = π ( i − )( n − ) , (B3)and next by calculating Z i,j = sin ( ω π X i,j cos θ − ω π Y i,j sin θ ) , (B4)where i, j = , . . . , n , with n being the ornament size, θ defining the stripes angle, and ω > n = ω = θ = 135 degrees. Previousworks [35, 40, 112] have also considered shuffled versionsof this periodic ornament, where a parameter controls thefraction of elements Z i,j that are randomly shuffled. Afunction implementing this geometric ornament is avail-able in ordpy ’s notebook. [1] H. Kantz and T. Schreiber, Nonlinear Time Series Anal-ysis (Cambridge University Press, New York, 2004).[2] E. Bradley and H. Kantz, Nonlinear time-series analysisrevisited, Chaos , 097610 (2015).[3] C. E. Shannon, A mathematical theory of communica-tion, The Bell System Technical Journal , 379 (1948).[4] C. Bandt and B. Pompe, Permutation entropy: A natu-ral complexity measure for time series, Physical ReviewLetters , 174102 (2002).[5] N. Nicolaou and J. Georgiou, Detection of epileptic elec-troencephalogram based on permutation entropy andsupport vector machines, Expert Systems with Appli-cations , 202 (2012).[6] L. Zunino, M. Zanin, B. M. Tabak, D. G. P´erez, andO. A. Rosso, Forbidden patterns, permutation entropyand stock market inefficiency, Physica A , 2854(2009).[7] J. Garland, T. Jones, M. Neuder, V. Morris, J. White,and E. Bradley, Anomaly detection in paleoclimaterecords using permutation entropy, Entropy , 931(2018).[8] R. Yan, Y. Liu, and R. X. Gao, Permutation entropy:A nonlinear statistical measure for status characteriza-tion of rotary machines, Mechanical Systems and SignalProcessing , 474 (2012).[9] M. Matilla-Garc´ıa and M. Ruiz Mar´ın, A non-parametric independence test using permutation en-tropy, Journal of Econometrics , 139 (2008).[10] M. Zanin, L. Zunino, O. A. Rosso, and D. Papo, Per-mutation entropy and its main biomedical and econo-physics applications: A review, Entropy , 1553(2012).[11] M. Riedl, A. M¨uller, and N. Wessel, Practical consider-ations of permutation entropy, The European PhysicalJournal Special Topics , 249 (2013).[12] J. M. Amig´o, K. Keller, and V. A. Unakafova, Ordi-nal symbolic analysis and its application to biomedicalrecordings, Philosophical Transactions of the Royal So- ciety A , 20140091 (2015).[13] K. Keller, T. Mangold, I. Stolz, and J. Werner, Permu-tation entropy: New ideas and challenges, Entropy ,134 (2017).[14] O. A. Rosso, H. A. Larrondo, M. T. Martin, A. Plastino,and M. A. Fuentes, Distinguishing noise from chaos,Physical Review Letters , 154102 (2007).[15] L. Zunino, D. P´erez, A. Kowalski, M. Mart´ın, M. Gar-avaglia, A. Plastino, and O. Rosso, Fractional Brownianmotion, fractional Gaussian noise, and Tsallis permuta-tion entropy, Physica A , 6057 (2008).[16] L. C. Carpi, P. M. Saco, and O. Rosso, Missing ordi-nal patterns in correlated noises, Physica A , 2020(2010).[17] U. Parlitz, S. Berg, S. Luther, A. Schirdewan, J. Kurths,and N. Wessel, Classifying cardiac biosignals using ordi-nal pattern statistics and symbolic dynamics, Comput-ers in Biology and Medicine , 319 (2012).[18] A. M. Unakafov and K. Keller, Conditional entropy ofordinal patterns, Physica D , 94 (2014).[19] Z. Liang, Y. Wang, X. Sun, D. Li, L. J. Voss, J. W.Sleigh, S. Hagihira, and X. Li, EEG entropy measuresin anesthesia, Frontiers in Computational Neuroscience , 16 (2015).[20] L. Zunino, F. Olivares, and O. A. Rosso, Permutationmin-entropy: An improved quantifier for unveiling sub-tle temporal correlations, EPL (Europhysics Letters) , 10005 (2015).[21] C. Bandt, A new kind of permutation entropy used toclassify sleep stages from invisible EEG microstructure,Entropy , 197 (2017).[22] H. V. Ribeiro, M. Jauregui, L. Zunino, and E. K.Lenzi, Characterizing time series via complexity-entropy curves, Physical Review E , 062106 (2017).[23] M. Jauregui, L. Zunino, E. K. Lenzi, R. S. Mendes, andH. V. Ribeiro, Characterization of time series via R´enyicomplexity-entropy curves, Physica A , 74 (2018).[24] W. Aziz and M. Arif, Multiscale permutation entropy of physiological time series, in (IEEE, 2005) pp. 1–6.[25] L. Zunino, M. C. Soriano, I. Fischer, O. A. Rosso,and C. R. Mirasso, Permutation-information-theory ap-proach to unveil delay dynamics from time-series anal-ysis, Physical Review E , 046212 (2010).[26] F. C. Morabito, D. Labate, F. L. Foresta, A. Bra-manti, G. Morabito, and I. Palamara, Multivariatemulti-scale permutation entropy for complexity analysisof Alzheimer’s disease EEG, Entropy , 1186 (2012).[27] L. Zunino, M. C. Soriano, and O. A. Rosso, Distinguish-ing chaotic and stochastic dynamics from time series byusing a multiscale symbolic approach, Physical ReviewE , 046210 (2012).[28] B. Fadlallah, B. Chen, A. Keil, and J. Pr´ıncipe,Weighted-permutation entropy: A complexity measurefor time series incorporating amplitude information,Physical Review E , 022911 (2013).[29] J. Xia, P. Shang, J. Wang, and W. Shi, Permutation andweighted-permutation entropy analysis for the complex-ity of nonlinear time series, Communications in Nonlin-ear Science and Numerical Simulation , 60 (2016).[30] H. Azami and J. Escudero, Amplitude-aware permuta-tion entropy: Illustration in spike detection and sig-nal segmentation, Computer Methods and Programs inBiomedicine , 40 (2016).[31] S. Chen, P. Shang, and Y. Wu, Weighted multiscaleR´enyi permutation entropy of nonlinear time series,Physica A , 548 (2018).[32] C. Bian, C. Qin, Q. D. Y. Ma, and Q. Shen, Modi-fied permutation-entropy analysis of heartbeat dynam-ics, Physical Review E , 021906 (2012).[33] D. Cuesta-Frau, M. Varela-Entrecanales, A. Molina-Pic´o, and B. Vargas, Patterns with equal values in per-mutation entropy: Do they really matter for biosignalclassification?, Complexity , 1324696 (2018).[34] H. V. Ribeiro, L. Zunino, E. K. Lenzi, P. A. Santoro,and R. S. Mendes, Complexity-entropy causality planeas a complexity measure for two-dimensional patterns,PLoS One , 1 (2012).[35] L. Zunino and H. V. Ribeiro, Discriminating image tex-tures with the multiscale two-dimensional complexity-entropy causality plane, Chaos, Solitons & Fractals ,679 (2016).[36] M. Small, Complex networks from time series: Captur-ing dynamics, in (2013) pp. 2509–2512.[37] M. McCullough, M. Small, T. Stemler, and H. H.-C. Iu,Time lagged ordinal partition networks for capturingdynamics of continuous dynamical systems, Chaos ,053101 (2015).[38] M. Small, M. McCullough, and K. Sakellariou, Ordinalnetwork measures — quantifying determinism in data,in (2018) pp. 1–5.[39] A. A. B. Pessa and H. V. Ribeiro, Characterizingstochastic time series with ordinal networks, PhysicalReview E , 042304 (2019).[40] A. A. B. Pessa and H. V. Ribeiro, Mapping imagesinto ordinal networks, Physical Review E , 052312(2020).[41] J. B. Borges, H. S. Ramos, R. A. F. Mini, O. A. Rosso,A. C. Frery, and A. A. F. Loureiro, Learning and distin- guishing time series dynamics via ordinal patterns tran-sition graphs, Applied Mathematics and Computation , 124554 (2019).[42] E. T. C. Chagas, A. C. Frery, O. A. Rosso, and H. S.Ramos, Characterization of SAR images with weightedamplitude transition graphs, in (2020) pp. 264–269.[43] The data deluge, Available: (2010),Accessed: 20 Oct 2020.[44] C. A. Mattmann, A vision for data science, Nature ,473 (2013).[45] D. M. Blei and P. Smyth, Science and data science, Pro-ceedings of the National Academy of Sciences , 8689(2017).[46] C. R. Harris, K. J. Millman, S. J. van der Walt, R. Gom-mers, P. Virtanen, D. Cournapeau, E. Wieser, J. Tay-lor, S. Berg, N. J. Smith, R. Kern, M. Picus, S. Hoyer,M. H. van Kerkwijk, M. Brett, A. Haldane, J. F. del R´ıo,M. Wiebe, P. Peterson, P. G´erard-Marchant, K. Shep-pard, T. Reddy, W. Weckesser, H. Abbasi, C. Gohlke,and T. E. Oliphant, Array programming with NumPy,Nature , 357 (2020).[47] J. M. Perkel, Pick up Python, Nature , 125 (2015).[48] H. Shen, Interactive notebooks: Sharing the code, Na-ture , 151 (2014).[49] T. Kluyver, B. Ragan-Kelley, F. P´erez, B. Granger,M. Bussonnier, J. Frederic, K. Kelley, J. Hamrick,J. Grout, S. Corlay, P. Ivanov, D. Avila, S. Abdalla,C. Willing, and Jupyter development team, Jupyternotebooks – a publishing format for reproducible com-putational workflows, in Positioning and Power in Aca-demic Publishing: Players, Agents and Agendas , editedby F. Loizides and B. Scmidt (IOS Press, 2016) pp. 87–90.[50] Y. Cao, W. wen Tung, J. B. Gao, V. A. Protopopescu,and L. M. Hively, Detecting dynamical changes in timeseries using the permutation entropy, Physical ReviewE , 046217 (2004).[51] D. Cuesta-Frau, J. P. Murillo-Escobar, D. A. Orrego,and E. Delgado-Trejos, Embedded dimension and timeseries length. Practical influence on permutation en-tropy and its applications, Entropy , 385 (2019).[52] L. Zunino, F. Olivares, F. Scholkmann, and O. A. Rosso,Permutation entropy based time series analysis: Equal-ities in the input signal can lead to false conclusions,Physics Letters A , 1883 (2017).[53] J. Amig´o, S. Zambrano, and M. A. Sanju´an, Combina-torial detection of determinism in noisy time series, EPL(Europhysics Letters) , 60005 (2008).[54] O. A. Rosso and C. Masoller, Detecting and quanti-fying temporal correlations in stochastic resonance viainformation theory measures, The European PhysicalJournal B , 37 (2009).[55] L. Zunino, M. Zanin, B. M. Tabak, D. G. P´erez, andO. A. Rosso, Complexity-entropy causality plane: Auseful approach to quantify the stock market ineffi-ciency, Physica A , 1891 (2010).[56] L. Zunino, A. Fern´andez Bariviera, M. B. Guercio,L. B. Martinez, and O. A. Rosso, On the efficiency ofsovereign bond markets, Physica A , 4342 (2012).[57] H. V. Ribeiro, L. Zunino, R. S. Mendes, and E. K. Lenzi,Complexity–entropy causality plane: A useful approach for distinguishing songs, Physica A , 2421 (2012).[58] H. Y. D. Sigaki, R. F. de Souza, R. T. de Souza, R. S.Zola, and H. V. Ribeiro, Estimating physical proper-ties from liquid crystal textures via machine learningand complexity-entropy methods, Physical Review E , 013311 (2019).[59] R. L´opez-Ruiz, H. L. Mancini, and X. Calbet, A statis-tical measure of complexity, Physics Letters A , 321(1995).[60] J. Lin, Divergence measures based on the Shannon en-tropy, IEEE Transactions on Information Theory ,145 (1991).[61] P. W. Lamberti, M. T. Martin, A. Plastino, and O. A.Rosso, Intensive entropic non-triviality measure, Phys-ica A , 119 (2004).[62] M. T. Martin, A. Plastino, and O. A. Rosso, General-ized statistical complexity measures: Geometrical andanalytical properties, Physica A , 439 (2006).[63] J. M. Amig´o, L. Kocarev, and J. Szczepanski, Orderpatterns and chaos, Physics Letters A , 27 (2006).[64] C. Bandt and F. Shiha, Order patterns in time series,Journal of Time Series Analysis , 646 (2007).[65] J. M. Amig´o, S. Zambrano, and M. A. Sanju´an, Trueand false forbidden patterns in deterministic and ran-dom dynamics, EPL (Europhysics Letters) , 50001(2007).[66] D. P. Landau and K. Binder, A Guide to Monte CarloSimulations in Statistical Physics (Cambridge Univer-sity Press, New York, 2009).[67] H. Y. D. Sigaki, M. Perc, and H. V. Ribeiro, History ofart paintings through the lens of entropy and complex-ity, Proceedings of the National Academy of Sciences , E8585 (2018).[68] Abstract Painting: Blue, Available: , Accessed: 24 Nov2020.[69] Abaporu, Available: , Accessed: 24 Nov2020.[70] Number 1 (Lavender Mist), Available: , Accessed: 24Nov 2020.[71] M. Zanin, Forbidden patterns in financial time series,Chaos , 013119 (2008).[72] K. Sakellariou, M. McCullough, T. Stemler, andM. Small, Counting forbidden patterns in irregularlysampled time series. II. Reliability in the presence ofhighly irregular sampling, Chaos , 123104 (2016).[73] M. McCullough, K. Sakellariou, T. Stemler, andM. Small, Counting forbidden patterns in irregularlysampled time series. I. The effects of under-sampling,random depletion, and timing jitter, Chaos , 123103(2016).[74] C. W. Kulp, J. M. Chobot, B. J. Niskala, and C. J.Needhammer, Using forbidden ordinal patterns to de-tect determinism in irregularly sampled time series,Chaos , 023107 (2016).[75] C. Tsallis, Possible generalization of Boltzmann-Gibbsstatistics, Journal of Statistical Physics , 479 (1988).[76] A. R´enyi, On measures of entropy and information, in Proceedings of the Fourth Berkeley Symposium on Math-ematical Statistics and Probability, Volume 1: Contribu- tions to the Theory of Statistics (University of CaliforniaPress, Berkeley, 1961) pp. 547–561.[77] C. Tsallis, Introduction to Nonextensive Statistical Me-chanics: Approaching a Complex World (Springer, NewYork, 2009).[78] T. van Erven and P. Harremos, R´enyi divergence andKullback-Leibler divergence, IEEE Transactions on In-formation Theory , 3797 (2014).[79] Y. Zou, R. V. Donner, N. Marwan, J. F. Donges, andJ. Kurths, Complex network approaches to nonlineartime series analysis, Physics Reports , 1 (2019).[80] M. Newman, Networks: An Introduction (Oxford Uni-versity Press, New York, 2010).[81] M. McCullough, M. Small, H. H. C. Iu, and T. Stemler,Multiscale ordinal network analysis of human cardiacdynamics, Philosophical Transactions of the Royal So-ciety A , 20160292 (2017).[82] R. M. Haralick, K. Shanmugam, and I. H. Dinstein, Tex-tural features for image classification, IEEE Transac-tions on Systems, Man, and Cybernetics , 610 (1973).[83] R. M. Haralick, Statistical and structural approaches totexture, Proceedings of the IEEE , 786 (1979).[84] T. P. Peixoto, The graph-tool Python library, figshare10.6084/m9.figshare.1164194 (2014).[85] A. A. Hagberg, D. A. Schult, and P. J. Swart, Explor-ing network structure, dynamics, and function usingnetworkX, in Proceedings of the 7th Python in ScienceConference , edited by G. Varoquaux, T. Vaught, andJ. Millman (Pasadena, 2008) pp. 11–15.[86] G. Csardi and T. Nepusz, The igraph software packagefor complex network research, InterJournal ComplexSystems , 1695 (2006).[87] C.-K. Peng, S. V. Buldyrev, S. Havlin, M. Simons, H. E.Stanley, and A. L. Goldberger, Mosaic organization ofDNA nucleotides, Physical Review E , 1685 (1994).[88] Centre for Research and Applications in Re-mote Sensing (CARTEL, University of Sher-brooke), Multiband texture database, Available: https://multibandtexture.recherche.usherbrooke.ca/normalized_brodatz.html , Accessed: 18 Dez 2020.[89] A. Safia and D.-C. He, New Brodatz-based imagedatabases for grayscale color and multiband textureanalysis, ISRN Machine Vision , 14 (2013).[90] M. Baker, 1,500 scientists lift the lid on reproducibility,Nature , 452 (2016).[91] D. Fanelli, Opinion: Is science really facing a repro-ducibility crisis, and do we need it to?, Proceedings ofthe National Academy of Sciences , 2628 (2018).[92] R. Hegger, H. Kantz, and T. Schreiber, Practical im-plementation of nonlinear time series methods: TheTISEAN package, Chaos , 413 (1999).[93] J. F. Donges, J. Heitzig, B. Beronov, M. Wiedermann,J. Runge, Q. Y. Feng, L. Tupikina, V. Stolbova, R. V.Donner, N. Marwan, H. A. Dijkstra, and J. Kurths, Uni-fied functional network and nonlinear time series analy-sis for complex systems science: The pyunicorn package,Chaos , 113101 (2015).[94] J. Alstott, E. Bullmore, and D. Plenz, powerlaw: APython package for analysis of heavy-tailed distribu-tions, PLoS One , 1 (2014).[95] R. M. May, Simple mathematical models with very com-plicated dynamics, Nature , 459 (1976).[96] L. L. Trulla, A. Giuliani, J. P. Zbilut, and C. L. Webber,Recurrence quantification analysis of the logistic equa- tion with transients, Physics Letters A , 255 (1996).[97] H. Sakai and H. Tokumaru, Autocorrelations of a cer-tain chaos, IEEE Transactions on Acoustics, Speech,and Signal Processing , 588 (1980).[98] H. G. Schuster and W. Just, Deterministic Chaos: AnIntroduction (Wiley, Weinheim, 2005).[99] O. R¨ossler, An equation for continuous chaos, PhysicsLetters A , 397 (1976).[100] S. H. Strogatz, Nonlinear Dynamics and Chaos: WithApplications to Physics, Biology, Chemistry, and Engi-neering (Perseus Books, New York, 1994).[101] P. Virtanen, R. Gommers, T. E. Oliphant, M. Haber-land, T. Reddy, D. Cournapeau, E. Burovski, P. Pe-terson, W. Weckesser, J. Bright, S. J. van der Walt,M. Brett, J. Wilson, K. Jarrod Millman, N. Mayorov,A. R. J. Nelson, E. Jones, R. Kern, E. Larson, C. Carey,˙I. Polat, Y. Feng, E. W. Moore, J. Vand erPlas, D. Lax-alde, J. Perktold, R. Cimrman, I. Henriksen, E. A. Quin-tero, C. R. Harris, A. M. Archibald, A. H. Ribeiro,F. Pedregosa, P. van Mulbregt, and SciPy 1.0 Contrib-utors, Scipy 1.0: Fundamental algorithms for scientificcomputing in Python, Nature Methods , 261 (2020).[102] A. F. Brito, J. A. Redinz, and J. A. Plascak, Dynamicsof rough surfaces generated by two-dimensional latticespin models, Physical Review E , 046106 (2007).[103] A. F. Brito, J. A. Redinz, and J. A. Plascak, Two-dimensional XY and clock models studied via the dy-namics generated by rough surfaces, Physical Review E , 031130 (2010).[104] B. B. Mandelbrot and J. W. V. Ness, Fractional Brow-nian motions, fractional noises and applications, SIAMReview , 422 (1968).[105] J. R. M. Hosking, Modeling persistence in hydrologi-cal time series using fractional differencing, Water Re-sources Research , 1898 (1984).[106] T. Diecker, Simulation of fractional Brownian motion ,Ph.D. thesis, University of Twente (2004).[107] T. Diecker, hosking.c, Available: , Accessed:26 Nov 2020.[108] R. F. Voss, 1/f (flicker) noise: A brief review, in (1979) pp.40–46.[109] N. J. Kasdin, Discrete simulation of colored noise andstochastic processes and 1 / f α power law noise genera-tion, Proceedings of the IEEE , 802 (1995).[110] J. Timmer and M. K¨onig, On generating power lawnoise, Astronomy and Astrophysics , 707 (1995).[111] F. Patzelt, colorednoise.py, Available: https://github.com/felixpatzelt/colorednoise , Accessed:26 Nov 2020.[112] A. Brazhe, Shearlet-based measures of entropy and com-plexity for two-dimensional patterns, Physical Review E97