David Coudert

André Nusser

Laurent Viennot

Enumeration of far-apart pairs by decreasing distance for faster hyperbolicity computation *

Keywords:

Hyperbolicity is a graph parameter which indicates how much the shortest-path distance metric of a graph deviates from a tree metric. It is used in various fields such as networking, security, and bioinformatics for the classification of complex networks, the design of routing schemes, and the analysis of graph algorithms. Despite recent progress, computing the hyperbolicity of a graph remains challenging. Indeed, the best known algorithm has time complexity O(n 3.69 ), which is prohibitive for large graphs, and the most efficient algorithms in practice have space complexity O(n 2 ). Thus, time as well as space are bottlenecks for computing the hyperbolicity.

In this paper, we design a tool for enumerating all far-apart pairs of a graph by decreasing distances. A node pair (u, v) of a graph is far-apart if both v is a leaf of all shortest-path trees rooted at u and u is a leaf of all shortest-path trees rooted at v. This notion was previously used to drastically reduce the computation time for hyperbolicity in practice. However, it required the computation of the distance matrix to sort all pairs of nodes by decreasing distance, which requires an infeasible amount of memory already for medium-sized graphs. We present a new data structure that avoids this memory bottleneck in practice and for the first time enables computing the hyperbolicity of several large graphs that were far out of reach using previous algorithms. For some instances, we reduce the memory consumption by at least two orders of magnitude. Furthermore, we show that for many graphs, only a very small fraction of far-apart pairs has to be considered for the hyperbolicity computation, explaining this drastic reduction of memory.

As iterating over far-apart pairs in decreasing order without storing them explicitly is a very general tool, we believe that our approach might also be relevant to other problems.

Introduction

This paper aims at computing the hyperbolicity of graphs whose size ranges from tens of thousands to millions of nodes. The hyperbolicity is a parameter of a metric space generalizing the idea of Riemannian manifolds with negative curvature. When considering the metric of a graph, the hyperbolicity measures, to some extent, how much the metric of the graph deviates from a tree metric. This parameter was first introduced by Gromov in the context of automatic groups [START_REF] Gromov | Hyperbolic groups[END_REF] in relation with their Cayley graphs.

Hyperbolicity has received great attention in computer science in the last decades as it seems to capture important properties of several large practical graphs such as Internet [START_REF] Shavitt | On the curvature of the internet and its usage for overlay construction and distance estimation[END_REF], the Web [START_REF] Munzner | Visualizing the structure of the world wide web in 3d hyperbolic space[END_REF] and databases relations [START_REF] Walter | On interactive visualization of high-dimensional data using the hyperbolic plane[END_REF]. It is also used to classify complex networks [START_REF] Abu | Metric tree-like structures in real-life networks: an empirical study[END_REF][START_REF] Alrasheed | Core-periphery models for graphs based on their δhyperbolicity: An example using biological networks[END_REF][START_REF] Kennedy | On the hyperbolicity of large-scale networks and its estimation[END_REF] and was proposed as a measure of how much a network is "democratic" [START_REF] Albert | Topological implications of negative curvature for biological and social networks[END_REF][START_REF] Borassi | Hyperbolicity measures democracy in real-world networks[END_REF]. Formal relationships between Gromov hyperbolicity and the existence of a core (a subset of nodes intersecting a constant fraction of all the shortest-paths) have also been investigated [START_REF] Chepoi | Core congestion is inherent in hyperbolic networks[END_REF]. Reciprocally, the existence of a core has been shown to be inherent to any hyperbolic network [START_REF] Chepoi | Core congestion is inherent in hyperbolic networks[END_REF]. Furthermore, small hyperbolicity has tractability implications and measuring hyperbolicity has applications in routing [START_REF] Boguñá | Sustaining the Internet with hyperbolic mapping[END_REF][START_REF] Chepoi | Additive spanners and distance and routing labeling schemes for hyperbolic graphs[END_REF][START_REF] Krauthgamer | Algorithms on negatively curved spaces[END_REF], approximating other graph parameters [START_REF] Chepoi | Diameters, centers, and approximating trees of delta-hyperbolic geodesic spaces and graphs[END_REF][START_REF] Bhaskar Dasgupta | Effect of Gromov-hyperbolicity parameter on cuts and expansions in graphs and some algorithmic implications[END_REF] and bioinformatics [START_REF] Chakerian | Computational tools for evaluating phylogenetic and hierarchical clustering trees[END_REF][START_REF] Dress | Basic Phylogenetic Combinatorics[END_REF]. See [START_REF] Abu | Metric tree-like structures in real-life networks: an empirical study[END_REF][START_REF] Feodor | Tree-like structures in graphs: A metric point of view[END_REF] for recent surveys.

Computing the hyperbolicity is often a prerequisite in the above applications. As hyperbolicity can be defined by a simple 4-point condition, it can be naively computed in O(n 4 ) time. As far as we know, the best theoretical algorithm [START_REF] Fournier | Computing the gromov hyperbolicity of a discrete metric space[END_REF] has time complexity O(n 3.69 ). Although its complexity is o(n 4 ), it is still supercubic and the algorithm appears to be impractical for graphs with a few tens of thousands of nodes. On the lower bounds side it was shown that under the Strong Exponential Time Hypothesis [START_REF] Impagliazzo | Which problems have strongly exponential complexity[END_REF] hyperbolicity cannot be computed in subquadratic-time, even for sparse graphs [START_REF] Borassi | Into the square: On the complexity of some quadratic-time solvable problems[END_REF][START_REF] Coudert | Recognition of C 4 -free and 1/2-hyperbolic graphs[END_REF][START_REF] Fournier | Computing the gromov hyperbolicity of a discrete metric space[END_REF].

The only practical algorithms that can manage larger graphs [START_REF] Borassi | On computing the hyperbolicity of real-world graphs[END_REF][START_REF] Cohen | On computing the gromov hyperbolicity[END_REF] enumerate all pairs of nodes by decreasing distance. For each pair, each 4-tuple obtained with a previous pair is tested with regard to the 4-point condition defining hyperbolicity. Each test of a 4-tuple provides a lower bound of hyperbolicity. This approach allows to stop the enumeration as soon as the distance of the scanned pair equals twice the best hyperbolicity lower bound found so far. As it scans a portion of all 4-tuples, its worst case complexity is O(n 4 ) but it appears much faster in practice as first scanning pairs with large distances allows to find good lower bounds early in the enumeration. A main optimization for further reducing the number of pairs scanned, consists in considering only far-apart pairs, that is, pairs such that no neighbor of one node is further apart from the other node, see Section 2 for a formal definition. It can be proven that the 4-point condition defining hyperbolicity holds on all 4-tuples if it holds on 4-tuples made up of two such far-apart pairs [START_REF] Noguès | δ-hyperbolicité et graphes[END_REF][START_REF] Abel | Quelques propriétés topologiques des graphes et applications à internet et aux réseaux[END_REF]. The main bottleneck of this method lies in its inherent quadratic space usage: all far-apart pairs and all-pair distances are stored in the current implementations of this approach. Computing the hyperbolicity of practical graphs with millions of nodes thus remains a great challenge.

Our approach

To make progress towards this challenge, we propose to enumerate far-apart pairs by decreasing distance without computing all-pair distances. The key of our approach is to first compute all eccentricities, that is, for each node, what is the largest distance from it. Computing all eccentricities is feasible in practice. Indeed, efficient algorithms for computing all eccentricities [START_REF] Takes | Computing the eccentricity distribution of large graphs[END_REF][START_REF] Feodor | Revisiting radius, diameter, and all eccentricity computation in graphs through certificates[END_REF][START_REF] Li | Exacting eccentricity for small-world networks[END_REF] came along in a line of research for improving diameter computation of real world graphs [START_REF] Akiba | An exact algorithm for diameters of large real directed graphs[END_REF][START_REF] Crescenzi | On computing the diameter of real-world undirected graphs[END_REF][START_REF] Borassi | Fast diameter and radius BFS-based computation in (weakly connected) real-world graphs: With an application to the six degrees of separation games[END_REF][START_REF] Feodor | Revisiting radius, diameter, and all eccentricity computation in graphs through certificates[END_REF].

• As the prime application of our algorithm for iterating over far-apart pairs, we significantly reduce the memory consumption when computing the graph hyperbolicity. The memory reduction is at least two orders of magnitude for some instances.

• This drastic memory reduction enables us to compute the hyperbolicity of many large graphs for the first time.

• Due to the significance of far-apart pairs in a graph (e.g., they are the defining nodes for radius, eccentricities, diameter, . . . ), we believe that our contribution of a far-apart pair iterator is also relevant in other settings.

Organization

Definitions and notations used in this paper are introduced in Section 2. We then present our far-apart pair iterator in Section 3. Section 4 is devoted to hyperbolicity. We recall its definition and review the best know algorithmic results on this parameter. We then present our memory efficient algorithm based on the proposed far-apart pair iterator. In Section 5, we report on the experimental evaluation of the algorithms presented in this paper on various graphs. In particular, we compare our algorithm for computing hyperbolicity with the state-of-the-art algorithm [START_REF] Borassi | On computing the hyperbolicity of real-world graphs[END_REF]. We conclude this paper in Section 6 with some directions for future research.

Definitions and notations

We use standard graph terminology [START_REF] John | Graph theory with applications[END_REF][START_REF] Diestel | Graph Theory[END_REF]. All graphs considered in this paper are finite, undirected, connected, unweighted and simple. The graph G = (V, E) has n = |V | nodes and m = |E| edges. The open neighborhood N G (S) of a set S ⊆ V consists of all nodes in V \ S with at least one neighbor in S.

Given two nodes u and v, a uv-path of length ℓ ≥ 0 is a sequence of nodes (u = v 0 v 1 . . . v ℓ = v), such that {v i , v i+1 } is an edge for every i. In particular, a graph G is connected if there exists a uv-path for all pairs u, v ∈ V , and in such a case the distance d G (u, v) is defined as the minimum length of a uv-path in G. When G is clear from the context, we write d (resp. N ) instead of d G (resp. N G ). The eccentricity ecc(u) of a node u is the maximum distance between u and any other node v ∈ V , i.e., ecc(u) = max v∈V d(u, v). The maximum eccentricity is the diameter diam(G) and the minimum eccentricity is the radius rad(G).

The notion of far-apart pairs of nodes has been introduced [START_REF] Abel | Quelques propriétés topologiques des graphes et applications à internet et aux réseaux[END_REF][START_REF] Noguès | δ-hyperbolicité et graphes[END_REF] to reduce the number of 4-tuples to consider in the computation of the hyperbolicity (see Section 4). Roughly, we say that two nodes u, v ∈ V are far-apart if for all w ∈ V neither u lies on a shortest path from w to v, nor v lies on a shortest path from u to w. More formally, we have:

Definition 1. In a graph G = (V, E), node u is far from node v, or v-far, if for any neighbor w of u, we have d(v, w) ≤ d(v, u). The pair u, v of nodes is far-apart if u is v-far and v is u-far.
The number of far-apart pairs in a graph can be orders of magnitude smaller than the total number of pairs. For instance, a p × q grid has only 2 far-apart pairs. On the other hand, all pairs in a clique graph are far-apart.

The set of all far-apart pairs can be determined in time O(nm) in unweighted graphs through breadth-first search (BFS), and the interested reader will find in Appendix A a discussion on several time and space complexity trade-offs for determining far-apart pairs. We now present some interesting properties of far nodes.

Lemma 1. For v ∈ V , node u ∈ V is v-far if and only if u is a leaf of all shortest path trees rooted at v.

Proof. Clearly, a non-leaf node u of a shortest path tree rooted at v has a neighbor at distance d(v, u) + 1 and cannot be v-far. Reciprocally, if u is not v-far, it has a neighbor w satisfying d(v, w) = d(v, u) + 1. Modifying any shortest path tree by setting u as the parent of w yields a valid shortest path tree where u is not a leaf.

From Lemma 1, we also get the following immediate results.

Corollary 1. For each u ∈ V , any v ∈ V such that d(u, v) = ecc(u) is u-far. Corollary 2. For any u, v ∈ V , if d(u, v) = diam(G), then the pair (u, v) is far-apart.
In particular, Corollary 1 implies that for any node u, the set of u-far nodes is non-empty. Interestingly, knowing the far nodes of a node u and their distance to u, it is possible to scan distant nodes in a sort of backward BFS as described in Appendix B.

Iterator over far-apart pairs

In this section, we engineer a data structure and algorithms to determine the set of far-apart pairs and return these pairs sorted by decreasing distances. Our objective is to provide an iterator that determines the next pair to yield on the fly, that postpones computations as much as possible, and which has an acceptable memory consumption.

The high-level idea of our data structure is that we maintain for each distance up to the diameter a dictionary that, for each node, contains one of the following:

• no entry: In this case, this node is not involved in any far-apart pair at this distance.

• a list of vertices: These are all the far nodes for this node at this distance. We only have to check whether it is also a far-apart pair.

• an empty list: This is a marker that we still have to compute the far nodes for this node and then proceed according to the previous case.

We can iterate over the dictionaries in decreasing distance and compute the far-apart nodes lazily. This lazy computation and the possibility to delete entries that are not needed anymore are the reasons why we are able to save memory in practice. In the remainder of this section, we present the details of our data structure.

Data structure. To store and organize data, we use an array F indexed by distances in range from 1 to diam(G), so that cell F d contains data related to far-apart pairs at distance d. More precisely, F d is a hash map associating to a node u ∈ V the subset F d u of u-far nodes at distance d from u. The subset F d u can be implemented using either a set data structure allowing to answer a membership query in O(1) time, or an array of size |F d u | whose elements are sorted according to any total ordering of the nodes to enable membership queries in time O(log |F d u |). We assume, as it is the case for most modern programming languages, that it is possible to visit the elements of a hash map in an arbitrary order (e.g., insertion order). We have the same assumption for set data structures.

Initialization. Recall that our aim is to iterate over far-apart pairs by non-increasing distances and to postpone computation as much as possible. To this end, we initially store in F , for each node u with eccentricity ecc(u), an empty set of u-far nodes in F d with d = ecc(u). This empty set will serve as an indicator to trigger the effective computation of the set of u-far nodes the first time it is needed (recall that by Corollary 1, this set is non-empty).

Clearly, this initialization procedure requires the knowledge of the eccentricities of all nodes. Fortunately, although the determination of the eccentricities has worst-case time complexity in O(nm), practically efficient algorithms have been proposed to perform this task [START_REF] Feodor | Revisiting radius, diameter, and all eccentricity computation in graphs through certificates[END_REF][START_REF] Takes | Computing the eccentricity distribution of large graphs[END_REF]. These algorithms maintain upper and lower bounds on the eccentricity of each node and improve these bounds by computing distances from a few well chosen nodes until all gaps are closed. In practice, these algorithms are orders of magnitude faster than a naive algorithm performing a BFS from each node.

Filling. Each time we compute distances from a node u that has not been considered before, the corresponding sets F d u for all d can be inserted in F . This is done in particular while computing the eccentricities during the initialization of the data structure.

Next. We define a function next(F ) that yields the next far-apart pair in the ordering. The use of a total ordering on the nodes allows to ensure that a far-apart pair (u, v) is returned only once, i.e., when u < v. During these operations, if we encounter an empty set F d x , we know from the initialization step that x has never been considered before and that ecc(x) = d. So, we compute distances from x and store the x-far nodes in F via the filling step. That is, we insert in F the corresponding sets

F d ′ x for each d ′ ≤ d.
Improvements. Depending on the usage of this data structure, some simple improvements can be done, such as:

1. When iterating over the nodes in F d u , each time a node w ∈ F d u is found such that u is not w-far, we can remove w from F d u . This avoids checking twice if a pair is far-apart, and, at the end of the iterations on F d u , this set will contain a node w only if u is wfar. Observe that this improvement is interesting only if an algorithm has to iterate more than once over the far-apart pairs. We use this improvement in our algorithm for computing hyperbolicity.

2. When all nodes in F d have been considered, one can delete F d to reduce the memory consumption if appropriate with the usage. Similarly, one can avoid storing F d if only far-apart pairs at distance d ′ > d are requested, as is the case for instance to enumerate the diameters of a graph or when computing hyperbolicity when a lower-bound has been found.

The time complexity for iterating over all far-apart pairs sorted by non-increasing distance is in expected time O(nm) when using sets for F d u (resp., O(nm + n 2 log n) when using arrays). Indeed, the algorithm computes BFS distances from each node of the graph (some during initialization, and others during queries). Furthermore, the query time to report all far-apart pairs involving a node

u is in O(n) since | ∪ ecc(u) d=1 F d u | = O(n) and checking if u ∈ F d w requires time O(1) (resp., O(log n)).
The sorting operation over the far-apart pairs is implicit and thus adds no extra cost. The space complexity of the data structure is in O(n 2 ). However, we will see in Section 5 that it is much smaller in practice. Figure 1: Hyperbolicity in a metric space (Fig. 1a) and in a tree (Fig. 1b).

Use Case: Gromov Hyperbolicity

We now present an interesting use case for our far-apart pair iterator with the computation of the Gromov hyperbolicity of a graph [START_REF] Gromov | Hyperbolic groups[END_REF].

Definitions

A metric space (V, d) is δ-hyperbolic if any triangle (u, v, x) in this metric space is δ-thin, that is if any point y along one of the geodesics (i.e., shortest path in a metric space) uv, ux or vx is at distance at most δ from one of the two other geodesics (see Fig. 1a for an illustration). This notion has been extended to graphs [START_REF] Gromov | Hyperbolic groups[END_REF]. More precisely, a metric space (V, d G ) is a tree metric if there exists a distance-preserving mapping from V to the nodes of an edge-weighted tree. If so, the graph G is said to be 0-hyperbolic because it satisfies the 4-points condition below with parameter δ = 0. Introducing some slack δ allows to measure how much the metric of a graph deviates from that of a tree. This slack δ is called the hyperbolicity of the graph: Definition 2 (4-points Condition, [START_REF] Gromov | Hyperbolic groups[END_REF]). Let G be a connected graph. For every 4-tuple u, v, x, y of nodes of G, we define δ(u, v, x, y) as half of the difference between the two largest sums among

S 1 = d(u, v) + d(x, y), S 2 = d(u, x) + d(v, y), and S 3 = d(u, y) + d(v, x).
The hyperbolicity of G, denoted by δ(G), is equal to max u,v,x,y∈V (G) δ(u, v, x, y). Moreover, we say that G is δ-hyperbolic whenever δ ≥ δ(G).

Other characterizations exist for 0-hyperbolic graphs and yield other definitions for the hyperbolicity δ of a graph that differ only by a small constant factor [START_REF] Bermudo | Gromov hyperbolic graphs[END_REF][START_REF] De | Sur les groupes hyperboliques d'après Mikhael Gromov[END_REF][START_REF] Gromov | Hyperbolic groups[END_REF].

Let us illustrate the behavior of the 4-points condition on trees. Without loss of generality, consider the example of Fig. 1b, where s u , s v , s x , s y and η denote the lengths of the corresponding paths in the figure. We have

S 1 = d(u, v) + d(x, y) = s u + s v + s x + s y + 2η, S 2 = d(u, x) + d(v, y) = s u + s v + s x + s y , and S 3 = d(u, y) + d(v, x) = s u + s v + s x + s y + 2η.
The two largest sums S 1 and S 3 are such that S 1 = S 3 , and so δ(u, v, x, y) = 0. The same holds for any quadruple in a tree and so trees are 0-hyperbolic.

From the 4-points condition above, it is straightforward to compute graph hyperbolicity in Θ(n 4 )-time. In theory, it can be decreased to O(n 3.69 ) by using a clever (max, min)-matrix product [START_REF] Fournier | Computing the gromov hyperbolicity of a discrete metric space[END_REF]; however, in practice, the best-known algorithms still run in O(n 4 )-time [START_REF] Borassi | On computing the hyperbolicity of real-world graphs[END_REF][START_REF] Cohen | On computing the gromov hyperbolicity[END_REF]. Graphs with small hyperbolicity can be recognized faster. In particular, 0-hyperbolic graphs coincide with block graphs, that are graphs whose biconnected components are complete subgraphs [START_REF] Bandelt | Distance-hereditary graphs[END_REF][START_REF] Howorka | On metric properties of certain clique graphs[END_REF]. Hence, deciding whether a graph is 0-hyperbolic can be done in time O(n + m). The latter characterization of 0-hyperbolic graphs follows from a more general result stating that the hyperbolicity of a graph is the maximum hyperbolicity of its biconnected components (see e.g. [START_REF] Cohen | On computing the gromov hyperbolicity[END_REF] for a proof). More recently, it has been proven that the recognition of 1 2 -hyperbolic graphs is computationally equivalent to deciding whether there exists a chordless cycle of length 4 in a graph [START_REF] Coudert | Recognition of C 4 -free and 1/2-hyperbolic graphs[END_REF]. The latter problem can be solved in deterministic O(n 3.26 )-time [START_REF] François | Faster algorithms for rectangular matrix multiplication[END_REF] and in randomized O(n 2.373 )-time [START_REF] Vassilevska Williams | Finding four-node subgraphs in triangle time[END_REF] by using fast matrix multiplication.

Several pre-processing methods for reducing the size of the input graph have been proposed. In particular, it has been proved [START_REF] Abel | Quelques propriétés topologiques des graphes et applications à internet et aux réseaux[END_REF] that the hyperbolicity of G is equal to the maximum of the hyperbolicity of the graphs resulting from both a modular [START_REF] Gallai | Transitiv orientierbare graphen[END_REF][START_REF] Habib | A survey of the algorithmic aspects of modular decomposition[END_REF] or a split [START_REF] William | A combinatorial decomposition theory[END_REF][START_REF] William | Decomposition of directed graphs[END_REF] decomposition of G. These decompositions can be computed in linear time [START_REF] Charbit | Linear time split decomposition revisited[END_REF] and algorithms with time complexity in O(mw(G) 3 • n + m) and O(sw(G) 3 • n + m) when parameterized respectively by the modular width mw(G) and the split width sw(G) have been proposed [START_REF] Coudert | Fully polynomial FPT algorithms for some classes of bounded clique-width graphs[END_REF]. Moreover, one can use modified versions of the atoms of a decomposition of G by clique-minimal separators [START_REF] Endre | Decomposition by clique separators[END_REF][START_REF] Berry | An introduction to clique minimal separator decomposition[END_REF][START_REF] Coudert | Revisiting Decomposition by Clique Separators[END_REF] as a pre-processing for hyperbolicity [START_REF] Cohen | Applying cliquedecomposition for computing gromov hyperbolicity[END_REF]. This decomposition can be obtained in time O(nm). Although we do not use these pre-processing methods in this work, they are complementary to our approach and it is thus possible to combine them with our approach.

Previous algorithm

In this section, we recall the state of the art algorithm for computing hyperbolicity [START_REF] Borassi | On computing the hyperbolicity of real-world graphs[END_REF]. This algorithm, referenced as Borassi et al. algorithm in the sequel, improves upon a former algorithm [START_REF] Cohen | On computing the gromov hyperbolicity[END_REF] by adding pruning techniques to further reduce the number of 4-tuples to consider. To the best of our knowledge, these algorithms are the only one enabling to compute the exact hyperbolicity of graphs with up to 50 000 nodes. These algorithms are based on the following lemmas.

Lemma 2 ( [START_REF] Abel | Quelques propriétés topologiques des graphes et applications à internet et aux réseaux[END_REF][START_REF] Cohen | On computing the gromov hyperbolicity[END_REF]). Let G be a connected graph. There exist two far-apart pairs (u, v) and (x, y) satisfying δ(u, v, x, y) = δ(G).

Lemma 3 ([23]

). For every 4-tuple u, v, x, y of nodes of a connected graph G, we have δ(u, v, x, y) ≤ min a,b∈{u,v,x,y} d(a, b). Furthermore, if S 1 = d(u, v) + d(x, y) is the largest of the sums defined in Definition 2 (which can be assumed w.l.o.g.), we have δ(u, v, x, y) ≤ 1 2 min {d(u, v), d(x, y)}. For the sake of completeness, we quickly give the proof of Lemma 3. Without loss of generality the second largest sum is

S 2 = d(u, x) + d(v, y). By triangle inequality, we have d(u, v) ≤ d(u, x) + d(x, y) + d(y, v) = S 2 + d(x, y), yielding S 1 -S 2 ≤ 2 d(x, y). We can similarly obtain S 1 -S 2 ≤ 2 d(u, v).
The key idea of the state of the art algorithms [START_REF] Borassi | On computing the hyperbolicity of real-world graphs[END_REF][START_REF] Cohen | On computing the gromov hyperbolicity[END_REF] is to visit the most promising 4-tuples first, that is, those made of pairs of far-apart nodes at largest distance, and to stop computation as soon as the bounds of Lemma 3 are reached. These algorithms thus need to iterate over far-apart pairs ordered by decreasing distances.

More precisely, Algorithm 1 below (see also [START_REF] Borassi | On computing the hyperbolicity of real-world graphs[END_REF]) iterates first over the far-apart pairs sorted by non-increasing distances. Then, given the i-th far-apart pair (x i , y i ), it iterates over the previous far-apart pairs (v, w), such that d(v, w) ≥ d(x i , y i ) in order to consider quadruples (v, w, x i , y i ) such that S 1 = d(v, w) + d(x i , y i ) is the largest sum. Note that such pairs (v, w) have been considered previously and satisfy (v, w) = (x j , y j ) for some j < i. This ensures by Lemma 3 that as soon Algorithm 1: Borassi et al. algorithm for computing the hyperbolicity [START_REF] Borassi | On computing the hyperbolicity of real-world graphs[END_REF] Input: F = ({x 1 , y 1 }, . . . , {x N , y N }): an ordered list of far-apart pairs; d: the distance matrix.

1 δ L ← 0 2 mates[v]← ∅ for each v 3 for i ∈ [1, N ] do 4 if d(x i , y i ) ≤ 2δ L then 5 return δ L 6 (acceptable, valuable) ← computeAccVal(x i , y i , δ L , d) 7 for v ∈ valuable do 8 for w ∈ mates[v] do 9 if w ∈ acceptable then 10 δ L ← max{δ L , δ(x i , y i , v, w)} 11 add y i to mates[x i ] 12 add x i to mates[y i ] 13 return δ L as d(x i , y i ) ≤ 2δ L ,
where δ L is the current best solution, we determined the hyperbolicity of the graph. So, the hyperbolicity δ L of the graph is then returned in Line 5. This ensures by Lemma 3 that as soon as d(x i , y i ) ≤ 2δ L , where δ L is the current best solution, we found the hyperbolicity of the graph. So, the hyperbolicity δ L of the graph is then returned in Line 5. To iterate over the pairs (v, w) such that d(v, w) ≥ d(x i , y i ), the algorithm maintains the mates of each node. Node w is a mate of v if (v, w) is a far-apart pair satisfying d(v, w) ≥ d(x i , y i ). In other words, (v, w) is a far-apart pair previously considered for some j < i such that (x j , y j ) = (v, w) or (x j , y j ) = (w, v).

To further prune the search space, we use the notions of skippable, acceptable and valuable nodes of the Borassi et al. approach that are computed by computeAccVal according to the definitions given below [START_REF] Borassi | On computing the hyperbolicity of real-world graphs[END_REF]. Algorithm 1 can be read before considering the details of this optimization. Its time complexity is in O(n 4 ) and its space complexity is in Θ(n 2 ). Indeed, it not only needs to store the list of far-apart pairs, but also the distance matrix, and the lists of mates.

Skippable, acceptable and valuable. Given a pair x, y of nodes and a lower bound δ L on hyperbolicity, Borassi et al. [START_REF] Borassi | On computing the hyperbolicity of real-world graphs[END_REF] propose a classification of the nodes to prune the nodes that cannot lead to any improvement of the lower-bound δ L known so far. For instance, a node v such that min{d(x, v), d(y, v)} ≤ δ L can be skipped, since by Lemma 3 we then have δ(x, y, v, w) ≤ δ L for any w. In Lemma 4 we summarize the conditions defining (x, y, δ L )-skippable nodes.

Lemma 4 ([12]

). A node v is (x, y, δ L )-skippable if and only if it satisfies any of the following conditions:

1. v does not belong to any far-apart pair considered before (x, y) (Lemma 5 in [START_REF] Borassi | On computing the hyperbolicity of real-world graphs[END_REF]); [START_REF] Borassi | On computing the hyperbolicity of real-world graphs[END_REF]). A node that does not satisfy any condition of Lemma 4 is defined as (x, y, δ L )-acceptable, and so it must be considered. This class is further refined by Borassi et al. [START_REF] Borassi | On computing the hyperbolicity of real-world graphs[END_REF] with the subset of c-valuable nodes, where c is any fixed node (a good choice is a node with small eccentricity or centrality) as specified in the following Lemma.

2. min{d(x, v), d(y, v)} ≤ δ L (Lemma 3); 3. 2 ecc(v) -d(x, v) -d(y, v) < 4δ L + 2 -d(x, y) (Lemma 8 in [12]); 4. ecc(v) + d(x, y) -3δ L -3 2 < max{d(x, v), d(y, v)} (Lemma 9 in
Lemma 5 ([12]). Let c be any fixed node.

A (x, y, δ L )-acceptable node v is c-valuable if 2 d(c, v) - 2δ L > d(x, v) + d(y, v) -d(x, y).
In Algorithm 1, the 4-tuples considered with far-apart pair (x i , y i ) are such that v is c-valuable, w is (x, y, δ L )-acceptable and (v, w) is a far-apart pair seen previously. Overall, the classification of the nodes is done in overall time O(n) for a given pair (x i , y i ) and lower bound δ L . Previous experiments [START_REF] Borassi | On computing the hyperbolicity of real-world graphs[END_REF] show that this classification leads to a significant reduction of the number of considered 4-tuples as well as computation time. Intuitively, the running time improvements using the classification into acceptable and valuable nodes stem from the following insights. If a node v is not valuable, then we avoid checking all mates of v. If a node v is valuable but one of its mates w is not acceptable, then we avoid the memory accesses to distances to compute δ(x, y, v, w).

Hub labeling

A main bottleneck of Algorithm 1 comes from the Θ(n 2 ) memory usage. This can be alleviated by using hub labeling [START_REF] Delling | Robust distance queries on massive networks[END_REF], a technique that allows to encode distances in a graph. The technique is also called two-hop labeling [START_REF] Cohen | Reachability and distance queries via 2-hop labels[END_REF]. It appears to give a very efficient space-time tradeoff in practice. We tried to use it as a replacement of the distance matrix in Algorithm 1. It perfectly fits in memory with all practical graphs we could test. However, we could not find a way to identify valuable nodes v for a pair x, y without computing the distances d(x, v) and d(y, v) for all v. It is then more efficient to perform BFSs from x and y. As the technique appeared as an inefficient way to extend Algorithm 1, we abandoned it.

Our algorithm

To improve upon Algorithm 1, and in particular to reduce the memory usage in practice, we do the following.

• We use the far-apart pair iterator presented in Section 3 to avoid the pre-computation and storage of the list of far-apart pairs. The use of Improvement 1 of Section 3 ensures that at the end of the visit of F d u for some u and d, it contains only the nodes w such that u is w-far, and so the pair (v, w) is far-apart.

• Let x, y be the currently considered far-apart pair. We design a function mates(F, v, d xy ) to iterate over the far-apart pairs at distance d ≥ d xy that involve v and that have previously been reported by next(F ). In the remainder we always use the function mates with the argument d xy set to d(x, y).

When d(x, y) < d ≤ diam(G), this function simply yields nodes from F d v since the order of operations of the algorithm when using Improvement 1 ensures that F d v contains only the nodes forming far-apart pairs with v.

Algorithm 2: New algorithm for computing the hyperbolicity

Input:

G = (V, E) 1 Initialize the far-apart pair iterator F 2 δ L ← lowerBoundHeuristic(G) 3 while has next(F ) do 4 (x, y) ← next(F )
// provides d(x, y)

5 if d(x, y) ≤ 2δ L then 6 return δ L 7 d x ← c x,y,δ L -prunedBFS(x) 8 d y ← c y,x,δ L -prunedBFS(y) 9 (acceptable, valuable) ← computeAccVal(x, y, δ L , G) 10 for v ∈ valuable do 11 for w ∈ mates(F, v, d(x, y)) do // provides d(v, w) 12 if w ∈ acceptable then 13 δ L ← max{δ L , δ(x, y, v, w)} 14 return δ L When d = d(x, y),
we have to ensure that a v-far node w is yielded if and only if the pair (v, w) has previously been reported by a call to next(F ) (so the pair (v, w) is far-apart). To do so, we modify the far-apart pairs iterator and its next(F ) function as follows. We use an extra hash map T (initially empty) associating to a node u the subset of u-far nodes at distance d from u that have previously been reported by next(F ). Then, when next(F ) yields a far-apart pair (x, y), we store y in T x and x in T y . This way, when d = d(x, y), function mates(F, v, d) simply has to yield nodes from T v . Finally, as soon as function next(F ) starts reporting pairs at distance d -1, we exchange hash maps T and F d (alternative implementation of Improvement 1), and proceed with a cleared hash map T .

• Instead of giving the distance matrix as input to the algorithm, we compute for each pair (x, y) the BFS distances from x and y before the call to computeAccVal(). Since the distance d(v, w) is obtained while extracting the mates of v from the far-apart pair iterator, we get all the needed distances to compute δ(x, y, v, w).

Although these repeated computations of BFS distances have no impact on the overall time complexity of the algorithm, which remains in O(n 4 ), they represent a significant computation time in practice. More precisely, the algorithm inspects all far-apart pairs (x, y) such that d(x, y) ≥ 2δ(G). Let Φ ≤ n 2 denote their number. The running time of the algorithm is then O(nm+Φm+Φ 2 ) as the initialization performs at most one BFS per node, the main loop runs two BFSs per inspected far-apart pair (x, y), and the number of (v, w) pairs inspected for (x, y) is at most the number of pairs considered so far, which is at most Φ. To reduce the impact on the overall computation time, we propose Optimizations 2 (Section 4.4.2) and 3 (Section 4.4.3) below.

See Algorithm 2 for the overall presentation of our algorithm. In the following, we present some optimizations aiming at reducing the computation time.

Optimization 1: Lower bound initialisation

The technique of acceptable and valuable nodes becomes more efficient when δ L is larger as Inequalities 3 and 4 of Lemma 4 and the inequality of Lemma 5 become stricter. For that reason, we first use the heuristic described in [START_REF] Cohen | On computing the gromov hyperbolicity[END_REF] to set an initial value to δ L . It is referenced as lowerBoundHeuristic in Algorithm 2.

Optimization 2: Cache of BFSs

We use a cache of BFSs to avoid recomputing a BFS that has recently been computed. This cache has a bounded capacity (e.g., 1000 BFSs) and stores the vectors of distances obtained through the most recent BFSs performed. This cache is used in two places: when computing distances from x and y for testing quadruples involving a pair (x, y), and when computing far-apart pairs within the function has next() of the far-apart pair iterator. Observe that even a cache of 2 BFSs is beneficial as function next(F ) reports successively all far-apart pairs at distance d involving a node x and such that x < y. The BFS from x is thus performed only once as long as we scan (x, y) pairs with same x node.

Optimization 3: Pruned BFS for searching acceptable nodes

When considering a pair (x, y), we perform BFS searches from both x and y to obtain distances from x and y and detecting both acceptable and valuable nodes. Lemma 4 allows to restrict both searches as follows.

First, we observe that if a node v satisfies ecc(v)d(x, v) < 3δ L -3 2d(x, y), then Lemma 4 applies and v is (x, y, δ L )-skippable. A (x, y, δ L )-acceptable node v must thus satisfy:

ecc(v) -d(x, v) ≥ c x,y,δ L , where c x,y,δ L = 3δ L - 3 2 -d(x, y). (1)
We then define the c x,y,δ L -pruned BFS search from x as a BFS search from x that visits only nodes satisfying Equation (1). More precisely, when visiting a node u, we enqueue only neighbors v of u that satisfy Equation [START_REF] Abu | Metric tree-like structures in real-life networks: an empirical study[END_REF]. Note that c x,y,δ L is constant given x, y and δ L . We can then safely replace the regular BFS from x by a pruned BFS as stated by the following lemma. Lemma 6. A c x,y,δ L -pruned BFS search from x visits all (x, y, δ L )-acceptable nodes.

Proof. The proof follows from two facts. First, any (x, y, δ L )-acceptable node must satisfy Equation [START_REF] Abu | Metric tree-like structures in real-life networks: an empirical study[END_REF]. Second, any node v ̸ = x satisfying Equation (1) must have some neighbor closer to x that satisfies Equation (1). This second condition proven below allows to easily prove by induction that all nodes satisfying Equation ( 1) are visited by a c x,y,δ L -pruned BFS search from x.

To prove the above second condition, consider a neighbor u of v at distance d(x, v) -1 from x. By triangle inequality, we have ecc(u

) ≥ ecc(v) -1, implying ecc(u) -d(x, u) ≥ ecc(v) -d(x, v).
As v satisfies Equation (1), so does u.

Notice that the set V x,y,δ L of nodes visited by a c x,y,δ L -pruned BFS search from x is larger than the set of (x, y, δ L )-acceptable nodes. Indeed, Conditions 1 to 3 of Lemma 4 cannot be used for the search. Furthermore, remark that the set V x,y,δ L depends on the distance d(x, y) and not on the precise node y. That is, Lemma 7. For any z such that d(x, y) = d(x, z), we have V x,z,δ L = V x,y,δ L .

The following results are direct consequences of Equation ( 1) and Lemma 7.

Corollary 3. For any z such that d(x, y) ≥ d(x, z), we have V x,z,δ L ⊆ V x,y,δ L . Corollary 4. For any δ > δ L , we have V x,y,δ ⊆ V x,y,δ L .

Lemma 7 and corollaries 3 and 4 enable the use of our pruned BFS in combination with a cache of BFSs. Indeed, during the execution of the algorithm both the considered distance d(x, y) decreases and the lower bound δ L increases. Hence, a cached c x,y,δ L -pruned BFS remains valid for future use. Hence, for Line 7 of Algorithm 2, we first check if a BFS from x is in the cache, and if so we retrieve it. Otherwise, we perform a c x,y,δ L -pruned BFS search from x and add it to the cache. We proceed similarly for y.

Observe also that to determine the sets of (x, y, δ L )-acceptable and c-valuable nodes, it suffices to consider nodes that have been visited by the c x,y,δ L -pruned BFS search, or by the c y,x,δ L -pruned BFS search if this set is smaller.

Parallelization

In this section, we discuss possibilities to parallelize our algorithm. The main issue that arises when parallelizing our algorithm is that we rely on knowing the mates, which in turn results from processing far-apart pairs. In a sequential implementation this is not an issue as a quadruple will be considered from one of the two pairs; however, if we process far-apart pairs in parallel, we might miss quadruples.

To overcome this issue, we can first build the set F d of far-apart pairs at distance d, thus precomputing all required mates. This can be done either sequentially, by iterating over the far-apart pairs at distance d, or in parallel, by first finding for each node u with eccentricity d the set of u-far nodes in parallel (we already know the u-far nodes of the nodes with larger eccentricities), and then consolidating F d to keep only far-apart pairs. Observe that each thread used to determine a set of u-far nodes requires to compute a BFS and thus uses space O(n).

Once the set F d has been consolidated, one can consider the far-apart pairs in F d in parallel, examining their associated quadruples. The only point of synchronization between threads is then the correct propagation of improvements of the lower bound δ L . Note however, that this only affects pruning and thus it is not relevant for the correctness of the algorithm but only for the performance. As we later see in the experiments (Figure 4) most of the time of the algorithm is spent visiting quadruples (and thus, to update δ L ). Hence, this two steps approach for a parallel version of our algorithm seems promising.

Experimental results

In this section we conduct experiments to test the performance of our new algorithm in comparison to the previous state-of-the-art one. We additionally test the impact that each proposed optimization has on the overall performance. To gain further insights into our main tool, the algorithm for efficiently enumerating far-apart pairs, we also conduct experiments to analyze the performance and the number of far-apart pairs that graphs exhibit.

Implementation notes

We have implemented all the algorithms in C++. Our code is publicly available (see [START_REF] Coudert | Hyperbolicity (version 1.0[END_REF]). In this section, we discuss some implementation choices.

We have implemented a cache of BFSs with bounded capacity κ that additionally maintains the age information of the data it stores. We use a counter τ , initialized to 0, that is increased by one each time a BFS that is already in the cache is accessed, or a BFS that is not in the cache is added. We associate to a cached BFS from a node x an age information a x , initialized to the value of τ at insertion time. Then, we set a x to the current value of τ each time the BFS from x is accessed. Hence, the last accessed BFS is such that a x = τ . The use of a hash map associating to a node the corresponding BFS and age information enables to decide in time O(1) if a BFS from x is in the cache, and if so to return a pointer on the corresponding data. Updating the age information is done in time O [START_REF] Abu | Metric tree-like structures in real-life networks: an empirical study[END_REF]. The insertion of a BFS in the cache takes expected time O(1) if the cache is not full, and time O(κ) as soon as it has reached its maximum capacity. Indeed, the insertion of a BFS when the cache is full requires to remove first the BFS with largest age information, that is, the one of the node x maximizing τa x , and so with smallest a x . Note that for κ much smaller than n, the time required for managing the cache is negligible with respect to the time required for a BFS. Observe that this cache will be accessed O(n 2 ) times by Algorithm 2, and more precisely at most twice the number of far-apart pairs at distance 2δ(G) or more.

Data & Hardware

To evaluate the performances of our algorithms, we use graphs from the BioGRID interaction database (BG-*) [START_REF] Oughtred | The biogrid interaction database: 2019 update[END_REF]; a protein interactions network (dip20170205) [START_REF] Salwinski | The database of interacting proteins: 2004 update[END_REF]; and graphs of the autonomous systems from the Internet (CAIDA as * and DIMES *) [START_REF]The Cooperative Association for Internet Data Analysis (CAIDA)[END_REF][START_REF] Shavitt | DIMES: Let the internet measure itself[END_REF]. We also test computer networks (Gnutella, Skitter), web graphs (notreDame), social networks (Brightkite, Epinions, Facebook, Slashdot), co-author graphs (ca-*, dblp), road networks (oregon2, FLA-t), a 3D triangular mesh (buddha), a graph from a computer game (FrozenSea), and grid-like graphs from VLSI applications (alue7065). The data is available from snap.stanford.edu, webgraph.di.unimi.it, www. dis.uniroma1.it/challenge9, graphics.stanford.edu, steinlib.zib.de, and movingai.com. Furthermore, we use synthetic inputs: grid300-10 and grid500-10 are square grids with respective sizes 301 × 301, and 501 × 501 where 10% of the edges were randomly deleted. Each graph is taken as an undirected unweighted graph and we consider only its largest (with respect to number of vertices) bi-connected component (available from [START_REF] Coudert | Hyperbolicity (version 1.0[END_REF]). See Table 1 for the characteristics of these graphs. This set of graphs has been selected as it is representative of the variety of networks (complex networks, road networks, etc.) that are considered in different fields. Most of these graphs have already been used to evaluate the performances of previously proposed exact and heuristic algorithms for hyperbolicity [START_REF] Borassi | On computing the hyperbolicity of real-world graphs[END_REF][START_REF] Cohen | Applying cliquedecomposition for computing gromov hyperbolicity[END_REF][START_REF] Cohen | On computing the gromov hyperbolicity[END_REF][START_REF] Kennedy | On the hyperbolicity of large-scale networks and its estimation[END_REF][START_REF] Abel | Quelques propriétés topologiques des graphes et applications à internet et aux réseaux[END_REF].

We used a computer equipped with two Intel Xeon Gold 6240 CPUs operating at 2.6GHz and 192G RAM to run our experiments. Note that our code uses a single thread.

Parameter choice

As there are instances which do not terminate in reasonable time or need unreasonable amounts of memory, we cap both resources at a fixed value to obtain a clear picture. More precisely, for each graph and each experiment, we kill the process as soon as it takes more than 6 hours or uses more than 192 GB of memory. We use the A symbol to indicate a killed process and, if applicable, put this symbol in the column (i.e., time or memory) that caused the process to be killed. Furthermore, for the cache described in Section 5.1, we use a size of 1000, unless mentioned otherwise.

Comparison to previous work

To evaluate the performance improvement over previous approaches, we compare to the practical state-of-the-art algorithm which is the Borassi et al. algorithm [START_REF] Borassi | On computing the hyperbolicity of real-world graphs[END_REF]. To this end, we measure the single-threaded computation time, the memory, and the best upper and lower bound found of our new algorithm as well as Borassi et al. algorithm. The results of this experiment are shown in Table 2.

There are several interesting observations that one can derive from these experiments. First, the new approach that we present in this paper needs significantly less memory. More precisely, the memory consumption is up to a factor of 28 times lower than for Borassi et al. algorithm (see the slashdot0902-d graph), only considering the graphs where both approaches stay within the limits. If we also consider the graphs where Borassi et al. algorithm runs out of memory, we use at least a factor 177 less memory (see grid300 -10). Mainly due to this significant reduction of memory consumption, we are able to compute the hyperbolicity for graphs which were previously out of reach due to excessive amounts of memory which would have been necessary. However, there are also graphs for which we can compute the hyperbolicity, which hit the timeout limit of 6 hours when using Borassi et al. algorithm. Note that in these two cases, the computation roughly takes between 4 and 5 hours, so Borassi et al. algorithm takes at least 1 to 2 hours more time to finish on these instances. In general, we note that all graphs for which Borassi et al. algorithm computes the hyperbolicity within the time and memory limits, also stay within the limits using our approach. Over all 40 benchmark instances, our new approach computes the hyperbolicity for 8 more graphs than Borassi et al. algorithm. If the hyperbolicity cannot be computed within the limits with our new approach, then this is due to running out of time instead of running out of memory. This again shows that we achieve a drastic reduction in memory consumption for computing the hyperbolicity.

While we are sometimes faster than Borassi et al. algorithm, one could also assume that sometimes it is the other way around. This is indeed the case, and while for almost all graphs we are slightly faster or slightly slower, there is one instance (namely, DIMES 201012) where we are a factor of 13 slower than Borassi et al. algorithm. We suspect that for this graph, there are many nodes from which we repeatedly perform BFSs to obtain distances. Indeed, as both algorithms follow a similar enumeration of pairs and quadruples, one could expect a slowdown with our algorithm since it needs to recompute some distances. However, the main bottleneck comes from the inspection of quadruples with regards to the 4-points condition (see Section 5.5): for each pair (x, y) inspected, we perform at most two BFSs with cost O(n + m), while the number of quadruples (x, y, v, w) inspected can be proportional to n 2 . The fact that both algorithms have comparable running times is thus explained by a marginal cost of distance computations, even with our approach that requires to recompute some of them several times. Note also that both algorithms use different orderings for visiting pairs and quadruples. This implies that one can improve its lower bound on hyperbolicity earlier than the other, greatly reducing the number of quadruples to be subsequently inspected. This is probably the main reason for the fluctuations of running time between the two algorithms.

We specifically want to highlight two grid graphs, grid300-10 and grid500-10, for which our approach only takes 10 seconds and 90 seconds, respectively, while Borassi et al. algorithm fails to compute the hyperbolicity due to exceeding the 192 GB memory limit. Note that for the former grid graph, only 1.08 GB of memory is needed for our approach and thus the memory usage is lower by at least a factor of 177. This drastic reduction of memory comes from the highly structured input: a perfect grid only has 2 far-apart pairs, the two pairs of opposing corners. Thus, the hyperbolicity computation only needs to consider these two pairs and computing the shortest paths between all pairs is hugely wasteful. Our approach includes this natural intuition to lazily compute the distances between the pairs of nodes to avoid huge running times and memory consumption in cases like these. For a better overview, we plotted a comparison of the time and memory usage of both algorithms, see Figures 2 and3. In Figure 2 we can see that the times indeed are very similar for both approaches except for the single outlier mentioned above. In Figure 3 we can see the drastic reduction in memory usage. Note again that the graphs included in this figure are only the graphs on which both algorithms terminate within the limits. The graphs where Borassi et al. algorithm exceeds the limits while our approach stays withing the limits are not shown. In particular, the memory reduction by at least a factor of 177 is not shown in this plot. Furthermore, we can see that with increasing memory consumption, also the advantage that our algorithm has over Borassi et al. algorithm with respect to memory usage becomes more and more pronounced.

Table 2 also lists the memory usage due to the BFS cache in our algorithm (last column). The intent of choosing a cache size of 1000 was to obtain better running times while having a reasonable memory consumption for big instances. This choice is questionable for some instances and we come back on that matter later on. Note however, that we obtain a much lower memory usage compared to Borassi et al. even with this choice of large cache.

Running time of parts

Additional to the comparison with the previous state-of-the-art algorithm, we also analyze the running time of the different parts of our algorithm. See Figure 4 for a visualization of how much running time as percentage of the total running time is spent in specific parts of our algorithm. First note that for each part of the algorithm, there exists and instance for which it is the bottleneck. However, we note that the initialization of our algorithm is only the bottleneck for instances that have a small total running time (see Table 2 for the running times). The other parts of our algorithm are also bottlenecks of harder instances. The most frequent bottleneck on our instances is the update of δ L , which is the part of the algorithm where we explicitly check quadruples. This suggests that the approach of enumerating far apart pairs using previous pruning techniques is at its limit and a speed-up has to come from an improved pruning of quadruples.

The next pair computation appears to be the bottleneck in instances where very few pairs are far-apart so that there are very few quadruples to check. Note also that Borassi et al. algorithm is faster on instances where the computation of distances from x and y is the bottleneck (see Table 2). This is expected as the time wasted in recomputing some distances then becomes critical.

Impact of different optimizations

We also conduct experiments to obtain insights into the impact of the different optimizations that we used in our algorithm. In particular, we want to know the impact of the lower bound initialization (see Section 4.4.1), the cache size (see Section 4.4.2), and the pruning (see Section 4.4.3). We conduct experiments with all our graphs with the normal cache size of 1000 with no heuristic and no pruning, heuristic but no pruning, pruning but no heuristic, and heuristic and pruning. To gain 10 -1 10 0 10 1 10 2 10 3 10 4

Borassi et al.

10 -1

Our Algorithm

Figure 2: Comparing the time in seconds for computing the hyperbolicity on all graphs that finish using both algorithms. The dashed line is the identity, i.e., where both algorithms would take the same time. Borassi et al.

Our Algorithm

Figure 3: Comparing the memory in MB for computing the hyperbolicity on all graphs that finish using both algorithms. The dashed line is the identity, i.e., where both algorithms would use the same amount of memory. The solid green line (2 x → 2 5.38+0.6x ) is a linear regression on the logarithmic coordinates. It suggests an exponential improvement in the memory consumption. more insight into the effect of the BFS cache, we also conduct the last experiment with a cache size of 2, which means that we only store the current two BFSs in memory. See Table 3 and Figure 5 for the results of these experiments.

The overall results of these experiments are that there is no significant benefit for the running time when using the heuristically computed initial lower bound (see Figures 5.a-b). However, the pruning has significant impact depending on the graph: sometimes it does not help much, but in several cases it does decrease the running time by a larger factor approaching 2 and up to a factor of 5 (see Figure 5.c and the DIMES 201204 in Table 3). The cache optimization can reduce the running time by a factor of 1.4, but it also sometimes slightly increases it (see Figure 5.d). The benefit of a general caching strategy thus seems marginal compared to just storing distances from x while scanning (x, y) pairs with same x. Note that this simple strategy is equivalent to our caching optimization with a cache size of c = 2. Furthermore, while the cache does not dominate the memory consumption for most instances, it still increases the memory usage significantly, see the last column of Table 2.

BFS cache size experiments

To gain further insights into how the BFS cache we use affects the running time behavior, we run experiments with different cache sizes on selected graphs, see Figure 6. We again see different behavior on different instances. While the times monotonously decrease with increasing cache size for two of the graphs, for the other two we have increase-decrease and decrease-increase patterns. One reason for this behavior might be the cache size of the processor. In particular, if the graph already fits in the CPU cache (e.g., L1), then computing a BFS is quite fast. Especially for large BFS cache sizes which might push the graph out of the CPU cache and also might reside in a higher level CPU cache, the computation can become slower. Overall, the benefit of the cache appears to be rather low in terms of running time while it can be quite costly in terms of memory usage, advocating for choosing cache size c = 2 in the end as stated in the previous paragraph.

Far-apart pairs iterator experiments

Finally, we perform experiments to only analyze the behavior of the far-apart pairs iterator. To this end, we let the far-apart pairs iterator run on all our benchmark graphs and measure the time as well as memory consumption. Additionally, we are interested in the number of far-apart pairs that the different instances have as well as how many pairs of an instance are necessary for the hyperbolicity calculation of our new algorithm. As the graphs are of different sizes, we put the number of pairs in the last two columns in relation to the total number of node pairs in the graph. The results of these experiments are shown in Table 4. For all except four graphs, the far-apart pairs iterator runs through in the given memory and time limits. This is one more tractable instance compared to hyperbolicity computation. Note, however, that there are graphs for which we can compute the hyperbolicity but the far-apart pairs iterator does not run through within the limits (e.g., dblp-2010). This is explained by the fact that to compute the hyperbolicity, we can stop at some point of iterating through the far-apart pairs and do not have to compute them all. More specifically, we show in column "hyp pairs" what percentage of all pairs are necessary to compute the hyperbolicity with our algorithm on each instance. As we only consider mates as third and fourth nodes (i.e., previous far-apart pairs) in a quadruple, counting the number of far apart pairs until our algorithm terminates indeed counts all the pairs of nodes that we consider as first and (a) cache on all graphs that finish: "cache" refers to a cache of size c = 1000 (instead of c = 2), "prune" refers to pruning, and "heur" refers to the lower-bound heuristic. The light-gray dashed lines (y = 2 * x and y = x/2) correspond to ratio 2. second tuple in the hyperbolicity computation. Observe that for the dblp-2010 graph, only around 0.27% of pairs are relevant for the hyperbolicity computation while 48.38% of the pairs are farapart. In general, most instances only have a single-digit percentage or less of relevant pairs for the hyperbolicity computation while the fraction of far-apart pairs is above 24% for all graphs except grid-like graphs. Indeed, grid-like graphs (grid300-10, grid500-10 and z-alue7065) have less than 0.1% of far-apart pairs and even less relevant pairs. The grid graphs, grid300-10 and grid500-10, have such a small number of relevant pairs, that they are rounded to 0 in the precision that we choose for the numbers in the table. These low numbers of relevant pairs give the main explanation of the drastic memory reduction that we achieve with our algorithm. A second explanation comes from the fact that we store far-apart pairs whose distance remains greater than the hyperbolicity lower-bound δ L . As we compute and store more far-apart pairs, the lower-bound increases and some previously computed far-apart pairs are cleared from memory. Figure 7 illustrates the maximum memory consumption of our hyperbolicity computation compared to the total number n 2 of pairs. We use a cache size c = 2 so that memory usage comes mainly from the storage of far-apart pairs. All in all, we can see that the memory consumption is lower than the size of the full distance matrix by roughly three orders of magnitudes for most graphs. We can also see that for many graphs, the far-apart pairs iterator is very fast, while the hyperbolicity computation takes a long time, which is explained by the fact that we might spend quadratic time per pair to compute the hyperbolicity. Considering the percentage of far-apart pairs, we see that most graphs have roughly 30% to 70% far-apart pairs. The graphs with grid structure are extreme outliers, exhibiting a very low percentage of far-apart pairs. This can be explained by the grid structure, for which -considering a grid without missing edges -only the two pairs of opposing corners of the grid are far-apart. Furthermore, the facebook combined graph has a very large percentage of far-apart pairs. This is explained by the fact that it has, by far, the highest average degree of all the graphs we consider. For such a well-connected graph, BFS trees have a very large number of leaves as shortest paths are not extended further from most nodes, which is necessary to produce such a large number of far-apart pairs.

To further gain insights into the distribution of far-apart pairs, i.e., how they are distributed over various distances, we put all the distances between far-apart pairs into a histogram for four selected graphs, see Figure 8. While the distribution looks somewhat reasonable and expected for three of the four graphs, the notreDame-d graph shows an interesting distribution with two modes. Note that anomalies with respect to coreness were already found in this graph [START_REF] Shin | Patterns and anomalies in k-cores of real-world graphs with applications[END_REF] where a special structure around a propeller-shaped subgraph was identified. We suspect that this structure connects a large fraction of pairs and could be related to this bi-modal observation. Finally, we also marked the part of the histogram that contains the far-apart pairs that are relevant for the hyperbolicity calculation. We observe that this region occurs after the peak of the histogram in all four instances that we show. In the notreDame-d instance, it even completely excludes the first, much larger mode.

Conclusion

In this work, we developed a fundamental algorithm to iterate over all far-apart pairs in nonincreasing distance. As primary application we consider the computation of graph hyperbolicity. Our new algorithm enables us to compute, for the first time, the hyperbolicity of some graphs with more than a hundred thousand nodes with non trivial structure (e.g., notreDame-d, loc-Table 4: Time and memory consumption of the far-apart pairs iterator only. We additionally show the percentage of far-apart pairs ("far pairs") and also the percentage of pairs (both with respect to the number of all node pairs in the graph) which have to be considered during the hyperbolicity computation of our algorithm ("hyp pairs"). Figure 8: Histograms of distribution of far-apart pairs for selected graphs. On the x-axis we plot the distance and on the y-axis the number of far-apart node pairs that have this distance. The green area shows the distance range for pairs that have to be evaluated in our hyperbolicity algorithm. We can see that this area always excludes the peak of the histogram in the shown examples.

Graph

gowalla edges, com-dblp.ungraph). We reduce the memory usage significantly, while not compromising on performance. Non-trivial graphs with more than five hundred thousands nodes unfortunately still remain out of reach with our method. We thus plan to investigate alternative approaches in future work in order to get closer to the million nodes barrier. Furthermore, we believe that iterating over far-apart pairs in non-increasing distance is such a fundamental task, that our work will enable faster algorithms also in other settings.

  [START_REF] Albert | Topological implications of negative curvature for biological and social networks[END_REF]-points condition on a tree.

  BG-S-Affinity Capture-RNA BG-S-Affinity Capture-Western BG-S-Biochemical Activity BG-S-Dosage Rescue BG-S-Synthetic Growth Defect BG-S-

Figure 4 :

 4 Figure 4: Percentage of time spent in different parts of our algorithm. Each row corresponds to the execution of the algorithm that was performed in Table2on the respective graph. The columns refer to the different parts of Algorithm 2. More precisely, they correspond to the following lines in the algorithm: initialize (l. 1-2), next pair (l. 3-4), (d x , d y )-BFS (l. 7-8), computeAccVal (l. 9), update (l. 10-13).

4 cache+prune+heurFigure 5 :

 45 Figure 5: Comparison of various combinations of optimizations in terms of running time (in seconds)on all graphs that finish: "cache" refers to a cache of size c = 1000 (instead of c = 2), "prune" refers to pruning, and "heur" refers to the lower-bound heuristic. The light-gray dashed lines (y = 2 * x and y = x/2) correspond to ratio 2.

Figure 6 :

 6 Figure 6: Plot for four graphs showing the running time development depending on the BFS cache size.

8

 8 memory usage in bytes (cache size c = 2)

Figure 7 :

 7 Figure 7: Memory consumption of hyperbolicity computation compared to the size n 2 of the distance matrix. We use a cache size of c = 2. The dashed line corresponds to identity while the light-gray dashed line corresponds to a factor 1/1000.

  consider these sets in decreasing distance d from v, starting with d = ecc(v), and stopping when d = c -1. For each node u ∈ L d , add each w ∈ N (u) to L d-1 for which w ̸ ∈ L d ∪ L d+1 . At the end of this procedure, we have that if d(v, u) = d, then u ∈ L d and all the sets are disjoint, i.e. L d ∩ L d ′ = ∅ for any c ≤ d, d ′ ≤ ecc(v) with d ̸ = d ′ . The claimed time bound follows immediately from the description of the algorithm.

  Observe that, since F is used to iterate over the far-apart pairs by non-increasing distances, when starting to consider far-apart pairs at distance d, we know that all pairs at distance d ′ > d have already been considered and that at initialization, we have added in F d an empty set for each node with eccentricity d. Hence, all nodes involved in a far-apart pair at distance d are known. If no node w is found, we go to the next node in F d through a call to next(F d ). When all nodes in F d have been considered, we start considering far-apart pairs at distance d -1.

	Therefore, we can define a function next(F d ) that returns either the next item (u, F d u ) in the arbitrary ordering on the items stored in the hash map F d , or Stop when all items have been
	considered. Similarly, we define a function next(F d u ) that returns either the next node in the ordering defined over the nodes in F d u , or Stop when all nodes have been considered.
	For a distance d and a node u, we repeat calls to next(F d u ) until finding a node w such that u is w-far (by testing if u is in F d w ), in which case the far-apart pair (u, w) at distance d is
	returned.

Table 1 :

 1 Some graph parameters of all the graphs that we use in our experiments. Note that for all graphs, we extracted the largest biconnected component and restrict to this subgraph in our experiments.

	Graph	#nodes #edges radius diameter mean ecc.
	BG-MV-Physical	9 851	45 558	11	14.45
	BG-S-Affinity Capture-MS	17 793	174 210	6	9.02
	BG-S-Affinity Capture-RNA	3 339	10 408	3	4.00
	BG-S-Affinity Capture-Western	9 971	44 331	10	12.60
	BG-S-Biochemical Activity	29 44	10 444	7	9.20
	BG-S-Dosage Rescue	1 521	4 143	9	12.21
	BG-S-Synthetic Growth Defect	3 013	21 341	4	5.72
	BG-S-Synthetic Lethality	2 258	12 187	4	5.45
	dip20170205	13 969	60 621	10	12.95
	CAIDA as 20000102	4 009	10 101	4	5.62
	CAIDA as 20040105	10 424	27 061	4	5.98
	CAIDA as 20050905	12 957	33 541	5	6.11
	CAIDA as 20110116	23 214	89 783	4	6.04
	CAIDA as 20120101	25 614	109 180	4	6.10
	CAIDA as 20130101	27 454	124 672	5	7.21
	CAIDA as 20131101	29 432	143 000	5	6.46
	DIMES 201012	18 764	84 851	4	5.16
	DIMES 201204	16 907	66 489	4	5.07
	p2p-Gnutella09	5 606	23 510	5	6.31
	gnutella31-d	33 812	119 127	6	7.41
	notreDame-d	134 958	833 732	18	20.99
	ca-CondMat	17 234	84 595	6	8.44
	ca-HepPh	9 025	114 046	6	7.83
	ca-HepTh	5 898	20 983	7	8.63
	com-dblp.ungraph	211 409	883 570	8	10.92
	dblp-2010	140 610	572 873	10	11.99
	email-Enron	20 416	163 257	5	6.54
	epinions1-d	36 111	365 253	5	6.36
	facebook combined	3 698	85 963	4	5.26
	loc-brightkite	33 187	188 577	6	7.95
	loc-gowalla edges	137 519	887 929	6	7.91
	slashdot0902-d	51 528	473 218	5	6.04
	oregon2 010331	7 602	27 870	4	5.89
	t.FLA-w	691 175	941 893	890	1378.52
	buddha-w	543 652 1 631 574	244	360.44
	froz-w	749 520 2 895 228	812	1130.38
	grid300-10	90 211	162 152	300	450.50
	grid500-10	250 041	449 831	500	750.49
	z-alue7065	34 040	54 835	213	319.43

Table 2 :

 2 Comparing the memory and time consumption for our algorithm and Borassi et al. algorithm. If there is a timeout, we state the best lower and upper bounds on the hyperbolicity that we obtained. The last column states the cache usage by our algorithm. For an analysis of the tradeoff between improved running time and increased memory consumption of the cache, see Sections 5.6 and 5.7. We highlight the values in the table where one algorithm evidently dominates the other algorithm with respect to running time or memory usage, respectively.

	Graph	Borassi et al. [12]			Our Algorithm
		time (s)	memory hyperb. time (s)	memory	hyperb. cache usage
	BG-MV-Physical BG-S-Affinity Capture-MS	10 429.0 A	1.03 GB -[2.0, 4.0] 4.5	6 725.0 166.33 MB A -	4.5 [2.0, 4.0]	115.83 MB -
	BG-S-Affinity Capture-RNA	0.8 157.18 MB	2.0	0.7	57.29 MB	2.0	39.22 MB
	BG-S-Affinity Capture-Western	7 827.0	1.05 GB	4.0	5 255.0 154.66 MB	4.0	116.82 MB
	BG-S-Biochemical Activity	8.4 104.90 MB	3.0	16.3	46.16 MB	3.0	34.68 MB
	BG-S-Dosage Rescue	12.0	30.54 MB	4.0	14.3	24.43 MB	4.0	17.82 MB
	BG-S-Synthetic Growth Defect	1.7 108.76 MB	2.0	2.9	43.82 MB	2.0	35.22 MB
	BG-S-Synthetic Lethality dip20170205	1.1 A	77.09 MB -[4.5, 5.0] 18 318.0 2.0 1.8	33.40 MB 339.30 MB	2.0 4.5	26.60 MB 163.96 MB
	CAIDA as 20000102	1.3 260.87 MB	2.5	1.2	56.40 MB	2.5	47.04 MB
	CAIDA as 20040105	18.4	1.90 GB	2.5	25.9 166.40 MB	2.5	122.06 MB
	CAIDA as 20050905	16.9	2.37 GB	3.0	21.8 203.72 MB	3.0	151.99 MB
	CAIDA as 20110116 CAIDA as 20120101	13 806.0 A	8.41 GB -[2.0, 2.5] 2.0 13 491.0 556.72 MB A -	2.0 [2.0, 2.5]	271.38 MB -
	CAIDA as 20130101 CAIDA as 20131101	2 961.9 A	10.09 GB -[2.0, 2.5] 16 108.0 2.5 3 535.8	1.47 GB 593.03 MB	2.5 2.5	321.09 MB 344.76 MB
	DIMES 201012	465.8	4.86 GB	2.0	6 043.0 482.41 MB	2.0	219.88 MB
	DIMES 201204	66.3	4.34 GB	2.0	56.2 304.48 MB	2.0	197.75 MB
	p2p-Gnutella09	2.9 381.64 MB	3.0	2.4	98.38 MB	3.0	65.45 MB
	gnutella31-d notreDame-d	139.5 -	13.13 GB A	3.5 -	109.5 4 514.0	1.51 GB 53.02 GB	3.5 8.0	395.77 MB 1.58 GB
	ca-CondMat	112.8	3.38 GB	3.5	172.9 281.18 MB	3.5	201.81 MB
	ca-HepPh	37.0	1.17 GB	3.0	86.3 152.86 MB	3.0	105.79 MB
	ca-HepTh com-dblp.ungraph dblp-2010	4.0 407.78 MB -A -A	4.0 --	3.0 5 449.0 6 904.0	91.54 MB 14.20 GB 7.51 GB	4.0 5.0 5.5	69.00 MB 2.47 GB 1.64 GB
	email-Enron	754.1	5.36 GB	2.5	897.7 404.38 MB	2.5	239.05 MB
	epinions1-d	696.7	18.60 GB	2.5	2 331.3 759.20 MB	2.5	422.24 MB
	facebook combined	6 919.0 248.00 MB	1.5	4 551.0 136.65 MB	1.5	44.99 MB
	loc-brightkite loc-gowalla edges	910.4 -	12.81 GB A	3.0 -	910.1 812.98 MB 14 478.0 4.04 GB	3.0 3.5	388.69 MB 1.61 GB
	slashdot0902-d	9 560.0	37.54 GB	2.5	4 718.0	1.34 GB	2.5	602.96 MB
	oregon2 010331 t.FLA-w buddha-w froz-w grid300-10 xgrid500-10	73.0 980.82 MB -A -A -A -A -A	2.0 -----	65.1 133.76 MB A -A -A -[367.5, 633.5] 2.0 [81.0, 835.5] [93.0, 221.0] 10.1 1.08 GB 280.0 99.6 2.99 GB 463.0	88.80 MB ---1.05 GB 2.92 GB
	z-alue7065	35.0	9.07 GB	138.0	33.5 431.18 MB	138.0	398.36 MB

Table 3 :

 3 Times for computing the hyperbolicity with different optimizations enabled. All entries are in seconds. The columns with "heur" use the lower bound initialization presented in Section 4.4.1, and the columns with "prune" use the pruning of Section 4.4.3. The value of c in the second row gives the size of the BFS cache presented in Section 4.4.2.

	Graph	-	heur	prune	heur & prune
		c = 1000 c = 1000 c = 1000	c = 2 c = 1000
	BG-MV-Physical BG-S-Affinity Capture-MS BG-S-Affinity Capture-RNA	7 040.0 A 0.7	7 047.0 A 0.7	6 745.0 A 0.7	6 867.0 A 0.7	6 725.0 A 0.7
	BG-S-Affinity Capture-Western	5 464.0	5 450.0	5 236.0	5 297.0	5 255.0
	BG-S-Biochemical Activity	16.5	16.5	16.2	20.2	16.3
	BG-S-Dosage Rescue	14.9	14.7	14.6	16.5	14.3
	BG-S-Synthetic Growth Defect	3.2	3.3	3.0	4.6	2.9
	BG-S-Synthetic Lethality	2.1	2.2	1.8	2.5	1.8
	dip20170205	19 088.0	19 144.0	18 492.0 18 473.0	18 318.0
	CAIDA as 20000102	1.2	1.3	1.2	1.7	1.2
	CAIDA as 20040105	39.9	40.0	25.7	35.0	25.9
	CAIDA as 20050905	35.6	36.9	21.7	23.8	21.8
	CAIDA as 20110116 CAIDA as 20120101 CAIDA as 20130101	16 079.0 A 3 292.8	16 132.0 A 3 313.8	13 700.0 13 412.0 A A 3 594.9 3 848.0	13 491.0 A 3 535.8
	CAIDA as 20131101	15 785.0	15 766.0	16 036.0 16 683.0	16 108.0
	DIMES 201012	6 144.0	6 139.0	6 032.0	6 574.0	6 043.0
	DIMES 201204	294.1	292.5	56.0	68.7	56.2
	p2p-Gnutella09	2.0	2.4	1.9	2.3	2.4
	gnutella31-d	172.2	191.2	91.1	104.7	109.5
	notreDame-d	4 317.0	4 315.0	4 472.0	4 795.0	4 514.0
	ca-CondMat	386.0	387.4	177.2	197.9	172.9
	ca-HepPh	144.7	145.0	87.8	103.4	86.3
	ca-HepTh com-dblp.ungraph dblp-2010	3.8 A 14 677.0	4.1 A 5 352.0 5 713.0 5 449.0 2.8 4.7 3.0 14 675.0 6 610.0 6 560.0 6 904.0
	email-Enron	1 121.2	1 108.8	906.5	1 184.8	897.7
	epinions1-d	2 626.6	2 607.0	2 331.4	3 563.1	2 331.3
	facebook combined	4 677.0	4 623.0	4 552.0	4 609.0	4 551.0
	loc-brightkite loc-gowalla edges slashdot0902-d	1 696.6 A 4 492.0	1 688.6 A 14 259.0 13 169.0 14 478.0 928.5 1 108.5 910.1 4 463.0 4 775.0 6 766.0 4 718.0
	oregon2 010331 t.FLA-w buddha-w froz-w	124.2 A A A	124.5 A A A	65.6 A A A	79.8 A A A	65.1 A A A
	grid300-10	10.2	10.4	10.1	8.7	10.1
	xgrid500-10	102.2	102.3	99.8	83.1	99.6
	z-alue7065	36.5	37.0	33.8	33.6	33.5

This work has been supported by the French government, through the UCA jedi Investments in the Future project managed by the National Research Agency (ANR) with the reference number ANR-15-IDEX-01, the ANR project Multimod with the reference number ANR-17-CE22-0016 and the ANR project Distancia with reference number ANR-17-CE40-0015.

A Time and space trade-offs for determining all far-apart pairs

In this section, we present algorithms offering different time and space trade-offs for the problem of computing all far-apart pairs. The space complexity considered here is the working memory, hence excluding the space needed to store the result, unless needed during computations.

A first algorithm to determine the set of far-apart pairs is to: 1) determine for each node u ∈ V the set F u of u-far nodes using BFS, and then 2) check for each node v ∈ F u if u is v-far (i.e., if u ∈ F v ). This can be done in time O(nm) and space O(n 2 ) since

Another algorithm is to execute two steps for each node u ∈ V : 1) determine the set F u of u-far nodes using BFS, and then 2) for each node v ∈ F u , check if there is w ∈ N (u) such that d(u, v) < d(w, v). The second step requires to compute distances from w and so the time complexity of this algorithm is

Observe that during the processing of each node u, this algorithm stores the set F u , the distances from u, and the distances from one neighbor w of u. Hence, the space complexity is O(n). The second method can be improved using the bit-parallel BFS proposed in [START_REF] Akiba | Fast exact shortest-path distance queries on large networks by pruned landmark labeling[END_REF]. Indeed, this algorithm computes simultaneously distances from u and b of its neighbors in time O(n + m) and space O(n), assuming that b is a constant and using bitwise operations on bit vectors of size b (typically 32 or 64). Hence, the number of BFSs to perform is reduced to u∈V ⌈|N (u)|/b⌉.

Let us now show how to modify the above algorithm to obtain an algorithm with time complexity in O(nm) and space complexity in O(n pw(G)), where pw(G) denotes the pathwidth of G [START_REF] Díaz | A survey of graph layout problems[END_REF][START_REF] Petit | Addenda to the survey of layout problems[END_REF]. The main idea is to compute the distances from u only once and to store them during as few iterations as possible. For that, let π : V → [n] be a linear ordering of the nodes (i.e., a bijective mapping) related to the pathwidth as described later and let Π(V ) be the set of all such orderings. The algorithm iterates over the nodes in the order π -1 (1), π -1 (2), . . . , π -1 (n). The distances from u are used for the processing of node u, and for the processing of each neighbor w ∈ N (u). Let w L := arg min w∈N (u) π(w) and w R := arg max w∈N (u) π(w). Hence, distances from u must be computed at iteration π(w L ) and stored until iteration π(w R ). Consequently, at iteration i, we have stored distances from all the nodes in

Now observe that the pathwidth [START_REF] Díaz | A survey of graph layout problems[END_REF][START_REF] Petit | Addenda to the survey of layout problems[END_REF] of a graph G is defined as pw(G) = min π∈Π(V ) p(G, π), where

Hence, at iteration i we store distances from at most 2p(G, π) + 1 nodes and there is an ordering ensuring to store distances from at most 2 pw(G) + 1 nodes. Consequently, the space complexity of this algorithm is in O(n pw(G)) and its time complexity is O(nm) assuming that the ordering π such that p(G, π) = pw(G) is given. However, the problem of determining an ordering π such that p(G, π) = pw(G) is NPcomplete [START_REF] Díaz | A survey of graph layout problems[END_REF][START_REF] Petit | Addenda to the survey of layout problems[END_REF]. Nonetheless, efficient heuristic algorithms have been proposed [START_REF] Coudert | Experimental evaluation of a branchand-bound algorithm for computing pathwidth and directed pathwidth[END_REF][START_REF] Castillo | A New Heuristic Algorithm for the Vertex Separation Problem[END_REF].

Finally, recall that the bandwidth of a graph G is defined as bw(G) = min π∈Π(V ) b(G, π), where b(G, π) = max uv∈E |π(u)π(v)| [START_REF] Díaz | A survey of graph layout problems[END_REF][START_REF] Petit | Addenda to the survey of layout problems[END_REF]. Therefore, distances from u are stored for at most 2b(G, π) + 1 iterations, and there is an ordering ensuring that this number of iterations is at most 2 bw(G) + 1. However, the problem of determining such an ordering is NP-hard [START_REF] Garey | Computers and Intractability: A Guide to the Theory of NP-Completeness[END_REF].

B Retrieving distances from far nodes

First note the following corollary which is a direct consequence of Lemma 1: Corollary 5. For any v ∈ V , let F v be the set of v-far nodes. The number of leaves of any shortest path tree rooted at v is at least |F v |.

Observe however that the lower bound of Corollary 5 does not imply the existence of a shortest path tree with |F v | leaves. For instance, consider the 4-cycle (u 1 , u 2 , u 3 , u 4 ). We have

but all shortest path trees rooted at u 1 have 2 leaves (either {u 2 , u 3 } or {u 3 , u 4 }).

We now show that, given the v-far nodes and a constant c, one can determine the nodes at distance at least c from v. Lemma 8. Let v ∈ V and let S be the set of all nodes u for which d(u, v) ≥ c. Given the set F v of v-far nodes and the distances from v to each of these nodes, one can compute the set S in time

Furthermore, a neighbor w ∈ N (u) with d(v, w) = d(v, u) -1 is not v-far, and if a neighbor w ′ ∈ N (u) is v-far then d(v, u) = d(v, w ′ ).

To determine all the nodes that are at distance at least c from v, it suffices to perform a reverse breadth-first search, starting from far nodes. More precisely, let {L d } d∈{c,...,ecc(v)+1} be a set family where L d is initialized with the v-far nodes at distance d from v, and let L ecc(v)+1 = ∅. Then,


A Medvedeva

F Clairet

S Hacquin

G Dif-Pradalier

Investigation of the sheared flows and staircases with a full-wave 2D synthetic diagnostic FeDoT

Introduction

Poloidal flows and avalanches are proved to participate in the tokamak plasma auto-organization and confinement through the turbulence suppression in sheared flow regions and therefore the establishment of weak transport barriers. A set of radially periodic transport barriers or E×B staircase was predicted by GYSELA simulations and indirectly observed in the Tore Supra tokamak by means of swept reflectometry [1]. However, the staircase detection remains challenging, since the shear regions size is of the order of a few millimeters, a spatial resolution inaccessible for the most part of measurement systems.

The Ultra-Fast-Swept Reflectometer (UFSR) system with probing frequencies of 50-105 GHz and the sweep time of 1 µs [START_REF] Clairet | Proc. 14th Intnl. Refl. Wksh. IRW[END_REF] has demonstrated a unique spatial resolution and a high signal-to-noise ratio during its operation on Tore Supra, ASDEX Upgrade and WEST tokamaks. These capacities have allowed to measure turbulence characteristics, such as the turbulence level, correlation and coherence lengths, at a millimeter scale covering core end edge plasma simultaneously [START_REF] Medvedeva | [END_REF]. Among other things, the UFSR measurements have provided a precise description of centimeter-size edge coherent modes (ECMs) in the sheared flow region and have shown an asymmetry in the signal frequency spectra [4]. High-frequency edge coherent modes have been detected during the experiments dedicated to the L-H transition studies on ASDEX Upgrade. They appear in different plasma confinement regimes, such as the I-phase or the H-mode with the mode frequency in the range of 40-200 kHz. UFSR frequency spectra with ECMs are highly asymmetrical in negative and positive frequencies which could be an indication of a sheared flow. This asymmetry was explained by a sawtooth-like perturbation of the cutoff [START_REF] Clairet | 45th EPS Conference on Plasma Physics[END_REF] and demonstrated this effect with synthetic signals reflected by periodic perturbations with positive and negative tilt.

The main objective of the current study is to analyze the effect of flows on the UFSR signal, frequency spectra and coherence using a reflectometry synthetic diagnostic FeDoT. If this effect is measurable, the analysis routines used in FeDoT will be integrated into an automatic tool for the staircase detection in tokamak experiments.

Reflectometry synthetic diagnostic FeDoT

The synthetic diagnostic FeDoT has been developed in particular for the UFSR X-mode data interpretation, although its modular structure is universal enough to be applied for other reflectometry systems. FeDoT includes Python routines which treat various experimental needs: from cutoff frequency visualization and analytical WKB calculation of the signal's phase to 2D full-wave simulations. As an input it uses a 2D electron density map in (R, z) coordinates which is taken either from an IMAS (ITER Integrated Modelling & Analysis Suite) shot file or a .hdf5 output of plasma simulation. The input preparation routine extracts a rectangular part of plasma poloidal cross-section of indicated size, typically of about 20 cm. In order to treat simplified models of turbulent structures, a fluid simulator has been developed, being able to study periodic structures evolving in an artificial radially localized flow. The synthetic phase can be calculated with the analytical WKB formula along the line of sight, with a 1D wave propagation code or a 2D wave propagation code FDTD-FW2D [START_REF] Medvedeva | Proc. 14th Intnl. Refl. Wksh. IRW[END_REF]. Further analysis consist of routines for frequency spectra, coherence and correlation computation as well as a standard Bottolier-Curtet profile reconstruction algorithm. The wave propagation 2D code FW2D-FDTD was developed and commissioned in 2018-2020. Based on the Yee discrete scheme of the Finite Difference Time Domain method (FDTD) [START_REF] Silva | Proc. 11th Intnl. Refl. Wksh. IRW[END_REF] the code simulates the propagation of an electromagnetic wave in 2D with X-mode polarization in the 2D computational grid. At each time step, the source beam at the antenna plane is added to the calculated field with forming a so-called soft source. The code efficiency has been studied and proved in terms of propagation, boundary absorption and cutoff reflection. The solution is easily destabilized when a strong density fluctuation is added or when the first cutoff is closer than 3 wave lengths to the wave source. The comparison between the simulated phase and an analytical solution for a linear density profile has shown a precision of 0.1 rad for a typical variation of the order of 1 000 rad along the plasma radius, which is acceptable for turbulence studies. Optimization of the numerical calculation through the Numba Python library led to a factor of hundred improvement of the simulation speed.

To mimic better the UFSR and to reduce the calculation time, we attempted to introduce a sweep of the frequency into the Yee scheme. Figure 1 shows an example of simulations for a probing wave of f 0 = 62 GHz which is stabilized at the time step 15 000 and then increasing its frequency by δf = 1, 10 or 100 MHz. The stationary regime is assumed reached when the variations of the signal amplitude at the antenna center become smaller than 1%. Then the phase and the amplitude of the reflectometer signal are extracted from the field values integrated over the antenna pattern. The largest possible frequency step is 10 MHz, as for higher frequency difference the signal stabilizes during more than 15 000 time steps and therefore the simulation takes longer than two parallel simulations with f 0 and f 0 + δf probing frequencies. The sweep used for the UFSR diagnostic (25 GHz sweep in 1µ s with 1 000 frequency points) would require 37 500 000 time steps with a change of frequency by 10 MHz every 15 000 steps. Moreover, for a fixed frequency the amplitude also starts to diverge after 25 000 time steps. Therefore, it was chosen to mimic a frequency sweep by a series of parallel simulations for 1 000 fixed frequencies. A typical grid of the order of 1 000×1 000 spatial points and 15 000 time steps is used for each frequency and turbulence map realization, which consumes about 20 CPU minutes per frequency. By the end of 2019, the FW2D-FDTD code, initially adapted to the CEA server environment, was modified to work on MARCONI HPC systems with 1024 parallel processors. The use of the MARCONI environment facilitates the coupling with gyrokinetic simulations, since GYSELA code runs on the same machines and exploits the shared data storage.

Results

The study of the sheared flow signatures in the reflectometer signal has been performed applying FeDoT to the GYSELA gyrokinetic simulation with staircases described in [1]. First we investigated the coherence length which can be reconstructed from the complex signal of the UFSR. Figure 2a depicts one of 1024 time slices of the density fluctuation map used as input for FW2D-FDTD simulation. The area marked in red corresponds to the cutoff positions of 28 probing frequencies in the range of 70-72.8 GHz, defined as the maximum of the wave electric field around the reflection region. This area contains one period of the radial electric field structure (Fig. 2c) with strong shear around 1.90 and 1.93 m, the coordinates R and Z are scaled to ASDEX Upgrade typical plasma dimensions. The coherence length is computed as a half width at half maximum of the coherence function both for the GYSELA density at Z = 0 and the cutoff position and the synthetic signal integrated over the antenna size with a Gaussian weight. Although the reflectometer coherence length (Fig. 2b in blue) is slightly lower than that of the turbulent structures, it follows the fine radial structure of the turbulent eddies size and can be used to detect the staircases. Same synthetic measurements were used to produce frequency spectra, trying to retrieve the flow signatures observed in the experiments. While a Doppler shift of the order of a kHz can exist in the UFSR signal and indicate a background flow direction, these indications are limited to the cases of the shift of the plasma in the vertical direction and consequently of the non-normal incidence of the wave. However, sheared flows might be detected through the asymmetry of the frequency spectra in negative and positive frequencies as it was shown for the ECMs. To prove this point, we performed simulations with a linear density profile and two ranges of modes (δn/n = 10%) with positive and negative tilt (Fig. 3a) moving downwards along the Z axis. The resulted frequency power spectra of the synthetic reflectometer signal (Fig. 3b) are compared for frequencies of 64 and 67 GHz probing the areas with positive and negative tilt angles. The asymmetry of the spectra of about 10-20% of the amplitude indicates the direction of the tilt and appears for the wide range of mode's frequency (1-400 kHz). The asymmetry has been quantified by introducing a spectral chirality which is the difference of a spectrum for positive and negative frequencies. A positive chirality, as for the probing wave of 64 GHz (Fig. 3b, solid blue), corresponds to a positive tilt angle counting from the Z axis. Therefore, the conventional reflectometer signal contains not only the information about the radial size of turbulent structures, but also an indication of the average tilt angle, which can be used to detect the relative flow velocity and the staircases. Quantitative analysis of the spectral chirality is a subject of future investigation. The chirality and the coherence length together might sufficiently describe the self-organized weak transport barriers and will be applied to the existing database of the UFSR measurements.

Conclusions

The synthetic diagnostic FeDoT is developed for an X-mode swept reflectometer and being updated to match ITER and community needs. The next steps will include a complete adaptation to IMAS, a simplified user interface and a possibility to use different reflectometer systems parameters through the ITER Machine Description database. An automatic tool will be created for the experimental data analysis, combining coherence and chirality. It has been shown that the coherence length periodic change created by staircase flows can be measured by a swept reflectometer. The chirality of frequency spectra indicates the direction of the tilt of turbulent structures and can be used to describe flows and staircases. The staircase dynamics will be studied with ASDEX Upgrade and WEST UFSR data in the nearest future. The identification of flows with a swept frequency reflectometer will allow to exploit fully the large database of reflectometry measurements and study the influence of the flows on the turbulent transport.

Figure. 1 :

 1 Figure. 1: Signal amplitude at the antenna center simulated with the FW2D-FDTD code for a fixed frequency of 62 GHz with a jump of frequency at t = 15000 by 1, 10 and 100 MHz

Figure. 2 :

 2 Figure. 2: The density fluctuation input for FW2D-FDTD code from GYSELA simulation for one time slice, the red area corresponds to one staircase period (a), the coherence length computed for the synthetic signal and GYSELA density around the cutoff position for frequencies 70-72.8 GHz (b), the average radial electric field of the GYSELA input (c)

Figure. 3 :

 3 Figure. 3: The density input for FW2D-FDTD code with cutoff frequencies and a contour plot of the electric field of the propagating wave at 64 GHz (a), frequency power spectra for the simulated reflectometer signal in the cases of positive (64 GHz, dashed blue) and negative (67 GHz, solid red) tilts of the turbulent structures and the chiral part for a positive tilt angle in solid blue (b)

Acknowledgments

This work has received funding from the Euratom research and training programme 2014-2018 and 2019-2020 under grant agreement No 633053. The views and opinions expressed herein do not necessarily reflect those of the European Commission. The authors wish to especially thank Garrard Conway and Ulrich Stroth for fruitful discussions on the experimental data, Robin Marcille for participating in the development of the FW2D-FDTD code and Stephane Heuraux for valuable expert comments on the numerical scheme and wave propagation.