A study of visibility graphs for time series representations

Universitat Politècnica de Catalunya Facultat d’Informàtica de Barcelona Final Degree Project A study of visibility graphs for time series representations Carlos Bergillos Varela Bachelor Degree in Informatics Engineering Specialization in Computing supervised by Argimiro Arratia Quesada Computer Science Department - UPC June 22, 2020 Abstract [EN] In this project we study and present the field of visibility graphs as a way to represent and characterize time series as introduced by Lucas Lacasa et al. in 2008. With visibility graphs, time series are converted into graphs that inherit in their topology some of the structural properties of the time series, allowing the novel analysis of time series via graph theory and complex network theory tools. We study and analyze the computational complexity of a selection of different visibility graph algorithms and subsequently develop and present the ts2vg Python package that provides efficient and easy-to-use implementations for some of these visibility graph algorithms. Additionally, in this project we explore and discuss some of the potential practical applications of these visibility graphs. [ES] En este proyecto estudiamos y presentamos el campo de los grafos de visibilidad (visibility graphs) como una forma de representar y caracterizar series temporales tal y como fue introducido por Lucas Lacasa et al. en 2008. Con los grafos de visibilidad las series temporales se convierten en grafos que heredan en su topología algunas de las propiedades estructurales de las series temporales, permitiendo así el novedoso análisis de series temporales a través de herramientas propias de la teoría de grafos y la teoría de redes complejas. Estudiamos y analizamos la complejidad computacional de algunos algoritmos para obtener grafos de visibilidad y, posteriormente, desarrollamos y presentamos el paquete de Python ts2vg que proporciona implementaciones eficientes y fáciles de usar para algunos de estos algoritmos. Además, en el proyecto exploramos y discutimos algunas de las posibles aplicaciones prácticas de estos grafos de visibilidad. [CA] En aquest projecte estudiem i presentem el camp dels grafs de visibilitat (visibility graphs) com una forma de representar i caracteritzar sèries temporals tal com va ser introduït per Lucas Lacasa et al. el 2008. Amb els grafs de visibilitat les sèries temporals es converteixen en grafs que hereten en la seva topologia algunes de les propietats estructurals de les sèries temporals, permetent així el nou anàlisi de sèries temporals a través d’eines pròpies de la teoria de grafs i la teoria de xarxes complexes. Estudiem i analitzem la complexitat computacional d’alguns algoritmes per obtenir grafs de visibilitat i, posteriorment, desenvolupem i presentem el paquet de Python ts2vg que proporciona implementacions eficients i fàcils d’usar per a alguns d’aquests algoritmes. A més, al projecte explorem i discutim algunes de les possibles aplicacions pràctiques d’aquests grafs de visibilitat. Keywords— visibility graphs, time series, graph theory, complex networks, Python Contents 1 Introduction 3 1.1 Context .......................................... 3 1.2 Descriptionoftheproject ................................ 3 1.2.1 Goals and scope of the project . . . . . . . . . . . . . . . . . . . . . . . . . 3 1.3 Motivations........................................ 4 1.4 State of the art and alternatives . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 2 Definitions and Theoretical Background 6 2.1 Timeseries ........................................ 6 2.1.1 Wiener process (Brownian motion) . . . . . . . . . . . . . . . . . . . . . . . 7 2.1.2 Fractional Brownian motion . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 2.2 Graphs........................................... 9 2.2.1 Types ....................................... 9 2.2.2 Measures ..................................... 9 2.3 Real-worldnetworks ................................... 11 2.4 Graphclustering ..................................... 12 2.5 Visibilitygraphs ..................................... 12 2.5.1 Horizontal visibility graphs . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 2.6 Powerlawfitting ..................................... 15 3 Algorithms and Complexity 16 3.1 Naivestrategy ...................................... 17 3.2 Slopestrategy....................................... 19 3.3 Divide-and-conquer strategy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 3.4 Streamingalgorithm ................................... 26 4 Technologies and Implementation 27 4.1 Programminglanguage.................................. 27 4.2 Time series implementation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 4.3 Visibility graph algorithms implementation . . . . . . . . . . . . . . . . . . . . . . 28 4.3.1 Backend optimizations with C and Cython . . . . . . . . . . . . . . . . . . 28 4.3.2 Iterative algorithm for recursion . . . . . . . . . . . . . . . . . . . . . . . . 29 4.3.3 Performanceanalysis............................... 29 4.4 ts2vgpackage....................................... 29 4.4.1 Potential future considerations . . . . . . . . . . . . . . . . . . . . . . . . . 30 1 Contents 4.5 Auxiliaryfunctionalities ................................. 31 4.5.1 Graph analysis functionalities . . . . . . . . . . . . . . . . . . . . . . . . . . 31 4.5.2 Generating Wiener process and fractional Brownian motion series . . . . . 31 4.5.3 Estimating power-law exponents . . . . . . . . . . . . . . . . . . . . . . . . 31 5 Applications and Experiments 32 5.1 Performance evaluation of the different algorithms . . . . . . . . . . . . . . . . . . 32 5.2 Periodic series result in regular visibility graphs . . . . . . . . . . . . . . . . . . . . 33 5.3 Fractal series result in scale-free visibility graphs . . . . . . . . . . . . . . . . . . . 35 5.4 Estimating the Hurst exponent via the visibility graph . . . . . . . . . . . . . . . . 36 5.5 Segmentation of time series via graph clustering . . . . . . . . . . . . . . . . . . . . 38 6 Discussion and Future Work 40 Bibliography 42 Appendices 46 A Conditions and Specifications 47 B Algorithm’s Source Code 48 B.1 Naive strategy implementation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49 B.2 Slope strategy implementation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50 B.3 Divide-and-conquer strategy implementation . . . . . . . . . . . . . . . . . . . . . 51 B.4 Streaming algorithm implementation . . . . . . . . . . . . . . . . . . . . . . . . . . 52 B.5 ts2vgpackage....................................... 53 C Project Management 54 C.1 Tasks ........................................... 54 C.2 Timeplanning ...................................... 56 C.3 Budget........................................... 58 C.3.1 Humanresources................................. 58 C.3.2 Hardwareresources................................ 58 C.3.3 Softwareresources ................................ 59 C.3.4 Indirectcosts................................... 59 C.3.5 Contingencies................................... 59 C.3.6 Unexpecteddeviations.............................. 59 C.3.7 Totalbudget ................................... 60 C.4 Economicalimpact.................................... 60 C.5 Environmentalimpact .................................. 60 C.6 Socialimpact....................................... 61 C.7 Regulations ........................................ 61 C.8 Technicalcompetences.................................. 61 2 1Introduction 1.1 Context This work is developed by Carlos Bergillos Varela as a Final Degree Project for the Bachelor Degree in Informatics Engineering at the Barcelona School of Informatics (FIB) of the Universitat Politècnica de Catalunya (UPC). The director of the project is Argimiro A. Arratia Quesada, Professor from the Department of Computer Science at UPC. The project is framed within the Computing specialization of the mentioned degree, and as such, aims to display and make use of technical competencies associated with this specialization. 1.2 Description of the project Time series are an extensively used and simple data representation (in most cases they can be thought of as just an ordered list of numbers) for which a new frame of study has recently emerged, with what is known as visibility graphs. In this project we will explore the techniques introduced by Lucas Lacasa et al. 2008 [30], which allow the study of time series through a completely new perspective aided by graph theory and complex network theory ideas. 1.2.1 Goals and scope of the project The goals of the project are multiple and can be grouped and summarized as: •to serve as an opportunity for the author to research, study, and later present the field of visibility graphs, its algorithms, its potential applications, and associated complex network theory ideas •to study and define the computational complexity, both in terms of time and memory, of a selection of different algorithmic strategies useful for obtaining and working with visibility graphs •to implement correct and efficient versions of the different studied algorithms with the endgoal of developing an easy-to-use and efficient software framework for visibility graphs •to reproduce known experiments and applications of visibility graphs using the implemented algorithms and assess the obtained results and behaviors •to explore and document new possibilities and applications in the field of visibility graphs 3 Chapter 1 Introduction 1.4 1.3 Motivations Time series (or signals) are simple data representations used extensively throughout lots of different human and natural phenomena, found in all kinds of fields like economics, physics, chemistry, biology, engineering, etc. Thus, time series are an extremely simple representation idea but used to capture what sometimes are very complex systems. The field of visibility graphs allows transforming time series data into graphs via an abstract mapping, thus building a bridge between time series analysis and graph theory [43]. Time series analysis, although a mature field, has some limitations when it comes to studying complex nonlinear systems and behaviors. Graph theory is also a very mature field in mathematics and computer science with many decades of study, with hundreds of associated studies, theories, and algorithms. With this, visibility graphs now open many new possibilities by allowing to study time series from a completely different and powerful nonlinear perspective, whose full potential is still not fully known and is worth exploring. More traditional time series analysis techniques should and will continue to be used in many applications, as they have proven to be very effective for many specific tasks, and perhaps more importantly, are more mature and mathematically well-founded. Time series analysis via visibility graphs should therefore obviously not be a replacement for other well established and more direct time series analysis techniques when these techniques already produce successful results. But research and exploration of the mathematical properties and potential applications of visibility graphs is certainly worth being continued and could result in benefits not yet discovered, by potentially providing new tools effective at solving problems that might be very hard to approach otherwise. With this project we intend to gain significant knowledge on the field of visibility graphs and their associated algorithms, and to study some of their useful properties and practical applications. We will also focus on studying and evaluating the computational complexity of different algorithms used to obtain visibility graphs from an input time series, worth taking into consideration when working with algorithms that will be in charge of processing potentially very large amounts of data. This project also emerges as a response for the need to have available tools to compute and work with visibility graphs, as no adequate open tools that fulfill all the potential requirements in time series research have been found. The algorithm implementations that result from this project can be a useful tool for future research in the field. 1.4 State of the art and alternatives The study of visibility graphs applied to time series analysis is a relatively recent field and a decent amount of studies have already been presented that attempt to establish mathematical foundations for visibility graphs and their use [30, 31, 35, 43], but further work on strengthening these mathematical ideas and building new bridges with other areas of mathematics is still required. Due to their novelness, many researchers have also started to explore the use of visibility graphs applied to time series analysis in a wide variety of many different scientific fields and applications with different levels of success. For example, the use of visibility graphs has been explored in medical applications [48, 56], the study of earthquakes series [28], periodicity of series [45], forecasting of series [9], air quality patterns in cities [55], turbulent fluid flows [25, 47], hydrogeological patterns [27], sunspots and solar wind activity [49, 57], financial series [8, 42], image processing [26], and even human behavior [15] and touristic activity [3], among others. Computing and working with visibility graphs presents a unique and interesting computational challenge, and despite the amount of research on visibility graphs presented, research on the computational aspects and efficiency of the algorithms used for visibility graphs work is rarer. 4 Chapter 1 Introduction 1.4 Due to the apparent demand for visibility graphs research, the availability of efficient and accessible software implementations of the algorithms prove to be of big importance for researchers. Some available functions and implementations of visibility graphs for different programming languages already exist, such as the one provided by Lucas Lacasa in the Fortran programming language [29], and other implemented functions for Python [19] and MATLAB [24]. These implementations might be limited in terms of usability, performance, and future extensions. In this project we have developed and presented a software library and implementation that attempts to improve the situation by providing an open, accessible, easy-to-use and very efficient tool to obtain and work with visibility graphs in the Python programming language. 5 2Definitions and Theoretical Background 2.1 Time series A time series is a sequence of data points listed in time order, usually taken in regular equally spaced points in time. In other words, a time series is any data sequence taken in discrete time. The data could, for example, come from real-life measurements (like the temperature at a certain location) or be derived from a mathematical formula or model (e.g. see Brownian Motion in Section 2.1.1). Time series are commonly visualized using line plots, area plots, or bar plots where the x axis corresponds to time, and the height (yaxis) corresponds to the value of the data for the corresponding point in time, see Figure 2.1 for an example. We can mathematically formalize any time series Tusing a set of pairs: T=(t1, y1),(t2, y2),...,(tn, yn) where each pair (ti, yi)corresponds to one observation (or data point) in which yiis the data value associated to the point in time ti, and such that t1< t2, . . . , < tn. For simplicity, we will sometimes refer to a given observation (ti, yi)as just data point i. We will use the term length to refer to the number of observations of a given time series. 2015 2016 2017 2018 2019 2020 Date (year) 8000 9000 10000 11000 12000 IBEX 35 index Figure 2.1: Example of a time series plot. Values of the IBEX 35 stock market index from the year 2015 until the beginning of 2020. Data provided by ICE Data Services and accessed via Yahoo Finance [53, 54]. 6 Chapter 2 Definitions and Theoretical Background 2.1 Additionally, if the series observations are equally spaced in time, we can define a time series as just an ordered list of values: T= (y1, y2, . . . , yn) assuming that any relevant relative time information is implicitly accounted for in the order of the elements. 2.1.1 Wiener process (Brownian motion) In addition to studying real-world time series, it is also very useful to explore some known mathematical series that might help understand other more complex real-world behaviors. One such mathematical series is the one resulting from the widely studied Wiener process, also known as Brownian motion in physics and other areas of science. The Wiener process has large practical and theoretical significance, it can be used as a model for many different real-world phenomena. For example, it can be used to model the movement of physical particles or to model the values of the stock market [14]. The Wiener process W(t),t∈R,t≥0, is a mathematical continuous-time stochastic process with the following properties [32]: 1.W(0) = 0 2.W(t)is a continuous function in t 3.Whas independent increments For 0≤s<t<u<v, the increments W(t)−W(s)and W(v)−W(u)are independent. 4.Whas Gaussian increments For 0≤s<t, the increment W(t)−W(s)is normally distributed with mean 0 and variance t−s:W(t)−W(s)∼ N(0, t −s). The Wiener process as defined is a continuous-time signal. We can easily construct a discretetime version Wt,t∈N, using an iterative random-walk approach defined as: W0= 0 Wi=Wi−1+ξi (2.1) Where ξ1, ξ2, . . . are independent and identically distributed (i.i.d.) random variables with mean 0 and variance 1. So, for any i:ξi∼ N(0,1). When n→ ∞,Wnbehaves similarly to the continuous Wiener process described before [32]. From now on, we will use the name Wiener process to refer to this discrete version. With this method we can construct time series as formalized in Section 2.1: TW=(1, W1),(2, W2),...,(n, Wn) This iterative approach can be easily implemented programmatically, thus allowing us to generate Wiener process time series of any length nwe want (only limited by practical computer limitations). Because Wtis a stochastic process (the stochastic nature emerges from the use of randomized increments) we can generate infinitely many different time series TW. Another way to look at Wiener process series is that they are constructed by integrating a Gaussian series with distribution N(0,1), integrating what is known as Gaussian white noise (see Figure 2.2 for an example). 7 Chapter 2 Definitions and Theoretical Background 2.6 11 12 13 14 15 16 17 18 19 20 21 22 (a) 1 2 3 4 5 6 7 8 9 10 11 12 (b) 1 2 3 4 5 6 7 8 9 10 11 12 (c) Figure 2.7: Visibility graphs are invariant to any affine transformations of the time series data. For example, and using the time series in fig. 2.6 as the source, the same visibility graph is maintained after horizontal translation (a) vertical rescaling (b) and addition of a linear trend (c) to the original data. Lossy The affine transformation invariance just described also implies that the visibility graph mapping is not bijective, and two different time series can have the same visibility graph. Thus, we are losing information in the conversion (information such as magnitudes and absolute values of the time series) and recovering the original time series from its visibility graph is not possible. 2.5.1 Horizontal visibility graphs Visibility graphs as described in Section 2.5 (which are also sometimes known as natural visibility graphs to avoid ambiguity) are not the only way to map time series into graphs. Another of such techniques is the one known as horizontal visibility graphs [35, 43]. In this case, a graph of nnodes is constructed (again, corresponding to the nobservations of the input time series), and any two of such nodes are connected if and only if there is obstruction-free horizontal visibility between the two data points (see an illustrated example in Figure 2.8). Mathematically, two observations (ta, ya)and (tb, yb), with ta< tb, are connected via an edge (a, b)in a horizontal visibility graph if and only if: ya, yb> ycfor all csuch that ta< tc< tb(2.14) If two nodes have horizontal visibility they will also have direct natural visibility (but not the other way around), therefore the horizontal visibility graph of a given time series will always be a subgraph of its natural visibility graph. 1 2 3 4 5 6 7 8 9 10 11 12 0.0 0.2 0.4 0.6 0.8 → 1 2 3 4 5 6 7 8 9 10 11 12 → Figure 2.8: Bar plot for an example time series (left), horizontal obstruction-free lines traced (center), and resulting horizontal visibility graph (right). The same example time series is used as in Figure 2.6. Notice that the horizontal visibility graph is a subgraph of the previous visibility graph. 14 Chapter 2 Definitions and Theoretical Background 2.6 2.6 Power law fitting As previously stated, a power-law distribution is one that follows the form: P(k) = ck−γ(2.15) The constant γis called the exponent of the power law and characterizes the slope of the distribution in a log-log plot. The constant cis not as interesting to work with, since, once γis fixed, cis determined by the requirement that the distribution P(x)sums to 1. When given a set of observations x1, . . . , xn(for example, the sequence of degrees of the n nodes of a graph) and believing that its probability distribution follows a power law (following a straight line in a log-log plot as a necessary but not sufficient condition) we can use various techniques to try to estimate the value for γ. This is, finding the value for the exponent γthat results in the ideal power-law distribution that most accurately models the distribution of our empirical observations. Note that in real data the power-law behavior often only manifests in the tail of the distributions, this is, only for values of xlarger than some xmin. One of the most reliable methods to obtain γis to use Maximum Likelihood Estimation (MLE), in this case given by the following formula [10, 41]: γ= 1 + n"n X i=1 log xi xmin #−1 (2.16) where xi,i= 1, . . . , n are the observed values of xsuch that xi≥xmin. 15 3Algorithms and Complexity We explored different algorithms to build visibility graphs. The input to all the algorithms should be a time series as described in Section 2.1, and they should return an undirected and unweighted graph complying with the visibility graph definition given in Section 2.5. Therefore, for all of the upcoming algorithms (unless otherwise specified) the input and output will be of the following form: Input Time series T:T=(t1, y1),(t2, y2),...,(tn, yn) Output Graph G:G= (V, E) As a quick summary, the computational cost in time of three algorithms is provided in Table 3.1. Strategy Best Case Worst Case Average Case Naive O(n2)O(n3)O(n3) Slope O(n2)O(n2)O(n2) D&C O(nlog n)O(n2)O(nlog n) Table 3.1: Asymptotic computational time complexity of three different later described algorithms to obtain the visibility graph from an input time series of length n. These algorithms are described in detail in the upcoming sections in this chapter. 16 Chapter 3 Algorithms and Complexity 3.1 3.1 Naive strategy Algorithm 1 Visibility Graph - Naive Strategy 1: function VisibilityGraphT=(t1, y1),...,(tn, yn) 2: V←{1, . . . , n} 3: E←∅ 4: for each (ta, ya),(tb, yb)∈Twith ta< tbdo 5: if NoObstructions(T,a,b)then 6: add edge (a, b)to E 7: return G= (V, E) 8: end function 9: 10: function NoObstructionsT,a,b 11: for each (tc, yc)∈Twith ta< tc< tbdo 12: if yc≥yb+ (ya−yb)tb−tc tb−tathen 13: return False 14: return True 15: end function A naive and straightforward approach to obtain the visibility graph consists on exhaustively traversing all pairs of nodes (all potential edges), and for each one, check if any obstruction is found between the two corresponding data points by checking if the inequality (2.13), in page 13, is broken by any in-between point and add the corresponding edge when appropriate. For a given pair of nodes, as soon as one obstructing in-between data point is found (see Figure 3.1 for an example) we can stop checking the rest of the in-between points since at that moment we are already certain that no edge should be present between the two nodes. On the other hand, for two nodes to be connected in the visibility graph, we need to guarantee that all in-between points have been checked and no obstruction has been found. 1234567 t Figure 3.1: Nodes corresponding to observations at t1and t6in this example will not be connected in their visibility graph because the data point at t4lays in the obstruction shaded area. Correctness The algorithm exhaustively visits all potential edges in the graph and checks inequality (2.13) (in page 13) for all of them, therefore directly complying with the visibility graph definition given in Section 2.5. Time complexity For any given time series of size nand its resulting graph with nnodes, we need to check all potential edges in the graph, and to do so, we are traversing all possible pairs of data points. We know that the number of such potential edges (and the number of such pairs 17 Chapter 3 Algorithms and Complexity 3.1 of data points) is n 2=n(n−1) 2. This results in n(n−1) 2iterations (O(n2)) of the external for-loop found in line 4. For each iteration of this loop with a given potential edge (a, b)we need to check whether all intermediate points c(with ta< tc< tb) satisfy the required inequality. Note that, and as explained before, the number of inequality checks that need to be performed is not fixed and depends on the particular structure of the time series (as soon as one intermediate point is found to not satisfy the inequality we can stop checking for that edge). In the worst case, for all pairs (a, b), all intermediate points cwill need to be checked (b−a−1intermediate points in each case, in the order of O(n)asymptotically), resulting in a total algorithm time complexity of O(n3). On the other hand, in the best case (e.g. in a flat time series, for this algorithm), each pair of nodes gets checked in constant time, as for all potential edges (a, b), either aand bare immediate neighbors (without intermediate points c, and therefore connected), or aand bare not immediate neighbors but the first intermediate point cthat gets checked will already present an obstruction not satisfying the inequality. In such time series checking for any given potential edge (a, b)will always require constant time, resulting in a total algorithm time complexity of O(n2)in this best case. The average time complexity of this algorithm is difficult to estimate, as there are infinitely many possible input time series, and the exact running time behavior depends on the structural properties of the input series. With that said, we can make the assumption that on average, for each potential edge (a, b), an unknown fixed fraction of the b−a−1intermediate points will get checked, and therefore the average time complexity does not change asymptotically, remaining at O(n3). Space complexity The algorithm does not require any additional data structures apart from the input time series of size n, and the output graph consisting of nnodes and, at most, n(n−1) 2edges, O(n2)asymptotically. Therefore the cost of the space required in memory for this algorithm will be dictated by the size of the output graph, in the order of O(n2). 18 Chapter 3 Algorithms and Complexity 3.2 3.2 Slope strategy Algorithm 2 Visibility Graph - Slope Strategy 1: function VisibilityGraphT=(t1, y1),...,(tn, yn) 2: V←{1, . . . , n} 3: E←∅ 4: for a←1to n−1do 5: slopemax ←−∞ 6: for b←a+ 1 to ndo 7: slopea,b ←(yb−ya)/(tb−ta) 8: if slopea,b >slopemax then 9: add edge (a, b)to E 10: slopemax ←slopea,b 11: return G= (V, E) 12: end function We can improve the previously described naive strategy by benefiting from information provided by the slope of the visibility lines as explained in [33]. Mathematically, for a time series T, a pair of data points a, b ∈T(with ta< tb) are connected nodes in the visibility graph of T if and only if the slope of the visibility line from ato bis larger than the maximum slope found between aand all intermediate data points cwith ta< tc< tb(see proof of this in the upcoming Theorem 3.2.1). This way we can find all edges for the node corresponding to time t1in a single pass through the other points by just storing the value for the maximum slope found so far (see Figure 3.2 for an example), and similarly for the remaining edges for the node at t2, etc. 1234567 t Figure 3.2: The slope of the visibility line is strictly increasing for the connected points as t increases. Correctness In order to prove that this approach complies with the visibility graph definition proposed in Section 2.5, we will first prove that the slope idea is more than just an intuition. To do so, we will define and prove the following: Theorem 3.2.1. For a given time series Tand (ta, ya),(tb, yb)∈Twith ta< tb, the edge (a, b)belongs to the visibility graph GTof Tif and only if the slope of the visibility line from ato bis larger than the maximum slope found between a and all other observations cwith ta< tc< tb. Proof. Let C0be the set of the intermediate points: C0={c|(tc, yc)∈Tand ta< tc< tb}. Let s(i, j)be the slope of the line that passes through some given data 19 Chapter 3 Algorithms and Complexity 3.2 points (ti, yi)and (tj, yj), calculated as: s(i, j) = yj−yi tj−ti (3.1) With this, and using the visibility graph definition and its inequality (2.12) (in page 13) we can obtain: (a, b)is an edge in the visibility graph ⇐⇒ ∀c∈C0:yc<yb−ya tb−ta (tc−ta) + ya ⇐⇒ ∀c∈C0:yc< s(a, b)·(tc−ta) + ya ⇐⇒ ∀c∈C0:yc−ya tc−ta < s(a, b) ⇐⇒ ∀c∈C0:s(a, c)< s(a, b) ⇐⇒ max c∈C0s(a, c)< s(a, b) By storing and updating the value for largest slope encountered so far (slopemax) when linearly traversing the data points to the right of a given data point a, we know that if the slope from ato a new data point bis larger than slopemax then the slope from ato bis larger than all the slopes from ato the previous points between aand band therefore there is direct line-of-sight visibility from ato band an edge (a, b)is added to the graph. In doing this linear pass starting from all possible nodes a, we will have checked the visibility for all potential edges (all pairs of nodes) in the graph. Time complexity We can quickly see that this algorithm makes use of two for-loops, one nested inside the other (lines 4 and 6). Both of these for-loops do a number of iterations in the order of O(n), resulting in a total number of iterations of the inner loop of O(n2). More precisely, the number of times that the block of code (of constant time) found inside the inner loop (lines 7 to 10) is executed is described by the following formula: n−1 X a=1 n X b=a+1 1 = n−1 X a=1 (n−a) = (n−1)n− n−1 X a=1 a = (n−1)n−(n−1)(n) 2 =n2−n 2 Which, not surprisingly, also corresponds to the number of potential edges in a graph with n nodes. This algorithm has a running time that only depends on the size of the input nand does not depend on the structure and values of the input time series like in the previous naive algorithm. Therefore its best-case, worst-case, and average-case time complexity is always O(n2). 20 Chapter 3 Algorithms and Complexity 3.2 Space complexity This algorithm only requires one extra auxiliary variable (slopemax), of constant size in memory O(1). Therefore, as with the previous naive algorithm, the memory required is mainly dictated by the size of the input time series and the size of the output graph, again, in the order of O(n2). 21 Chapter 3 Algorithms and Complexity 3.3 3.3 Divide-and-conquer strategy Algorithm 3 Visibility Graph - Divide-and-Conquer Strategy 1: function VisibilityGraphT=(t1, y1),...,(tn, yn) 2: V←{1, . . . , n} 3: E←VisibilityGraphRecursive(T,1,n) 4: return G= (V, E) 5: end function 6: 7: function VisibilityGraphRecursiveT,left,right 8: E0←∅ 9: if left <right then 10: imax ←index i,left ≤i≤right, with highest yi 11: for a←left to right except imax do 12: if node imax ‘sees’ node athen .i.e. using the slope strategy locally 13: add edge (imax, a)to E0 14: E0←E0∪VisibilityGraphRecursive(T,left,imax −1) 15: E0←E0∪VisibilityGraphRecursive(T,imax + 1,right) 16: return E0 17: end function When looking for an even more optimal strategy we can see that we can benefit from the fact that very large data points in the time series act as a kind of barrier or wall between the data points found to their left and those found to their right [33], preventing most nodes from the two different sides from being connected in the visibility graph. We can make use of this property to prune many of the candidate edges, substantially reducing the number of potential edges that need to be checked. In particular, the highest data point in the time series will prevent any data point to its left from having line-of-sight visibility with any node to the right of this highest point. We can, therefore, divide the input time series into two segments split by this highest data point, and consider the two resulting smaller time series as two independent visibility graph subproblems, without having to check for edges corresponding to nodes that belong to different sides, as we already know that none will exist. The visibility graph of the original time series will be the result of taking these two smaller visibility graphs together and adding the node of the highest data point and its edges (which, in this case, can indeed reach both sides, and will need to be checked using, ideally, the slope strategy described in Section 3.2). This idea can be recursively applied to obtain the visibility graph of each of the two mentioned smaller time series, whose highest data point can again split the segment into two even smaller subproblems. This can continue until the time series segments cannot be further divided. See Figure 3.3 for an example. Correctness This algorithm builds on top of the previously described algorithms and greatly reduces the number of potential edges that need to be checked. Let (tmax, ymax)be the observation with the highest value of the input time series T(such that ymax ≥yifor all (ti, yi)∈ T). Then we know for certain that there cannot be an edge connecting two nodes a,bif they are on different sides of tmax (e.g. ta< tmax < tb). As, by definition, all points other than (tmax, ymax) have a value not higher than ymax and thus will not have line-of-sight visibility with any point at the other side of (tmax, ymax)because (tmax, ymax)will always be an intermediate data point causing an obstruction. 22 Chapter 3 Algorithms and Complexity 3.3 1 2 3 4 5 6 7 8 9 10 11 . & 1 2 3 4 5 7 8 9 10 11 . & & 1 2 4 5 9 10 11 Figure 3.3: Divide-and-conquer strategy demonstrated in an example time series. In this case, the maximum value of the time series corresponds to y6, so we connect its corresponding node with all other nodes in sight, and then we recursively repeat the same strategy on both smaller time series formed by the points to the left and to the right of t6. Let TLbe the smaller time series formed by all data points to the left of (tmax, ymax), and equivalently TRbe the smaller time series formed by all data points to the right (note that (tmax, ymax)is not included in either). No node from TLis connected to a node in TRin the visibility graph of T, and therefore the visibility graph of Tcan be constructed by the union of the nodes and edges of the visibility graphs of TLand TR, with the addition of the node corresponding to (tmax, ymax)and its edges. Both TLand TRcan be considered independent visibility graph subproblems and we can therefore recursively apply the same algorithmic idea again to obtain their visibility graphs. The recursion stops when the time series segments become too small and cannot be further divided. Time complexity We explored the different behaviors experienced in the worst and bestcase scenarios for this divide-and-conquer strategy. Due to the nature of the algorithm, the time required for its execution can vary significantly even when comparing different time series of the same input size n. The running time of the algorithm will depend on the placement of the maximum values of the time series at each recursion step and how balanced each split happens to be. In any case, the running time of each subproblem of input size n0is O(n0). This is, the function in line 7 in Algorithm 3, excluding the two subsequent recursive calls (excluding lines 14 and 15), when given an input time series of size n0has cost O(n0). This arises from the linear cost of finding the maximum value of the time series (in line 10), and the also linear cost of finding all the appropriate edges attached to the node corresponding to this maximum value (in line 11), which can be achieved in linear cost by taking advantage of the slope strategy as described in Section 3.2. 23 Chapter 4 Technologies and Implementation 4.5 The package benefits from an efficient C backend for its functions (with the use of Cython) while still providing native integration in the Python environment. Therefore, ts2vg can easily work with input data from many sources and providers as available with existing Python tools, as well as allowing the study and analysis of the resulting visibility graphs with any of the many existing graph analysis, data science, and visualization packages and tools available for Python. A minimum working example, useful for many use cases, of ts2vg as used in a Python script can be: from ts2vg import NaturalVisibilityGraph ts = [0.87 ,0.48,0.36,0.83,0.87,0.48,0.36,0.83] edges =NaturalVisibilityGraph(ts).edgelist() Similarly, an igraph object representation of the graph [12] (see Section 4.5.1) can be very easily obtained with: vg =NaturalVisibilityGraph(ts).as_igraph () The divide-and-conquer strategy is used in all cases by default. ts2vg can also be used as a stand-alone command line program useful to quickly obtain visibility graphs from time series without having to type any Python code. As a basic example, the following command can be used from the console to obtain the edge list (as saved to a file out.edg) of the visibility graph resulting from an input time series found in a local timeseries.txt file: >>> ts2vg ./ timeseries.txt −o out.edg Full source code, documentation, and installation instructions for the ts2vg Python package and program can be found at: https://github.com/CarlosBergillos/ts2vg 4.4.1 Potential future considerations The current state of our implementation and the associated ts2vg package, although stable, and successfully fulfilling the initial requirements, need not be final, it allows for future additions and improvements, some of which are out of the scope of the project or proved not viable due to technical or time limitations. It could be expanded and improved, for example, in the following ways: •Adding new functionalities by providing implementations for other visibility graph algorithms (e.g. horizontal visibility graphs) and/or enrichening and expanding the existing ones (e.g. directed and weighted versions of visibility graphs). •Exploring advanced performance improvements to the existing algorithms, for example, with the use of parallel and multithreaded strategies or with the use of GPU accelerated implementations, exploring tools like CUDA [44]. The above are just presented as ideas worthy of consideration, but in no case are planned future developments and might even prove to not be practical. 30 Chapter 4 Technologies and Implementation 4.5 4.5 Auxiliary functionalities To conduct experiments and explore potential visibility graph applications (see Chapter 5) we also needed the help of many auxiliary functionalities, whose implementation was outside of the scope of this project, and in most cases pointless due to different already existing and well-established Python libraries suited for the tasks. 4.5.1 Graph analysis functionalities In order to conduct different graph analysis tasks, we used the Python igraph package [12]. This package provides many easy-to-use graph construction, manipulation, and analysis functionalities, such as, to name a few, obtaining their degree distribution, obtaining path lengths and many other measures, applying different clustering algorithms and plotting the graphs using different visualization strategies. The Python igraph package is (for the most part) only an interface for a core implementation of igraph written in C, again, benefiting from more efficient and portable methods. Full official documentation for igraph has been used and can be found in [50]. 4.5.2 Generating Wiener process and fractional Brownian motion series A major part of our experimentation with visibility graphs required the availability of stochastic time series as input for the visibility graph algorithms. Understanding the behavior of Wiener process series and fractional Brownian motion series can help in forming conclusions that might then be useful for real-world time series. We used the Python stochastic package [17] to construct these series. The implementations of fractional Brownian motion series provided by this package use the methods described by Hosking 1984 [22] and by Davies and Harte 1987 [13]. 4.5.3 Estimating power-law exponents We used our own Python implementation (with the support of numpy functions) of the formula for estimating the exponent γof power-law data using MLE (equation (2.16) shown in Section 2.6 and based on [41]). For more complex distribution estimations we also studied and considered using the Python powerlaw package [1]. 31 5Applications and Experiments As part of our research on visibility graphs we will empirically test some of the ideas described in previous chapters, we will further study some of the properties of visibility graphs and we will explore some of their potential applications. We will be using our visibility graph algorithms implementations (using the Python programming language and Cython, as described in Chapter 4). Detailed information on relevant software packages and hardware used for all the experiments conducted and described in this chapter can be found in Appendix A. 5.1 Performance evaluation of the different algorithms In Chapter 3 we described and mathematically evaluated the time and space complexity of three different algorithm strategies that can be used to compute visibility graphs from an input time series. To verify such theoretical analysis, we decided to test and measure the time performance of our implementations of the algorithms. In Table 5.1 we show the execution time of the different visibility graph algorithms used on Gaussian white noise time series of different length. As expected, there is a clear difference in the running time of the different algorithmic strategies, the divide-and-conquer strategy being the fastest one, and, on the other hand, the naive algorithm the slowest one. It is important to consider that the visibility graph algorithms can require drastically different execution times for different input time series, even when of the same length n. For example, the execution time of the divide-and-conquer algorithm depends both on the time series length n, and on the distribution of maximum values through the time series (how balanced the time series is). Gaussian white noise time series can be considered a very well balanced time series, and thus greatly benefit from the divide-and-conquer algorithm. Table 5.2 shows the execution time of the same algorithms but this time with Wiener process time series as input, a less balanced group of time series, but maybe closer to more realistic use cases of visibility graphs. In this case, as expected the slope algorithm shows no significant time difference when compared to the previous white noise case, but both the naive and the divide- and-conquer algorithms require noticeable longer execution times. From these results and the theoretical analysis from previous sections, we can see and conclude that there is no real reason to prefer the naive or the slope algorithms over the divide-and-conquer algorithm for any practical application. The divide-and-conquer algorithm shows significantly better time performance at the expense of a very small additional memory penalty. 32 Chapter 5 Applications and Experiments 5.2 Execution time on Gaussian white noise time series n(×104) Naive Algorithm (s) Slope Algorithm (s) D&C Algorithm (s) 1 0.6141 0.0964 0.0066 2 2.6288 0.3736 0.0140 3 6.1971 0.8321 0.0213 4 11.3760 1.4703 0.0289 5 18.3866 2.2898 0.0366 6 26.6579 3.2924 0.0439 7 36.8757 4.4717 0.0511 8 48.6662 5.8315 0.0588 9 62.4979 7.3732 0.0668 10 78.4211 9.0952 0.0750 100 - - 0.7977 Table 5.1: Execution times of the different algorithmic strategies (naive, slope, and divide-and-conquer) when used to obtain the visibility graph of Gaussian white noise time series of different lengths n(averaged over 10 samples each). Execution time on Wiener process time series n(×104) Naive Algorithm (s) Slope Algorithm (s) D&C Algorithm (s) 1 4.4416 0.1099 0.0246 2 23.6012 0.3977 0.0565 3 66.0205 0.8730 0.0866 4 134.4352 1.5266 0.1242 5 227.7911 2.3682 0.1526 6 417.1478 3.3837 0.2018 7 653.0515 4.5799 0.2308 8 703.4591 5.9542 0.2841 9 907.9010 7.5241 0.3352 10 - 9.2722 0.4089 100 - - 8.4202 Table 5.2: Execution times of the different algorithmic strategies (naive, slope, and divide-and-conquer) when used to obtain the visibility graph of Wiener process time series of different lengths n(averaged over 10 samples each). 5.2 Periodic series result in regular visibility graphs Periodic time series are those formed by the repetition of one same subsequence over and over: T= (y1, y2, . . . , yp, y1, y2, . . . , yp, y1, y2, . . . , yp, . . . ) The length pof this subsequence that gets repeated is called the period of the series. Two examples of periodic time series are shown in Figure 5.1. The visibility graph resulting from a periodic time series will be formed by the many repetitions of one same graph motif, as each instance of the repeated subsequence results in the same equivalent subgraph. Therefore, the resulting visibility graphs present a very regular degree distribution, as there will be at most punique degree values for all nnodes in the graph (excluding the ones from the nodes in the two endpoint segments, but the effect of those on the degree distribution is negligible 33 Chapter 5 Applications and Experiments 5.2 0 100 200 300 400 500 t −0.5 0.0 0.5 1.0 y A 0 100 200 300 400 500 t 0 1 2 y B Figure 5.1: First 500 observations of two example periodic time series, both of repeating period 32. when nis large, with pn). It can also be easily seen that the degree distribution will stay essentially constant independently of the length nof the periodic time series (for a fixed, repeated subsequence of length p, with pn). Additionally, these degrees will be quite small as the highest data point in the sequence appears repeatedly every pobservations, effectively segmenting the whole time series into repeated independent regions. In particular, each of the nodes corresponding to the highest data point in the segments will have degree at most 2p, and all other non-maximum nodes will have degree at most p. Clearly, there are no hub nodes, no power-law distribution, and the diameter and average shortest path length are very large (on the order of n/p). In Figure 5.2 we show the degree distributions of the visibility graphs obtained from the two time series shown in Figure 5.1 (when continued up to length 105). Arguably, the discrete degree distributions resulting from periodic time series can be considered as a sort of fingerprint of the time series periods [30]. 34 Chapter 5 Applications and Experiments 5.3 0 10 20 30 degree k 0.00 0.05 0.10 0.15 0.20 P(k) A 0 10 20 30 degree k 0.00 0.05 0.10 0.15 0.20 P(k) B Figure 5.2: Degree distributions of the two periodic time series, Aand B, shown in Figure 5.1, continued up to length 105. Note that for both time series no node has degree larger than 32 (the period). The distributions show clearly distinct discretized structures that could be used as a sort of fingerprint of the periodic time series. This distribution will stay essentially constant with other varying lengths nof the series, as long as the repeated period subsequence stays the same. 5.3 Fractal series result in scale-free visibility graphs An interesting finding in the field is that visibility graphs resulting from fractal series (like Wiener process series, fractional Brownian motion series, and other mathematical and real-world series of similar fractal behavior) present small-world and real-world network properties [30]. We generated different samples of Wiener process series (such as the one shown in Figure 5.3), obtained their visibility graph, and confirmed that they match the properties previously presented in Section 2.3, usually associated to real-world networks. 0 200 400 600 800 1000 t −20 −10 0 10 20 y Figure 5.3: First 1000 observations of an example Wiener process time series (also known as Brownian motion). These properties are: Small diameter These visibility graphs present diameter lengths and average shortest path lengths that are orders of magnitude smaller than the number of nodes in the graphs. In Figure 5.4 we show empirical measurements for these values as a function of the time series length n. Both the diameter and the average shortest path lengths seem to grow logarithmically with respect to the graph size n. We observe very small path lengths joining any two nodes, i.e. paths consisting of less than 10 edges even in graphs with millions of nodes. This is possible, in part, thanks to the presence of some nodes with a very high degree, commonly known as hub nodes that allow quick navigation through the network. High clustering coefficient In the case of the Wiener time series shown in Figure 5.3, it presents a remarkably high average clustering coefficient (also known as transitivity), 35 Chapter 5 Applications and Experiments 5.4 103104105 n 4 6 8 10 12 14 16 18 Path length (nº of edges) Diameter Average shortest path Figure 5.4: Medians of the measured diameters of the visibility graphs obtained from 15 sampled Wiener process time series (blue squares) and their average shortest path length (black circles) as a function of the network size (and time series length) n. The two measures seem to grow logarithmically with respect to the network size for these time series (note that the xaxis is using alog scale). A best fitted logarithmic function line has also been obtained and is provided for each case. with values of 0.685 ±0.01 for many different lengths nconsidered for the series. This finding also holds for other Wiener process series studied. Power-law degree distribution As shown in Figure 5.5 the Wiener process time series presents a degree distribution clearly compatible with a power-law distribution (starting after a certain xmin value, and getting progressively noisy at the end of its tail). From this, it is clear that most of the nodes in the graph have a very small degree, and only a small amount of nodes have a large degree. These nodes of large degree (hub nodes) might correspond to regional peaks in the input time series (data points with higher yvalue than their surroundings, and therefore having direct visibility with many other data points in the series). 0 200 400 600 800 degree k 0.00 0.02 0.04 0.06 0.08 0.10 P(k) 101102103 degree k 10−5 10−4 10−3 10−2 10−1 P(k) Figure 5.5: Linear plot (left) and log-log plot (right) of the degree distribution of the visibility graph obtained from the Wiener time series shown in Figure 5.3 and as continued to length 105. The tail of the degree distribution visually follows a distinct power-law shape. 5.4 Estimating the Hurst exponent via the visibility graph The Hurst exponent of time series is an important and useful measure used in time series analysis. It provides information on the volatility of the series, and can have practical applications, for example, on financial time series analysis, and it could aid in trying to predict the future behavior of the series. Visibility graphs obtained from time series have been shown to provide a new unique 36 Chapter 5 Applications and Experiments 5.4 way of estimating the Hurst exponent of such time series. In particular, the exponent of the power law that approximates the degree distribution of the visibility graph is used for this estimation. More so, a linear dependency between this power-law exponent and the Hurst exponent of the time series has been proven to exist [31, 42]. Fractional Brownian motion series (characterized by the Hurst exponent) can be used to model certain financial series and real-world time series. In Figure 5.6 we show the degree distributions of three visibility graphs obtained from fractional Brownian motion series characterized by different values of H. As it can be seen, the three visibility graphs exhibit scale-free behavior and can be characterized as power laws P(k) = ck−γof different shape. For each of the distributions, we have estimated its best-fitting exponent γusing Maximum Likelihood Estimation (MLE) as explained in Section 2.6 and using xmin = 8. And at first glance, we can observe that the values for γdecrease as Hincreases. 101102103 k 10−5 10−4 10−3 10−2 10−1 P(k) γ= 2.70 H= 0.2 101102103 k γ= 1.98 H= 0.5 101102103 k γ= 1.39 H= 0.8 Figure 5.6: Degree distributions (in log-log axes) of three visibility graphs each obtained from a different fBm series of length 105with H= 0.2(left), H= 0.5(center) and H= 0.8(right). The tails of the three degree distributions follow a power law P(k) = ck−γ, each one with a different shape and with noticeably different exponents γ. Lacasa et al. 2009 [31] showed that the exponent γof visibility graphs is linearly related to the Hurst exponent Hof the series in the form: γ(H)=3−2H(5.1) We conducted a concise experiment trying to replicate these findings. In Figure 5.7 we show the result of empirical tests of this relationship between γand Hand we compare it with the previous theoretical formula (5.1). The empirical tests seem to follow the theoretical formula quite accurately, but it can be observed that the values of γstart to deviate for values of Hclose to 0 and close to 1, and being the most accurate near H=1/2. These small deviations could be due to biases in the fBm series generated, limitations of working with series of finite length, and inaccuracies in the values of γobtained through MLE. 37 Chapter 5 Applications and Experiments 5.5 0 0.2 0.4 0.6 0.8 1 H 1.0 1.5 2.0 2.5 3.0 γ γ(H) = 3 −2H Figure 5.7: Numerical estimation (using MLE, xmin = 8) of the exponent γof the visibility graphs obtained from fBm series of length 104 with different values of H. For each value of H, 20 fBm samples have been used. The dashed line corresponds to the theoretical formula γ(H) = 3−2H. 5.5 Segmentation of time series via graph clustering Automatic segmentation of time series can also be a useful application of visibility graphs. Time series segmentation consists on dividing an input time series into a sequence of discrete and disjoint segments according to some notable underlying property. For example, we might be interested in segmenting an almost periodic time series (i.e. a periodic time series that might have been subjected to noise on the yand/or tdimension). This can be achieved on certain time series by applying clustering algorithms to their resulting visibility graphs, and extrapolating the discovered clusters (or communities) to the original time series data points. Different clustering algorithms can be used to solve the task, achieving potentially different results. A clustering algorithm based on modularity optimization presented by Blondel et al. 2008 [7] has been used in all our experiments presented in this section. In Figure 5.8-A we show the result of segmenting an almost-periodic time series using this technique. We also show the result of two other time series examples (B and C) whose periodicity is much more diffuse, but the technique still manages to divide the time series into segments of distinct behavior. Note that due to the nature of visibility graphs, clusters of nodes are likely to correspond to valleys of the time series (where direct visibility among data points is easier), and regional peaks of the series are likely to divide segments. Therefore, the resulting segmentation of the time series will also reflect this behavior. It is easy to see that peak-centered segments (separated by valleys) could also be obtained by multiplying the time series heights by −1before obtaining its visibility graph. It is worth noting that the clustering algorithms used have no notion of time in their graph analysis, but still sequential-in-time clusters naturally emerge in many practical cases. More complex techniques could be used that include a weighting of the edges based on the difference in time between two connected nodes, therefore aiming to avoid broken-in-time clusters as it happened in Figure 5.8-C. Additionally, it can be theorized that different possible segmentations of one same time series (corresponding to the potentially overlapped different frequencies of the time series) could be obtained by optimizing for clusters of different sizes in each case. 38 Chapter 5 Applications and Experiments 5.5 −20 0 20 y A 0 100 200 300 400 500 t −20 0 20 y −20 0 y B 0 100 200 300 400 500 t −20 0 y −10 0 10 y C 0 100 200 300 400 500 t −10 0 10 y Figure 5.8: Results of time series segmentation via visibility graph clustering used in three different input time series. Each color corresponds to a distinct cluster/segment obtained in the time series. An example almost-periodic time series is used in A, and two chosen Wiener-process-like series in Band C. 39 Appendices 46 AConditions and Specifications To conduct all the performance analysis, visibility graph tests, and experiments presented throughout this project we have used one same computer with the characteristics and specifications shown in Table A.1. Component Description CPU Intel®Core™i5-4690K, 4 Cores @ 3.50GHz RAM 16 GB Memory Crucial®MX500 1000GB SSD OS Manjaro Linux 19.0.2 Kyria (rolling release 03-2020) Table A.1: Details and specifications of the computer used to conduct the experiments presented in this project. Details on the release versions used for Python and for any other relevant software tools are shown in Table A.2. Component Version Python 3.8.2 Cython 0.29.17 gcc 9.3.0 numpy 1.18.2 igraph 0.8.0 stochastic 0.4.0 Table A.2: Details on the Python release version and on the associated relevant software tools and libraries used for the implementations and experiments presented in this project. Additionally, all diagrams and figures shown throughout this whole project have been created by the author with the help of Python and Matplotlib [23]. 47 BAlgorithm’s Source Code Source code for working implementations of the visibility graph algorithms that have been described and studied in Chapter 3 are presented here. These implementations, as provided here, have been used for all experiments conducted and presented in Chapter 5. The programming language used is Cython, a superset of the Python programming language with additional support for static typing and easy integration of C methods. Additionally, we later provide information on the source code for the ts2vg package. 48 Appendix B Algorithm’s Source Code B.1 B.1 Naive strategy implementation def visibility_graph_naive(np.float64_t[:] ts): """ Computes the visibility graph of a time series using a naive strategy. Args: ts (numpy 1d array): Time series data. Returns: g (igraph.Graph): The visibility graph of ‘ts ‘. """ cdef uint n=ts.size g=Graph () g.add_vertices(n) edges = [] cdef uint a,b for ain range(n−1): for bin range(a+1,n): if _no_obstructions(ts,a,b): edges.append((a,b)) g.add_edges(edges) return g cdef bint _no_obstructions(np.float64_t [:] ts,uint a,uint b): """ True if there is direct line−of−sight visibility (no intermediate obstructions) between data points at positions ‘a‘ and ‘b‘ in the time series ‘ts‘. False otherwise. """ cdef float y_a =ts[a] cdef float y_b =ts[b] cdef float y_c cdef uint c for cin range(a+1,b): y_c =ts[c] if y_c >= y_b + (y_a−y_b)∗(b−c)/(b−a): return False return True 49 Appendix B Algorithm’s Source Code B.2 B.2 Slope strategy implementation def visibility_graph_slope(np.float64_t[:] ts): """ Computes the visibility graph of a time series using a slope strategy. Args: ts (numpy 1d array): Time series data. Returns: g (igraph.Graph): The visibility graph of ‘ts ‘. """ cdef uint n=ts.size g=Graph () g.add_vertices(n) edges = [] cdef uint a,d cdef float y_a,y_b cdef float max_slope,slope for ain range(n−1): y_a =ts[a] max_slope =−INFINITY for din range(1,n−a): y_b =ts[a+d] slope = (y_b−y_a)/d# d = b−a if slope >max_slope: edges.append((a,a+d)) max_slope =slope g.add_edges(edges) return g 50 Appendix B Algorithm’s Source Code B.3 B.3 Divide-and-conquer strategy implementation def visibility_graph_dc(np.float64_t [:] ts): """ Computes the visibility graph of a time series using a divide−and−conquer strategy. Args: ts (numpy 1d array): Time series data. Returns: g (igraph.Graph): The visibility graph of ‘ts ‘. """ cdef uint n=ts.size g=Graph () g.add_vertices(n) edges = [] cdef uint left,right ,i,d cdef float y_i,y_a cdef float max_slope,slope cdef PairQueue queue =PairQueue() queue.push((0,n)) while not queue.is_empty(): (left,right) = queue.pop() if left+1<right: i=_argmax(ts,left,right) y_i =ts[i] max_slope =−INFINITY for din range(1,i−left+1): y_a =ts[i−d] slope = (y_a−y_i)/d# d = a−i if slope >max_slope: edges.append((i,i−d)) max_slope =slope max_slope =−INFINITY for din range(1,right−i): y_a =ts[i+d] slope = (y_a−y_i)/d# d = a−i if slope >max_slope: edges.append((i,i+d)) max_slope =slope queue.push((left,i)) queue.push((i+1,right)) g.add_edges(edges) return g 51 Appendix B Algorithm’s Source Code B.4 B.4 Streaming algorithm implementation class VisibilityGraphStream: """ Class representation for a visibility graph updated via a stream of data. The graph obejct can be accessed at any point via the variable ‘vg ‘. """ def __init__(self): """ Method executed on initialization. """ self.vg =Graph () self.ts =np.empty(0) self.n=0 self.__max_i =−1 self.__max_v =−INFINITY def update(self,new_y): """ Updates the visibility graph with one new time series observation. A new node is added to the graph , as well as all required new edges. """ self.vg.add_vertex() self.ts =np.append(self.ts,new_y) edges = [] max_slope =INFINITY limit =self.__max_i if new_y >= self.__max_v: if new_y >self.__max_v: limit =0 self.__max_i =self.n self.__max_v =new_y for ain reversed(range(limit ,self.n)): slope = (self.ts[a]−self.ts[self.n])/(a−self.n) if slope <max_slope: edges.append((self.n,a)) max_slope =slope self.vg.add_edges(edges) self.n+= 1 52 Appendix B Algorithm’s Source Code B.5 B.5 ts2vg package A selection of the previously shown functions and algorithms has been included and adapted in the ts2vg Python package that we have developed. ts2vg aims to be an easy-to-use and fast tool for anyone to obtain and work with visibility graphs both in Python and directly through a command-line interface. The full source code for the ts2vg package is public and can be found at: https://github.com/CarlosBergillos/ts2vg 53 CProject Management C.1 Tasks A summary of the tasks required in order to successfully complete this project is provided in Table C.1, a more detailed description of tasks and their subtasks is provided later. Task Description Est. Hours T.0 Project Management (GEP) Project planning, management and assessment as required for the GEP course. 90 T.1 Literature Review & Research Research and study on visibility graphs topics and all the related time series, complex networks, and data science ideas. 100 T.2 Algorithms Study Studying different visibility graphs algorithms, evaluating their correctness and computational complexity. 80 T.3 Algorithms Implementation Implementing and testing the visibility graph algorithms and related algorithms and functions. 100 T.4 Applications & Experiments Conducting case studies and exploring real-world applications for the visibility graph ideas. 100 T.5 Documentation Writing the final report of the project and any related documentation, including the creation of any required tables and figures. 60 T.6 Defense Preparation Preparing the required slides and material for the oral presentation, rehearsing it and preparing for potential questions. 10 Total: 540 Table C.1: Summary of the different tasks planned for the project and their estimated durations in hours. 54 Appendix C Project Management C.1 T.0 Project Management (GEP) [90 hours] Project planning, management, and assessment as required for the GEP course. T.0.1 Deliverable 1 Work on the first deliverable for the GEP course, this deliverable will focus on the context and scope of the project. T.0.2 Deliverable 2 Work on the second deliverable for the GEP course (this document), for this deliverable we will define the project tasks and resources and establish a viable time planning for the project. T.0.3 Deliverable 3 Work on the third deliverable for the GEP course, this deliverable will focus on the budget and sustainability aspects of the project. T.0.4 Deliverable 4 Work on the fourth deliverable for the GEP course, this deliverable will combine the work and improve on the three previous deliverables. T.1 Literature Review & Research [100 hours] The project requires deep research and study on many concepts and mathematical ideas. Before starting to write any code, we need to conduct research on visibility graphs topics and all the related time series, complex networks, and data science ideas. T.1.1 Visibility Graphs Theory Exploring literature for visibility graphs concepts. T.1.2 Time Series Theory & Models Exploring literature for time series ideas and models such as the Brownian Motion Series. T.1.3 Complex Networks Theory Exploring literature for complex networks theory ideas, their measures, and the properties of real-world networks. T.1.4 Data Science Theory Exploring literature and researching on different data science concepts such as regression models. T.2 Algorithms Study [80 hours] Studying different visibility graphs algorithms, evaluating their correctness and computational complexity. T.2.1 State-of-the-Art Algorithms Study and conceptualize different state-of-the-art visibility graph algorithms. T.2.2 Complexity Study Analyze and compare the computational complexity of the previously studied algorithms. T.3 Algorithms Implementation [100 hours] Implementing and testing the visibility graph algorithms and related algorithms and functions. T.3.1 Research on Tools & Libraries Properly choosing which programming language, tools, libraries and technologies to use for the implementation and programming development is a very important decision, as this cannot be easily changed further down the line once the development has started. For this reason, we are devoting substantial time to research and compare different options available to us. 55 Appendix C Project Management 6.8 plicable field, in particular in the fields related to computation, perception and operation in intelligent environments. In this project we have worked with time series as a formalization and representation of real-world measurements, and, on the other hand, we have also worked with graphs as a formalization and representation of the same data. Visibility graphs as presented provide a unique bridge between the two (traditionally independent) methods. •CCO2.4: To demonstrate knowledge and develop techniques about computational learning; to design and implement applications and system that use them, including these ones dedicated to the automatic extraction of information and knowledge from large data volumes. This project required a complex and diverse theoretical mathematical background regarding techniques suitable to extract information from potentially very large data sources (e.g. using data science ideas, time series analysis and complex network theory) as presented in Chapter 2 (Definitions and Theoretical Background) and later put into practice in Chapter 5 (Applications and Experiments). •CCO3.1: To implement critical code following criteria like execution time, efficiency and security. One of the initial goals of the project was to implement efficient solutions to compute visibility graphs even for very large input sizes. Execution time and efficiency of our code has therefore been an important requirement in our implementations, and as described in Chapter 4 (Technologies and Implementation), special attention has been payed to choose adequate tools and libraries and to develop optimized implementations of the algorithmic strategies chosen. 62