Visibility graphs of random scalar fields and spatial data
Lucas Lacasa Jacopo Iacovacci








School of Mathematical Sciences, Queen Mary University of London, Mile End Road, London E14NS (UK)




             


The family of visibility algorithms were recently introduced (Lacasa et al, PNAS 105 (2008)) as mappings between time series and graphs. Here we extend this method to characterize spatially extended data structures by mapping scalar fields of arbitrary dimension into graphs. After introducing several possible extensions, we provide analytical results on some topological properties of these graphs
associated to some types of real-valued matrices, which can be understood as the high and low disorder limits of real-valued scalar fields. In particular, we find a closed expression for the degree distribution of these graphs associated to uncorrelated random fields of generic dimension, extending a well known result in one-dimensional time series. As this result holds independently of the field's marginal distribution, we show that it directly yields a statistical randomness test, applicable in any dimension. We showcase its usefulness by discriminating spatial snapshots of two-dimensional white noise from snapshots of a two-dimensional lattice of diffusively coupled chaotic maps, a system that generates high dimensional spatio-temporal chaos.
We finally discuss the range of potential applications of this combinatorial framework, which include image processing in engineering, the description of surface growth in material science, soft matter or medicine and the characterization of potential energy surfaces in chemistry, disordered systems and high energy physics. An illustration on the applicability of this method for the classification of the different stages involved in carcinogenesis is briefly discussed.



                             

Visibility graphs of random scalar fields and spatial data
    Lucas Lacasa and Jacopo Iacovacci
    December 30, 2023
==========================================================



 




§ INTRODUCTION

Visibility and Horizontal Visibility Graphs are a family of mappings between ordered sequences and graphs <cit.>. Consider an ordered sequence { x(t)}_t=1^N, where x(t)∈ℝ^m,  m≥ 1. For m=1 a typical case of such sequence are time series describing the activity of some system, whereas for m>1 we consider multivariate time series or in general high-dimensional dynamical systems. In every case, a univariate time series of N data is mapped into a graph of N nodes such that two nodes are linked in the graph if a particular visibility criterion holds in the sequence (the multivariate setting has been explored recently <cit.>). This mapping enables the possibility of performing graph-theoretical time series analysis and builds a bridge between the theories of dynamical systems, signal processing and graph theory.


In recent years, this mapping has been used to provide a topological characterization of different routes to low dimensional chaos <cit.>, or different types of stochastic and chaotic dynamics <cit.>. From an applied angle, it is being widely used to extract in a simple and computationally efficient way informative features for the description and classification of empirical time series appearing in several areas of physics including optics <cit.>, fluid dynamics <cit.>, geophysics <cit.> or astrophysics <cit.>, and extend beyond physics in areas such as physiology <cit.>, neuroscience <cit.> or finance <cit.> to cite only a few examples. Whenever each element in a given classification task is naturally encoded as an ordered sequence,
one can map such sequence into a visibility graph and subsequently extract a certain set of topological properties of these graphs as the feature vector with which to train classifiers in supervised learning tasks.


Here we propose to extend this methodology from time series { x(t)}_t=1^N to  scalar fields h(x,y):ℝ^d→ℝ. This extension, which has only been scarcely explored <cit.> is conceptually closer to the original context of visibility graphs <cit.> and enables the possibility of constructing the visibility graphs of images, landscapes, and general large-scale spatially-extended surfaces. In what follows we will introduce the concept along with a few definitions and properties. In section III we provide analytical results on some topological properties of these graphs
associated to some types of real-valued matrices which can be understood as the high and low disorder limits of real-valued scalar fields. In particular, we find a closed expression for the degree distribution of these graphs associated to uncorrelated random fields of generic dimension, extending the result known for one-dimensional time series. As this result holds independently of the field's marginal distribution, we show that this result directly yields a statistical randomness test, applicable in arbitrary dimensions. In section IV we showcase its usefulness by discriminating two-dimensional white noise from two-dimensional lattice of diffusively coupled chaotic maps (a system that generated high dimensional spatio-temporal chaos). In section V we discuss the range of potential applications of this combinatorial framework and we further briefly illustrate its usefulness for characterizing the process of oncogenesis.




§ DEFINITIONS AND BASIC PROPERTIES

Definition (VG) Let S={x_1,…,x_N} be a ordered sequence of N real-valued, scalar datapoints. A Visibility Graph (VG) is an undirected graph of n nodes, where each node i∈ [1,N] is labelled by the time order of its corresponding datum x_i. Hence x_1 is mapped into node i=1, x_2 into node i=2, and so on.
Then, two nodes i and j (assume i<j without loss of generality) are connected by a link if and only if one can draw a straight line connecting x_i and x_j that does not intersect any intermediate datum x_k,  i<k<j. Equivalently, i and j are connected if the following convexity criterion is fulfilled:

    x_k< x_i + k-i/j-i[x_j-x_i], ∀ k: i<k<j

The same definition applies to a Horizontal Visibility Graph (HVG) but in this latter graph two nodes i, j (assume i<j without loss of generality) are connected by a link if and only if one can draw a horizontal line connecting x_i and x_j that does not intersect any intermediate datum x_k,  i<k<j. Equivalently, i and j are connected if the following ordering criterion is fulfilled:


    x_k<inf(x_i,x_j), ∀ k: i<k<j

From a combinatoric point of view, HVGs are outerplanar graphs with a Hamiltonian path <cit.>, i.e. noncrossing graphs as defined in algebraic combinatorics <cit.>.
Note that the former definitions focus on discrete sequences, such that the index labeling is such that i+1≡ i+Δ, where Δ is the spacing between data. Interestingly, both VG and HVG are invariant under changes in Δ. In particular, this enables to directly consider the continuous version of a discrete time series simply as the limit Δ→ 0. This invariance property will allow treating continous scalar fields as the Δ→ 0 limit of matrices as we will show later.
  

Extension classes. One can now extend the definition of visibility to handle two-dimensional manifolds, by simply extending the visibility criteria along one-dimensional sections of the manifold. The question is, in how many different ways one can do that? As a matter of fact, there exist several possibilities, here we consider just a few of them. We firstly consider manifolds of dimension d which have a natural Euclidean embedding and define two extension classes, labelled as canonical and FCC respectively. In the canonical extension class, the rule of thumb for extending the definition of a visibility graph to a manifold of dimension d will be by applying the VG/HVG to d orthogonal sections of the manifold (which define n=2d directions). In other words, at each point of the manifold one constructs the VG/HVG in the direction of the (canonical) Cartesian axis. On the other hand, the FCC extension class allows an additional number of sections in the direction of the main diagonals. Accordingly, in this second class the number of directions is n=2d+2^d directions (see figure <ref> for an illustration in the case d=2). Finally, a third extension class (which in this work will only be studied for d=2 flat surfaces) is defined by taking n directions in such a way that the set of n vectors make an homogeneous angular partition of the plane with constant angle 2π/n. This class is labelled as the order-n class. Obviously, the order-8 and order-4 classes coincide, when d=2, with the FCC and canonical classes respectively. These special classes are indeed of special relevance as they are the most natural algorithmic implementation for image processing <cit.>. We are now ready to give a more formal definition of visibility graphs in these extension classes.






Definition (IVG_n) Let I be a N× N matrix, where I_ij∈ℝ and N>0. For an arbitrary entry ij, make an angular partition of the plane into n directions, such that direction labelled as p makes an angle with the row axis of 2π(p-1)/n. The Image Visibility Graph of order n IVG_n is a graph with N^2 nodes, where each node is labelled by a duple ij in association with the indices of the entry I_ij, such that two nodes ij and i'j' are linked if


  * i'j' belongs to one of the n angular partition lines, and

  * I_ij and I_i'j' are linked in the VG defined over the ordered sequence which includes ij and i'j'.

The Image Horizontal Visibility Graph (IHVG_n) follows equivalently if in the second condition we make use of HVG instead of VG. Note that in the preceding definition, I can be understood as a two-dimensional square lattice, which is naturally embedded in ℝ^2 if we associate a certain lattice length Δ_p>0 to the separation between any two neighbors in each direction p. In the limit N→∞, Δ_p→ 0 this matrix I converges in some mathematically well-defined sense to a continuous scalar field h(x,y):ℝ^2→ R. Accordingly, the continuous version of these graphs can be obtained for n→∞, and in that case I(H)VG_∞ would be an infinite graph. In this work we keep n finite and from now on only consider finite discretizations of scalar fields, however the infinite case is certainly of theoretical interest and is left for future investigations.


For a given dimension d, one can define in a similar fashion the Visibility Graphs in the canonical extension class labelled IVG^c(d) by modifying condition (1): i'j' belongs to one of the d Cartesian axis which span ℝ^d and have origin in ij. Analogously, the Visibility Graphs in the FCC extension class IVG^FCC(d) are obtained by modifying again condition (1) appropriately to allow visibility in the main diagonals. Finally, again the Horizontal version follows equivalently if in the second condition we make use of HVG instead of VG.


A trivial but important remark is that ∀ I, I(H)VG_4=I(H)VG^c(2) and I(H)VG_8=I(H)VG^FCC(2). Note also that the special class IVG^c(2) has been explored recently under the name row-column visibility graph <cit.>.

Once any of these graphs has been extracted from a given matrix I, one can further compute standard topological properties on this graph using classical measures from Graph Theory <cit.> or recent metrics defined in Network Science <cit.>, which in turn might be used to provide a topological characterization of I. For instance, the degree k of a node is the number of links of that node. This allows to construct the degree matrix K∈ℕ^N× N, where K_ij is the degree of node labelled with the pair i,j. The degree distribution P(k) determines the probability of finding a node of degree k and can be straightforwardly computed from the degree matrix. In this work for concreteness we will only consider these metrics, however we should emphasize here that a large toolbox of measures could be used for feature extraction in context-dependent applications. Here we are motivated to use these very simple metrics as it has recently been proved that, in the one-dimensional case, the set of degrees is on bijection with the adjacency matrix and hence is indeed an optimal feature <cit.>.


In what follows we depict some exact results on the topology of these graphs associated to simple types of matrices which can be understood as the high order and high disorder limits of real images. From now on we only consider the Horizontal version of the visibility criteria, and we assume N→∞ to avoid border effects.



§ SOME EXACT RESULTS

Periodicity: monochromatic images and chess. 
We start by considering trivial configurations at the end of total order. For monochromatic images where I_ij=c, the IHVG_n is such that K_ij=n and thus P(n)=1 and P(k≠ n)=0. Then we can consider chess. This is a periodic lattice, where in each row the same periodic sequence is represented (black,white,black,…)≡ (1,-1,1,-1,…), except for a one-step translation in even rows. Accordingly, neglecting boundary conditions I_ij= 1 if i · j is odd and -1 otherwise. For IHVG_4 we find K_ij=8 if i · j is odd and 4 otherwise. For IHVG_8 we find K_ij=12 if i · j is odd and 8 otherwise. From this latter matrix the degree distribution is simply P(k)=1/2 for k=8,12 and zero otherwise. For other types of periodic structures it is easy to see that the degree matrix will inherit such periodicity and thus the degree distribution will only be composed by a finite number q of non-null probabilities, where q in turn is typically bounded by a function that depends on the period of the periodic structure.
 









Uncorrelated random fields. 
We then consider a limit configuration at the end of total disorder: a two-dimensional uncorrelated random field, i.e. white noise. Then, the following theorem holds for the degree distribution of IHVG_n:


Theorem.  Consider an N × N matrix with entries I_ij=ξ, where ξ is a random variable sampled from a distribution f(x) with continuous real support x∈(a,b).  Then, for n>0 and in the limit N→∞ the degree distribution of the associated IHVG_n converges to

    P(k)=     (1/n+1)(n/n+1)^k-n,    if k≥ n 
    
          0,    otherwise

A few comments are in order before presenting a proof is presented. First, note that this equation reduces, for n=2 (d=1), to the well-known result for time series of i.i.d. variables P(k)=(1/3)(2/3)^k-2 <cit.>. Second, in the specific class n=8 (equivalent to the FCC class in d=2, being this the selected version for image processing <cit.>), eq.<ref> yields

    P(k)=     (1/9)(8/9)^k-8,    if k≥ 8 
    
          0,    otherwise

Third, note that in the limit of large n we would have a continuous visibility scanning. The extension for any generic n can also be directly interpreted as a generalization to higher dimensional (discrete) scalar fields, so it is easy to show that eq.<ref> also applies to the degree distribution of (i) the canonical extension for dimension d=n/2 (i.e. only even values of n are allowed in this case), and (ii) the FCC extension for dimension d, where n=2d+2^d (i.e. for n=8,14,24,42,…). We are now ready to provide the proof of the theorem.


Proof.  The proof essentially makes use of the diagrammatic formalism introduced in <cit.> where, in the case of time series, the probability of each degree was expanded in a series of terms, each term associated to a different diagram and contributing with different amplitude.

Let us start by considering the concrete case n=8 (which describes the case implemented in our algorithm for image filtering) and we will generalize for all n thereafter. Using the jargon developed in <cit.>, a node chosen at random which has horizontal visibility of k others can be modeled as a seed (contributing with probability 𝔖) which has visibility of k-8 inner nodes (contributing with ℑ) distributed along the n=8 directions (such that direction i contributes with k_i inner nodes),
and whose visibility is finally bounded by 8 bounding nodes (contributing with probability 𝔅). 
The probability of this event can thus be formally expressed as

    P(k)=∑_{k_1,k_2… k_8}𝔖𝔅^8∏_i=1^8ℑ_k_i,

where the sum enumerates all admissible combinations of {(k_1,k_2,…,k_8)} such that ∑_i=1^8=k-8 (by construction, every node always has visibility of its boundary, here formed by n=8 nodes). It is easy to see that a possible enumeration is 

    k_i=0,1,…,k-8-∑_m=1^i-1k_m  for i=1,2,…,7;  k_8=k-8-∑_i=1^7 k_i.

Making use of the cumulative distribution  F(x)=∫_a^x f(x')dx' (with F(a)=0,  F(b)=1) and following <cit.>, geometrically it is easy to see that

    𝔖=∫_a^b f(x_0)dx_0;  𝔅=∫_x_0^b f(x)dx=1-F(x_0);


To describe the probability of finding p inner nodes ℑ_p, by construction we shall take into account that an arbitrary number r (from zero to an infinite amount) of hidden data can lie in between every pair of aligned inner nodes. Such arbitrary number of hidden data should contribute with the following amplitude

    ∑_r=0^∞∏_j=1^r ∫_a^x f(n_j)dn_j=1/1-F(x),

where we have used the properties of the cumulative distribution to find the last identity. Accordingly, the concatenation of p inner data which might have an arbitrary number of interspersed hidden data can be expressed as


    ℑ_p=∫_a^x_0f(x_1)dx_1/1-F(x_1)∏_j=1^n-1∫_x_j^x_0f(x_j+1)dx_j+1/1-F(x_j+1).

This latter calculation is easy but quite tedious. One proceeds
to integrate equation <ref> step by step and a recurrence quickly becomes evident.
One can easily prove by induction that

    ℑ_p=(-1)^p/p![ln(1-F(x_0))]^p.

We are thus ready to tackle eq. <ref>.
Taking advantage of the closure ∑_i=1^8 K_i=k-8, we first have

    ∏_i=1^8ℑ_k_i=(-1)^k-8[ln(1-F(x_0))]^k-8/∏_i=1^8(k_i)!,

so after some reordering, 

    P(k)=∑_{k_1,k_2… k_8}(-1)^k-8/∏_i=1^8(k_i)!∫_a^b f(x_0)(1-F(x_0))^8[ln(1-F(x_0))]^k-8dx_0.

Now, in this latter equation the integral is easy to compute:

    ∫_a^b f(x_0)(1-F(x_0))^8[ln(1-F(x_0))]^k-8dx_0=(-1)^k-8(k-8)!(1/9)^k-7

Consider finally the term 

    ∑_{k_1,k_2… k_8}(k-8)!/∏_i=1^8(k_i)!=(k-8)!∑_k_1=0^k-8∑_k_2=0^k-8-k_1…∑_k_7=0^k-8-∑_j=1^7 k_j1/k_1!1/k_2!…1/k_7!1/(k-8-∑_j=1^7 k_j)!=8^k-8


where the last indentity was found by iteratively applying the binomial theorem

∑_k=0^aa  k r^k=(1+r)^a. 


Altogether, we can write down explicitly for n=8

    P(k)=(1/9)(8/9)^k-8

for k≥ 8 and zero otherwise. This result is independent of f(x) as expected since HVG is an order statistic <cit.>, and coincides with eq. <ref> for n=8 (i.e. eq. <ref>).
 

We are now ready to generalize the whole derivation. For a generic n, trivially

    P(k)=∑_{k_1,k_2… k_n}(-1)^k-n/∏_i=1^n(k_i)!∫_a^b f(x_0)(1-F(x_0))^n[ln(1-F(x_0))]^k-ndx_0.

with

    ∫_a^b f(x_0)(1-F(x_0))^n[ln(1-F(x_0))]^k-ndx_0=(1/n+1)^k-n+1(-1)^k-n(k-n)!

such that

    P(k)=(1/n+1)^k-n+1∑_{k_1,k_2… k_n}(k-n)!/∏_i=1^n(k_i)!.

Finally since

    ∑_{k_1,k_2… k_n}(k-n)!/∏_i=1^n(k_i)!=n^k-n,

we find 

    P(k)=(1/n+1)^k-n+1n^k-n=(1/n+1)(n/n+1)^k-n,

what concludes the proof. ▪ 




Note that a similar result can be found much easily at the expense of using a non-rigorous heuristic argument. In the case n=8, the probability that the seed node has visibility of exactly k nodes can be expressed as the probability that there are k-8 nodes that are not bounding times the probability that after these, the boundary prevents larger visibility. Accordingly, we shall write

    P(k)=(1-P(8))^k-8P(8)

For k=8, k_i only take the value k_i=0  ∀ i=1… 8, hence this term is straightforward to compute

    P(8)=𝔖𝔅^8=∫_a^b f(x_0)[∫_x_0^b f(x)dx]^8 dx_0=1/9,  ∀ f

which then yields the correct shape for P(k):

    P(k)=(1-P(8))^k-8P(8)=(1/9)(8/9)^k-8

A similar argument can be used for a generic n, yielding

    P(k)=(1-P(n))^k-nP(n)=(1/n+1)(n/n+1)^k-n

for k≥ n and zero otherwise, in good agreement with eq. <ref>.






Finite size effects. To assess the convergence speed to eq. <ref> for finite N, we have estimated the degree distribution of IHVG_8 associated to N× N random matrices whose entries are i.i.d. uniform random variables U[0,1]. In figure <ref> we plot, in semi-log scales, the resulting (finite size) degree distributions, for different N=2^7,2^8,…,2^12. As we can see, the distributions are on excellent agreement with eq. <ref> for k≤ k_0, where the location of the cut-off value k_0 scales logarithmically with the system's size N as shown in the bottom of the figure. In other words, finite size effects only affects the tail of the distribution, which converges logarithmically fast with N. 




§ A SIMPLE APPLICATION

The results for uncorrelated random fields found in the previous section are indeed of practical interest because eq.<ref> holds independently of the noise marginal distribution f. Resorting to the contrapositive, if the degree distribution of IHVG_n deviates from eq.<ref> for some empirical field I, one can conclude that the field is not uncorrelated noise. This theorem thereby allows for the straightforward design of a randomness statistical test which would be applicable to data structure of arbitrary dimension d, where n(d)=2d if one uses the canonical extension class, or n(d)=2d+2^d in the case of FCC.



Coupled Map Lattices. To illustrate this we consider a simple application of discriminating noise from high-dimensional chaos. 
Chaotic processes display irregular and unpredictable behavior which is often confounded with randomness, however chaos is a deterministic process which indeed hides in some cases some patterns that can be extracted by appropriate techniques. The endeavor of distinguishing noise from chaos has been an area of intense research activity in the last decades <cit.> and applications have pervaded nearly every scientific discipline where complex, irregular empirical signals emerge. Here we consider spatially extended structures and thus we will be dealing with spatio-temporal chaos, i.e. chaotic behavior in space and in time, and we will explore whether if visibility graphs are able to distinguish such dynamics from simple randomness.
Let us define I(t) as a two dimensional square lattice of N^2 diffusively coupled chaotic maps which evolve in time <cit.>. In each vertex of this coupled map lattice (CML) we allocate a fully chaotic logistic map x_t+1=Q(x_t), Q(x)=4x(1-x), and the system is then spatially coupled as it follows:

    I_ij(t+1)=(1-ϵ)Q[ I_ij(t)] + ϵ/4∑_i',j'Q[ I_i'j'(t)],

where the sum extends to the Von Neumann neighborhood of ij (four adjacent neighbors)
The update is parallel and we use periodic boundary conditions. The coupling strength ϵ∈[0,1]. For ϵ=0 the system is uncoupled and the N^2 logistic maps evolve independently. For positive ϵ>0 there is a balance between the internal (chaotic) dynamics which drives a local tendency towards inhomogeneity and the diffusion term (in the right
hand of the equation one can easily recognize the discrete version of the Laplacian) which induces a global tendency
towards homogeneity in space. This balance is tuned by ϵ, acting as an effective viscosity constant, and the system evolves into different spatio-temporal dynamics as ϵ varies. For a small yet positive value of the coupling the system displays so-called Fully Developed Turbulence, a phase with incoherent spatiotemporal
chaos and high dimensional attractor <cit.>. In other words, the system evolves both temporally and spatially in a very irregular way, yet it is not totally uncorrelated. For illustration, in figure <ref> we plot, for N=200, grayscale snapshots of this system for ϵ=0 (uncoupled), ϵ=0.1 (weak coupling) and ϵ=0.7 (strong coupling) along with a 200× 200 matrix of U[0,1] i.i.d. random variables (white noise). Note that the snapshot of the uncoupled case reduces to a collection of independent and identically distributed chaotic variables with a marginal distribution that coincides with the invariant measure of the fully chaotic logistic map: the Beta distribution 𝔹(1/2,1/2)=π^-1x^-1/2(1-x)^-1/2. In other words, such a snapshot is indistinguishable from white, Beta-distributed noise, which should be then equivalent under the IHVG mapping to any type of white noise and should therefore fulfill our theorem. When ϵ>0 spatial correlations settle in and the snapshots are in theory statistically different, however this difference is only evident for large coupling.




Distinguishing noise from chaos. To explore such differences we can exploit our theorem as it follows: first, we estimate the degree distribution of the IHVG_8 of each snapshot, and compare against the theoretical equation for white noise. To account for finite size effects, it is necessary to compare the estimation of the chaotic case not just with eq.<ref> but also with a finite i.i.d. sample. 

We have generated 20 realizations of each process (random uniform noise, ϵ=0 and 0.1) and have extracted the degree distribution of IHVG_8 for each case. Sample results of these distributions can be shown in the left panel of figure <ref> along with the theoretical prediction for i.i.d (eq. <ref>). As expected, the distributions are apparently very well approximated by eq. <ref> in every case (there are strong deviations for k>35 but this is due to finite size effects as similar deviations take place for the i.i.d. white uniform noise case). To quantify potential deviations from the theory (which according to the theorem would imply non-randomness), for each case we have computed the χ^2 statistic 
    χ^2=N∑_k [P_th(k)-P_exp(k)]^2/P_th(k),

where we have taken k=8,9,…,44. Results are shown in the bottom panel of figure <ref>, showing now a clear separation between the uncorrelated cases (uncoupled chaotic maps and uniform white noise) and the weakly coupled system. This clear distinction is further confirmed in a principal component analysis (PCA) depicted in the right panel of the same figure, where each degree distribution P(k) has been projected in a two-dimensional space spanned by the first two principal components (this subspace accounts for 60% of the variability). One does not need to apply any clustering algorithm as the non-random matrices are very clearly clustered together and apart from the i.i.d. cases.




Phase diagram.  As mentioned previously, the spatio-temporal dynamics of the coupled map lattice show a rich phase diagram as we increase the coupling constant ϵ. An easy way of encapsulating and visualizing such richness in a single diagram is presented in the left panel of figure <ref>. For each ϵ, we compute the degree distribution of the associated IHVG_8. Then we compute the distance D between the degree distribution at ϵ and the corresponding result for ϵ=0 (eq. <ref>) D=∑_k |P(k)-(1/9)(8/9)^k-8|.
D acts as a scalar order parameter describing the spatial configuration of the CML, and interestingly, evidences sharp changes for the different phases, such as: for ϵ<0.12, the system develops Fully-Developed Turbulence (FDT) with weak spatial correlations. This regime shifts to a Periodic Structure (PS) for 0.12<ϵ<0.27. This regime then parsimoniously shifts into a phase with spatially Coherent Structures (CS), which ultimately break down for ϵ>0.88 in favor of periodic patterns. For 0.88<ϵ<1 the spatial structure shows a mix between CS and PS. We conclude that the degree distribution of the IHVG_8 captures this rich spatial structure, something confirmed via principal component analysis in the right panel of figure <ref>.





§ DISCUSSION

This framework allows the possibility of describing discretized scalar fields of arbitrary origin in a combinatorially compact fashion, and enables using the tools of graph theory and network science for the practical description and classification of spatially-extended data structures. For the sake of exposition and concreteness, in this work we have only used a couple of graph measures (degree matrix and degree distribution) which can be argued that were optimal in the one-dimensional case <cit.>, but it should be highlighted that this method is much more general and allows to extract from these graphs any desired property.

For d=1 the method was naturally designed for the task of time series analysis, and has been exploited accordingly and extensively in the last years -both from a theoretical point of view and for applications- as was acknowledged in the introduction section. Here we have presented a natural extension of these algorithms to deal with (discretized) scalar fields of arbitrary dimension, along with a few exact results on simple -yet relevant- cases. From a mathematical point of view, the task of characterizing the graphs in these extension classes provide a wide range of challenging open questions, which could parallel recent advancements in the one-dimensional case <cit.>. Now, what are the potential applications of this framework?




For d=2 (either using the canonical or FCC extension classes, or the order-n class), a plethora of applications emerge, here we only enumerate and discuss a few: (i) Image Processing: a (grayscale) image is just a discrete scalar field. Once we extract the visibility graphs of a given image, can we use the topological properties of this graph to build feature vectors which can feed automatic classifiers for several statistical learning tasks involving images <cit.>? Can we define the distance between two images using graph kernels <cit.> on the associated visibility graphs? 
(ii) Physics of Interfaces: can we provide a topological characterization of fractal surface growth <cit.>? Can we -for instance- account for spatial self-similar structures much in the same way the Hurst exponent of fractional Brownian motion was estimated with visibility graphs <cit.> (a preliminary analysis via row-column visibility graphs has partly addressed this issue recently <cit.>). Furthermore, can we apply this methodology in biologically-relevant problems and beyond, for instance to classify tumoral or calli surfaces? (iii) Urban Planning: can we automatically cluster cities by only resorting to combinatorial properties extracted from their visibility graphs? And can we link such emerging clusters with architectural, historical or cultural properties of cities? (iv) Random Matrix theory: Is there a visibility graph characterization of different random matrix ensembles?


To illustrate the potential applicability of the method to the case of tumor description, in the left panel of figure <ref> we plot the degree distribution of the IHVG_8 associated to three atomic force microscopy (AFM) images (94× 94 after grayscale preprocessing) of normal, immortal (premalignant) and cancer (malignant) cervical epithelial cells <cit.>. This very preliminary evidence suggests that the carcinogenesis transition normal → premalignant → cancer is paralleled in graph space by a systematic deviation of the degree distribution from the i.i.d. case. In the right panel of the same figure we plot the degree distribution associated to IVG_8, whose tails have been fitted to exponential functions ∼exp(-λ k) finding λ_normal<λ_immortal<λ_cancer <cit.>. These are of course very preliminary results given simply for illustration, and future research should confirm their accuracy and their potential use for carcinogenesis description and early detection.





The most exciting application for higher dimensions d≥ 2 is perhaps on describing the spatial structure of generic energy landscapes <cit.> V: x∈ℝ^d→ℝ, where d is the number of degrees of freedom. Typically, these fields describe an energy function whose minimum is associated to the macroscopic behavior of many-body systems, and play a major role in physics and chemistry. The structure of these fields is however rather messy. As a matter of fact, in spin glasses and other disordered systems their macroscopic properties do not necessarily relate directly to a configuration of minimal energy as the system gets trapped in local, metastable minima of this energy surface: in this sense the spatial distribution and overall structure of these minima (stationary points) gives valuable information on the system dynamical evolution. These energy surfaces are also of great interest in chemistry (Kramer's reaction rate theory for the thermally activated escape from metastable
states) and high energy physics (e.g. local minima of supersymmetric energy landscape corresponds to the field theory vacuum). The formalism presented here would enable the description of such energetic landscapes, opening a thread of questions such as: Can we classify different types of field theories only using combinatorial criteria on their energy landscapes? What is the spatial distribution of stationary points of different canonical disordered systems in the light of this new method?


To conclude, we hopefully made the case that to encode spatially extended structures in a combinatorial fashion is an enterprise that opens exciting theoretical questions as well as applications. The approach presented here is promising and there exist several possible avenues for future research, and we hope that these methods spark interest in some of these communities accordingly.




We thank I. Sokolov for granting permission to reproduce the images on normal, immortal and cancer cells. LL acknowledges funding from EPSRC Fellowship EP/P01660X/1.























10
GT B. Bollobas, Modern Graph Theory (Springer, 1998).
NS M. Newman, Networks: and introduction (Oxford University Press, 2010).
Luque_Theorem B. Luque, L. Lacasa, Canonical horizontal visibility graphs are uniquely determined by their degree sequence,
Eur. Phys. J. B Sp. Top. (in press).
libro_chaos H. Kantz, T. Schreiber, Nonlinear time series analysis (Cambridge University Press).
PNAS L. Lacasa, B. Luque, F.J. Ballesteros, J. Luque, and J.C. Nuno, From time series to complex networks: the visibility graph, Proc. Natl. Acad. Sci. USA 105, 13 (2008).
image_processing J. Iacovacci and L. Lacasa, Visibility graphs: a combinatorial framework for image processing (in preparation).
PRE B. Luque, L. Lacasa, J. Luque, F.J. Ballesteros, Horizontal visibility graphs: exact results for random time series, Phys. Rev. E 80, 046103 (2009).
severini S. Severini, G. Gutin, T. Mansour, A characterization of horizontal visibility graphs and combinatorics on words, Physica A 390, 12  (2011) 2421-2428.
flajo P. Flajolet and M. Noy, Analytic combinatorics of non-crossing configurations, Discrete Math. 204 (1999) 203-229.
nonlinearity L. Lacasa, On the degree distribution of horizontal visibility graphs associated to Markov processes and dynamical systems: diagrammatic and variational approaches, Nonlinearity 27, 2063-2093 (2014).
nonstationary L.Lacasa and R. Flanagan, Time reversibility from visibility graphs of non-stationary processes,
Phys. Rev. E 92, 022817 (2015).
EPL L. Lacasa, B. Luque, J. Luque and J.C. Nuno, The Visibility Graph: a new method for estimating the Hurst exponent of fractional Brownian motion, EPL 86, 30001 (2009). 
rowcolumn X Qin, P Xue, L Xin-Li, M Stephen, Y Hui-Jie, J Yan, W Jian-Yong, Z. Quin-Jung, Row-column visibility graph approach to two-dimensional landscapes, Chinese Physics B 23, 7 (2014).
original A. Turner, M. Doxa, D. O'sullivan, and A. Penn, From isovists to visibility graphs: a methodology for the analysis of architectural space. Environment and Planning B: Planning and design, 28(1), 103-121 (2001).
CML K. Kaneko, Overview of Coupled Map Lattices, Chaos 2, 3(1992).
kernel S.V. N. Vishwanathan, N.N. Schraudolph, R. Kondor and K.M. Borgwardt, Graph kernels, Journal of Machine Learning Research 11 (2010) pp.1201-1242.
barabasiA.L. Barabasi and H.E. Stanley, Fractal Concepts in Surface Growth (Cambridge University Press, 1995).
PEL1 D. Wales, Energy Landscapes : Applications to Clusters,
Biomolecules and Glasses (Cambridge University Press, 2004).




jns B. Luque, L. Lacasa, F. Ballesteros, A. Robledo,  Analytical properties of horizontal visibility graphs in the Feigenbaum scenario, Chaos 22, 1 (2012) 013109.
quasi B. Luque, A. Núñez, F. Ballesteros, A. Robledo, Quasiperiodic Graphs: Structural Design, Scaling and Entropic Properties, Journal of Nonlinear Science 23, 2, (2012) 335-342.
pre2013 A.M. Núñez, B. Luque, L. Lacasa, J.P. Gómez, A. Robledo, Horizontal Visibility graphs generated by type-I intermittency, Phys. Rev. E, 87 (2013) 052801.

physics3 A. Aragoneses, L. Carpi, N. Tarasov, D.V. Churkin, M.C. Torrent, C. Masoller, and S.K. Turitsyn, Unveiling Temporal Correlations Characteristic of a Phase Transition in the Output Intensity of a Fiber Laser, Phys. Rev. Lett. 116, 033902 (2016).
    
fluiddyn0 M. Murugesana and R.I. Sujitha1, Combustion noise is scale-free: transition from scale-free to order at the onset of thermoacoustic instability, J. Fluid Mech. 772 (2015).
fluiddyn1A. Charakopoulos, T.E. Karakasidis, P.N. Papanicolaou and A. Liakopoulos, The application of complex network time series analysis in turbulent heated jets, Chaos 24, 024408 (2014).
fluiddyn2 P. Manshour, M.R. Rahimi Tabar and J. Peinche, Fully developed turbulence in the view of horizontal visibility graphs, J. Stat. Mech. (2015) P08031.


physics2 RV Donner, JF Donges, Visibility graph analysis of geophysical time series: Potentials and possible pitfalls, Acta Geophysica 60, 3 (2012).





suyal V. Suyal, A. Prasad, H.P. Singh, Visibility-Graph Analysis of the Solar Wind Velocity, Solar Physics 289, 379-389 (2014)
Zou Y. Zou, R.V. Donner, N. Marwan, M. Small, and J. Kurths, Long-term changes in the north-south asymmetry of solar activity:
a nonlinear dynamics characterization using visibility graphs, Nonlin. Processes Geophys. 21, 1113-1126 (2014).

physio1 J.F. Donges, R.V. Donner and J. Kurths, Testing time series irreversibility using complex network methods, EPL 102, 10004 (2013).
meditation_VG S. Jiang, C. Bian, X. Ning and Q.D.Y. Ma, Visibility graph analysis on heartbeat dynamics of meditation training, Appl. Phys. Lett. 102 253702 (2013).




neuro M Ahmadlou, H Adeli, A Adeli, New diagnostic EEG markers of the Alzheimer's disease using visibility graph, J. of Neural Transm. 117, 9 (2010).


ryan1 R. Flanagan and L. Lacasa, Irreversibility of financial time series: a graph-theoretical approach,
Physics Letters A 380, 1689-1697 (2016)
multivariate L. Lacasa, V. Nicosia, V. Latora, Network Structure of Multivariate Time Series, Sci. Rep. 5, 15508 (2015)


motifs J. Iacovacci and L. Lacasa, Sequential visibility-graph motifs, Phys. Rev. E 93, 042309 (2016)

cancer_paper M.E. Dokukin, N.V. Guz, C.D. Woodworth, and I. Sokolov, Emerging of fractal geometry on surface of human cervical epithelial cells during progression towards cancer, New journal of physics 17,3 (2015).
cancer2 M.E. Dokukin, N.V. Guz, R. M. Gaikwad, C.D. Woodworth, and I. Sokolov, Cell Surface as a Fractal: Normal and Cancerous Cervical Cells Demonstrate Different Fractal Behavior of Surface Adhesion Maps at the Nanoscale
Phys. Rev. Lett. 107, 028101 (2011).





