J D Íaz

D Mitsche
email: dmitsche@unice.fr

A Dieudonné

G Perarnau

Jordi Girona

ON THE RELATION BETWEEN GRAPH DISTANCE AND EUCLIDEAN DISTANCE IN RANDOM GEOMETRIC GRAPHS

Keywords: Random geometric graphs, Graph distance, Euclidean distance, Diameter 2010 Mathematics Subject Classification: Primary 05C80 Secondary 68R10

Given any two vertices u, v of a random geometric graph G (n, r), denote by dE(u, v) their Euclidean distance and by dG(u, v) their graph distance. The problem of finding upper bounds on dG(u, v) conditional on dE(u, v) that hold asymptotically almost surely has received quite a bit of attention in the literature. In this paper, we improve the known upper bounds for values of r = ω( √ log n) (i.e. for r above the connectivity threshold). Our result also improves the best known estimates on the diameter of random geometric graphs. We also provide a lower bound on dG(u, v) conditional on dE(u, v).

Introduction

Given a positive integer n and a non-negative real function r = r(n), a random geometric graph G on n vertices and radius r is defined as follows. The vertex set V = V (G) is obtained by choosing n points independently and uniformly at random in the square S n = [-√ n/2, √ n/2] 2 (Note that, with probability 1, no point in S n is chosen more than once, and thus we assume |V | = n). For notational purposes, we identify each vertex v ∈ V with its corresponding geometric position v = (x v , y v ) ∈ S n , where x v and y v denote the usual x-and y-coordinates in S n . For every two points u, v ∈ S n , we write d E (u, v) for their Euclidean distance. Finally, the edge set E = E(G) is constructed by connecting each pair of vertices u, v ∈ V by an edge if and only if d E (u, v) ≤ r. We denote this model of random geometric graphs by G (n, r), and use the notation G ∈ G (n, r) (or often simply G (n, r)) to refer to a random outcome of this distribution. We will always assume that r ≤ √ 2n, as for r ≥ √ 2n the graph obtained is always a clique.

Random geometric graphs were first introduced in a slightly different setting by

Gilbert [START_REF] Gilbert | Random Plane Networks[END_REF] to model the communications between radio stations. Since then, several closely related variants of these graphs have been widely used as a model for wireless communication, and have also been extensively studied from a mathematical point of view. The basic reference on random geometric graphs is the monograph by Penrose [START_REF] Penrose | Random Geometric Graphs[END_REF] (see [START_REF] Walters | Random Geometric Graphs[END_REF] for a more recent survey).

The properties of G (n, r) are usually investigated from an asymptotic perspective, as n grows to infinity and r = r(n). Throughout the paper, we use the following standard notation for the asymptotic behavior of sequences of non-negative numbers a n and b n :

a n = O(b n ) if lim sup n→∞ a n /b n ≤ C < +∞; a n = Ω(b n ) if b n = O(a n ); a n = Θ(b n ) if a n = O(b n ) and a n = Ω(b n ); a n = o(b n ) if lim n→∞ a n /b n = 0; and a n = ω(b n ) if b n = o(a n ).
We also use a n b n and b n a n to denote a n = o(b n ). Finally, a sequence of events H n holds asymptotically almost surely (a.a.s.) if lim n→∞ Pr(H n ) = 1.

It is well known that r c = log n/π is a sharp threshold function for the connectivity of a random geometric graph (see e.g. [START_REF] Goel | Sharp thresholds for monotone properties in random geometric graphs[END_REF][START_REF] Penrose | The longest edge of the random minimal spanning tree[END_REF]). This means that for every ε > 0, if r ≤ (1 -ε)r c , then G (n, r) is a.a.s. disconnected, whilst if r ≥ (1 + ε)r c , then it is a.a.s.

connected.

Given a connected graph G, we define the graph distance between two vertices u and v, denoted by d G (u, v), as the number of edges on a shortest path from u to v. Observe first that any pair of vertices u and v must satisfy d G (u, v) ≥ d E (u, v)/r deterministically by the triangle inequality, since each edge of a geometric graph has length at most r. The goal of this paper is to provide upper and lower bounds that hold a.a.s. for the graph distance of two vertices in terms of their Euclidean distance and in terms of r (see Figure 1).

Related work.

This particular problem has risen quite a bit of interest in recent years. Given any two vertices u, v ∈ V , most of the work related to this problem has been devoted to study upper bounds on d G (u, v) in terms of d E (u, v) and r, that hold a.a.s. Ellis, Martin and Yan [START_REF] Ellis | Random geometric graph diameter in the unit ball[END_REF] showed that there exists some large constant K such that for every r ≥ (1 + ε)r c , G ∈ G (n, r) satisfies a.a.s. the following property: for

every u, v ∈ V such that d E (u, v) > r, d G (u, v) ≤ K • d E (u, v) r . (1.1)
Their result is stated in the unit ball random geometric graph model, but it can be easily adapted into our setting. This result was extended by Bradonjic et al. [START_REF] Bradonjic | Efficient broadcast on random geometric graphs[END_REF] for the range of r for which G (n, r) has a giant component a.a.s., under the extra condition that d E (u, v) = Ω(log 7/2 n/r 2 ). Friedrich, Sauerwald and Stauffer [START_REF] Friedrich | Diameter and broadcast time of random geometric graphs in arbitrary dimensions[END_REF] improved this last result by showing that the result holds a.a.s. for every u and v satisfying d

E (u, v) = ω(log n/r). They also proved that if r = o(r c ), a linear upper bound of d G (u, v) in terms of d E (u, v)/r is no longer possible.
In particular, a.a.s. there exist vertices u and

v with d E (u, v) ≤ 3r and d G (u, v) = Ω(log n/r 2 ).
The motivation for the study of this problem stems from the fact that these results provide upper bounds for the diameter of G ∈ G (n, r), denoted by diam(G), that hold a.a.s., and the runtime complexity of many algorithms can often be bounded from above in terms of the diameter of G. For a concrete example, we refer to the problem of broadcasting information (see [START_REF] Bradonjic | Efficient broadcast on random geometric graphs[END_REF][START_REF] Friedrich | Diameter and broadcast time of random geometric graphs in arbitrary dimensions[END_REF]).

One of the important achievements of our paper is to show that one can take the constant K for which (1.1) holds as K = 1 + o(1) a.a.s., provided that r = ω(r c ). By the aforementioned result in [START_REF] Friedrich | Diameter and broadcast time of random geometric graphs in arbitrary dimensions[END_REF], we know that the statement is false if r = o(r c ).

A similar problem has been studied by Muthukrishnan and Pandurangan [START_REF] Muthukrishnan | The bin-covering technique for thresholding random geometric graph properties[END_REF]. They The following theorem is the main result of our paper.

Theorem 1.1. Let G ∈ G (n, r) be a random geometric graph on n vertices and radius

0 < r ≤ √ 2n. A.a.s., for every pair of vertices u, v ∈ V (G) with d E (u, v) > r (as
otherwise the statement is trivial) we have:

(i) if d E (u, v) ≥ max 12(log n) 3/2 /r, 21r log n , then d G (u, v) ≥ d E (u, v) r 1 + 1 2 (rd E (u, v)) 2/3 ; (ii) if r ≥ 224 √ log n, then d G (u, v) ≤ d E (u, v) r 1 + γr -4/3 , where γ = γ(u, v) = max 1358 3r log n r + d E (u, v) 2/3 , 4 • 10 6 log 2 n r 8/3 , 30000 2/3 .
In order to prove (i), we first observe that all the short paths between two points must lie in a certain rectangle. Then we show that, by restricting the construction of the path on that rectangle, no very short path exists. For the proof of (ii) we proceed similarly. We restrict our problem to finding a path contained in a narrow strip. In this case, we show that a relatively short path can be constructed. We believe that the ideas in the proof can be easily extended to show the analogous result for d-dimensional random geometric graphs for all fixed d ≥ 2.

Remark 1.1.

(1) Note that the condition d E (u, v) ≥ max{12(log n) (2) Similarly, the constant 224 in the condition r ≥ 224 √ log n of (ii) (as well as those in the definition of γ) is not optimized either, and could be made slightly smaller.

However, our method as is cannot be extended all the way down to r ≥ log n/π = r c .

(3) Moreover, the error term in part

(i) is 2(rd E (u, v)) 2/3 -1 = O(1/ log n) = o(1).
(4) Finally, the error term in (ii) is

γr -4/3 = Θ max log n r 2 + rd E (u, v) 2/3 , √ log n r 4 , r -4/3 , which is o(1) iff r = ω( √ log n) = ω(r c ).
Hence, for r = ω(r c ), statement (ii) implies that a.a.s.

d G (u, v) ≤ (1 + o(1)) d E (u, v) r ,
thus improving the result in [START_REF] Ellis | Random geometric graph diameter in the unit ball[END_REF].

diam(G) ≤ √ 2n r 1 + γr -4/3 , (1.3)
where

γ = Θ r log n √ n 2/3 + log 2 n r 8/3 + 1 .
(In fact, (1.3) holds for all r ≥ (1 + ε)r c as a consequence of (1.2)). From straightforward computations, one can check that (1.3) improves (1.2) provided that r = Ω log 5/8 n (log log n) 1/8 . On the other hand, for the lower bound on the diameter, observe the following: for any function ω growing arbitrarily slowly with n, we can a.a.s. find two vertices u and v, each at distance at most ω from one corner (opposite from each other) of the square S n . For such two vertices, we trivially (and deterministically) have

diam(G) ≥ d G (u, v) ≥ √ 2n 1 -Θ(ω/ √ n) r (1.4)
Assuming that log 3 n/n r √ n/ log n, our bound from Theorem 1.1 applied to these vertices gives that a.a.s.

diam(G) ≥ d G (u, v) ≥ √ 2n 1 -Θ(ω/ √ n) r 1 + Θ(r -2/3 n -1/3 ) ≥ √ 2n 1 -Θ(ω/ √ n) + Θ(r -2/3 n -1/3 ) r -1 . (1.5)
Assuming the additional constraint r n 1/10 , we have that r -2/3 n -1/3 r/ √ n, and also r -2/3 n -1/3 ω/ √ n (for ω tending to infinity sufficiently slowly). In this case, our bound in (1.5) improves upon the trivial lower bound (1.4), and can be written as

diam(G) ≥ √ 2n r 1 + Θ(r -2/3 n -1/3 ) . (1.6)
Note that this is a.a.s. still valid if we drop the constraint r log 3 n/n, since for r = O log 3 n/n the random geometric graph is a.a.s. disconnected (and has infinite diameter). Hence, by (1.3) and (1.6), we obtain the following corollary:

Corollary 1. Let G ∈ G (n, r
) be a random geometric graph on n vertices and radius

0 < r ≤ √ 2n. A.a.s. we have: (i) if r ≥ (1 + ε)r c , then diam(G) ≤ √ 2n r 1 + γr -4/3 , where γ = Θ r log n √ n 2/3 + log 2 n r 8/3 + 1 . (ii) if r n 1/10 , then diam(G) ≥ √ 2n r 1 + Θ(r -2/3 n -1/3 ) .

Proof of Theorem 1.1

In order to simplify the proof of Theorem 1.1 we will make use of a technique known as de-Poissonization, which has many applications in geometric probability (see for ex. [START_REF] Penrose | Random Geometric Graphs[END_REF] for a detailed account of the subject). Here we sketch it.

Consider the following related model of a random geometric graph G with two

distinguished vertices u, v. The vertex set of G is V = V (G) = {u, v} ∪ V ,
d E (u 1 , u 2 ) ≤ r. We denote this new model by G u,v (n, r).
The main advantage of defining V = V \ {u, v} as a Poisson point process is motivated by the following two properties: the number of points of V that lie in any region A ⊆ S n of area a has a Poisson distribution with mean a; and the number of points of V in disjoint regions of S n are independently distributed. Moreover,

conditional on N = n-2, the distribution of G u,v (n, r) is the one of G (n, r). Therefore, since Pr(N = n -2) = Θ(1/ √ n), any event holding in G u,v (n, r) with probability at least 1-o(f n ) must hold in G (n, r) with probability at least 1-o(f n √ n).
We make use of this property throughout the article, and do all the analysis for a graph G ∈ G u,v (n, r)

or related models of Poisson point processes.

We will need the following concentration inequality for the sum of independently and identically distributed exponential random variables. For the sake of completeness we provide the proof here.

Lemma 2.1. Let X 1 , . . . , X m be independent exponential random variables of parameter λ > 0 (i.e. expectation 1/λ) and let

X = X 1 + • • • + X m .
Then, for every ε > 0 we have

Pr X ≥ (1 + ε)E(X) ≤ 1 + ε e ε m ,
and for any 0 < ε < 1 we have

Pr X ≤ (1 -ε)E(X) ≤ (1 -ε)e ε m ≤ e -ε 2 m/2 .
Proof. Let us prove the bound for the upper tail. The bound for the lower tail is proved in a similar way and its proof is omitted.

We have EX = mEX 1 = m/λ. By Markov's inequality, for every 0 < β < λ and every ε > 0 we have

Pr(X ≥ (1 + ε)EX) = Pr(e βX ≥ e β(1+ε)m/λ ) ≤ E(e βXi ) e β(1+ε)m/λ = (ϕ X1 (β)) m e -β(1+ε)m/λ ,
where ϕ X1 (β) = E(e βX1 ) = λ λ-β is the moment-generating function of an exponentially distributed random variable with parameter λ. Thus,

Pr(X ≥ (1 + ε)EX) ≤ λ λ -β m e -β(1+ε)m/λ Now we set β = ε 1+ε λ to obtain Pr(X ≥ (1 + ε)EX) ≤ 1 + ε e ε m .

Proof of statement (i)

In this subsection we prove the lower bound in Theorem 1.1. For every t ≥ 0, we introduce the following model of infinite random geometric graphs. The vertex set is constructed by adding two vertices u = (0, 0) and v = (t, 0) to a homogeneous Poisson point process of intensity 1 in the infinite plane R 2 . We denote this new model by

G ∞ u,v (r, t).
The main task in the sequel is to show that, for any t ≥ max 12(log n) 3/2 /r, 21r log n , the lower bound in part (i) of Theorem 1.1 holds with probability at least 1 -o(n -5/2 ) in G ∞ u,v (r, t), for the distinguished vertices u = (0, 0) and v = (t, 0). Combining this with an appropriate de-Poissonization argument will allow us to conclude the desired result for G (n, r).

Our next lemma shows that short paths connecting u and v are contained in small strips. The lemma is stated in the more general context of a deterministic geometric graph G = (V, E) of radius r, where the vertex set V is an arbitrary subset of points in R 2 (containing u and v) and edges connect (as usual) every pair of vertices at Euclidean distance at most r. For a given r > 0, for every k ∈ N and α > 0, consider the rectangle

R(k, α) = -α 2 kr , kr × [-α, α] .
Lemma 2.2. Given any r, t, α > 0 and any k ∈ N satisfying t ≥ kr -2α 2 kr , let G be a geometric graph of radius r in R 2 , and suppose that u = (0, 0) and v = (t, 0) are two vertices of G. Then all paths of length at most k from u to v are contained in R(k, α).

Proof. If there is no path of length at most k from u to v, the statement of the lemma is trivially true. Thus, we suppose that P = (u = z 0 , z 1 , . . . , z = v) is a path of length ≤ k, where z i = (x i , y i ) for every 0 ≤ i ≤ . Also, note that it suffices to prove the lemma for α satisfying t = kr -2α 2 kr , since this trivially implies the statement for larger α. In particular, we have α 2 < (kr) 2 /2.

Write x + = max{x i : 0 ≤ i ≤ } and x -= min{x i : 0 ≤ i ≤ }. It is clear that
x + ≤ r ≤ kr since every edge has length at most r. Let 0 ≤ j ≤ be such that

x j = x -and observe that x -≤ x 0 = 0. Then kr ≥ d E (u, z j ) + d E (z j , v) ≥ -x -+ (t -x -) ≥ kr -2 x -+ α 2 kr ,
and we obtain x -≥ -α 2 /(kr).

Now write y + = max{y i : 0 ≤ i ≤ } and y -= min{y i : 0 ≤ i ≤ }. We will show that y + ≤ α and that y -≥ -α. For every 0 ≤ i ≤ , we have

kr ≥ d E (u, z i ) + d E (z i , v) = x 2 i + y 2 i + (t -x i ) 2 + y 2 i ≥ t 2 + 4y 2 i ,
where we used the fact that the left-hand side of the last inequality is minimized at

x i = t/2. Using that t ≥ kr -2α 2 kr , we obtain (kr) 2 ≥ t 2 + 4y 2 i ≥ kr - 2α 2 kr 2 + 4y 2 i .
Thus, for every 0 ≤ i ≤ , we have |y i | ≤ α 1 -α 2 /(kr) 2 ≤ α, and so in particular

-α ≤ y -≤ y + ≤ α.
Using the bounds on x + , x -, y + and y -, we conclude that P is contained in R(k, α).

Proposition 2.1. For every t > r, let G ∈ G ∞ u,v (r, t) be a random geometric graph on R 2 with additional vertices u = (0, 0) and v = (t, 0). Then, for every 0 < δ < 2 -1/3 , we have that

Pr d G (u, v) ≤ t r 1 + δ (tr) 2/3 ≤ (1 + o(1))t r exp -δ/2(tr) 2/3 + exp -(1 - √ 2δ 3 -o(1)) 2 t 2r .
(2.1) Proof. We first set

k = t r 1 + δ (tr) 2/3
and α = δ 2

k 3 r 2 t 1/3 . Observe that since t > r, we have k ≥ 1. Let A 1 the event that d G (u, v) ≤ k; that is,
there exists a path P in G ∞ u,v (r, t) from u to v of length at most k. Our goal is to show that the probability of A 1 is small. Let x 1 be the largest x-coordinate of the vertices inside the rectangle R 1 = [0, r] × [-α, α] (possibly x 1 = 0 if u is the only vertex in R 1 ). Define the random variable a 1 = r -x 1 . We proceed similarly for every 2 ≤ i ≤ k. We define x i as follows: if

R i = (x i-1 + a i-1 , x i-1 + r] × [-α, α]
is non-empty, let x i be the largest x-coordinate of the vertices inside R i ; otherwise, set x i = x i-1 + a i-1 (see Figure 2). Define then

also a i = x i-1 + r -x i . Claim: If A 1 holds, then t ≤ x k .
Proof of the claim. Suppose that A 1 holds, let P = (u = z 0 , z 1 , . . . , z = v) be one such path and for every 0 ≤ i ≤ and let xi be the x-coordinate of z i . We will prove by induction on i that we have xi ≤ x i . In particular, this implies t = x ≤ x ≤ x k , and proves the claim.

Observe that t ≥ kr Therefore, our choice of x 1 and the fact that

1 + δ (tr) 2/3 ≥ kr 1 - δ (tr) 2/3 = kr - 2α 2 kr . ( 2
z 1 ∈ R × [-α, α] imply that x1 ≤ x 1 . So,
the statement holds for i = 1. Now we inductively assume that xi-1 ≤ x i-1 . We must have xi ≤ xi-1 +r (since z i-1 and z i are adjacent vertices), and therefore xi ≤ x i-1 +r.

Similarly as before, since z i ∈ R × [-α, α] and by the choice of x i , we conclude that xi ≤ x i , as desired. This completes the proof of the claim.

Thus, it suffices to show that x k ≥ t with very small probability. We first study the random variables a i . Define a 0 = 0. By the choice of x i , for every 1

≤ i ≤ k we have that 0 ≤ a i ≤ r -a i-1 . Recall that R i = (x i-1 + a i-1 , x i-1 + r] × [-α, α]. Since G ∈ G ∞ u,v (r, t
), the number of vertices from V inside a region of R 2 is a Poisson random variable with mean equal to the area of that region. So, for every 2 ≤ i ≤ k we have

Pr(a i ≥ β) =      Pr((x i-1 + r -β, x i-1 + r] × [-α, α] empty) = e -2αβ if 0 ≤ β ≤ r -a i-1 0 if β > r -a i-1 . (2.3)
Thus, a i is stochastically dominated by an exponentially distributed random variable ãi of parameter 2α. We assume that a i and ãi are coupled together in the same probability space, so that a i = min{ã i , r -a i-1 } ≤ ãi .

Moreover, since the regions R 1 , R 2 , . . . , R k that define the random variables a i are disjoint, the joint distribution of a 1 , a 2 , . . . , a k is stochastically dominated by the joint distribution of ã1 , ã2 , . . . , ãk , that is, the distribution of k i.i.d. exponentially distributed random variables of parameter 2α. Expanding recursively from the relations x i = x i-1 + r -a i and x 1 = r -a 1 , we get

x k = k i=1 (r -a i ) = kr -a.
Let us consider the event A 2 defined by ãi ≤ r/2 for all 1 ≤ i ≤ k. Since we aim to bound the probability that x k is large (or equivalently, that a is small), we cannot use the fact that the joint distribution of the a i 's is stochastically dominated by the ones of ãi 's. Nevertheless, note that conditional on A 2 , we have a i = ãi for all 1 ≤ i ≤ k;

if a i-1 ≤ r/2 then a i ≤ r/2 ≤ r -a i-1
, and from (2.3), a i = ãi . In other words, for

every β ≥ 0 Pr(a ≤ β, A 2 ) = Pr(ã ≤ β, A 2 ) .
Since each ãi is exponentially distributed with parameter 2α and stochastically dominates a i , we can bound the probability that A 2 does not occur:

Pr(A 2 ) ≤ k i=1 Pr (a i > r/2) ≤ k i=1
Pr (ã i > r/2) = ke -αr .

(2.4)

Therefore, using the bound on t given in (2.2), we have

Pr(x k ≥ t) ≤ Pr(A 2 ) + Pr(x k ≥ t, A 2 ) ≤ ke -αr + Pr(kr -a > t, A 2 ) ≤ ke -αr + Pr a ≤ 2α 2 kr , A 2 = ke -αr + Pr ã ≤ 2α 2 kr , A 2 ≤ ke -αr + Pr ã ≤ 2α 2 kr . (2.5)
Thus it remains to give a good upper bound on the lower tail of ã. Notice that

E(ã) = k i=1 E(ã i ) = k 2α .
We use the definition of k and α, as well as Lemma 2.1 with

ε = 1 - √ 2δ 3/2 -o(1) to show Pr ã ≤ 2α 2 kr = Pr ã ≤ 4α 3 k 2 r • E(ã) ≤ Pr ã ≤ √ 2δ 3/2 + o(1) E(ã) ≤ e -ε 2 k/2 . (2.6)
Finally, we use (2.5), (2.6) and the definition of k, α and ε to obtain

Pr (x k ≥ t) ≤ (1 + o(1))t r exp -δ/2(tr) 2/3 + exp -(1 - √ 2δ 3 -o(1)) 2 t 2r .
Since the event A 1 implies x k ≥ t, Pr(A 1 ) ≤ Pr(x k ≥ t) and the proposition follows.

Proposition 2.2. Let G ∞ u,v (r, t) be a random geometric graph in R 2 with two additional distinguished vertices u = (0, 0) and v = (t, 0) such that

t = d E (u, v) ≥ max 12(log n) 3/2 /r, 21r log n .
(2.7)

Then we have

d G (u, v) ≤ t r 1 + 1 2(rt) 2/3
, with probability at most o(n -5/2 ).

Proof. Set δ = 1/2. Since t ≥ 12(log n)

Pr d G (u, v) ≤ t r 1 + 1 2(rt) 2/3 = o(n -5/2 ) .
The same conclusion in Proposition 2.1 must be true (for t ≤ √ 2n) if we restrict

G ∞ u,v (

Proof of statement (ii)

In this subsection we complete the proof of Theorem 1.1. To derive the bound on the upper tail on the graph distance between u, v ∈ V , we first assume that u = (0, 0) and v = (t, 0) (for some 0 < t ≤ √ 2n), and analize G ∞ u,v (r, t) restricted to a suitable rectangle. Our goal is to find a path P from u to v inside of that rectangle that gives an appropriate upper bound on d G (u, v). Then, we will use similar ideas to those at the end of Subsection 2.1 to derive the desired conclusion about G (n, r).

For every measurable set S ⊆ R 2 containing u and v, let G S,u,v (r, t) denote the random geometric graph obtained as the subgraph of G ∞ u,v (r, t) induced by the vertices contained in S. Observe that G S,u,v (r, t) can also be constructed by taking as the vertex set a Poisson point process of intensity 1 in S, adding the vertices u = (0, 0) and v = (t, 0), and connecting any two vertices by an edge if they are at Euclidean distance at most r.

For every 0 < α ≤ r, we define the rectangle

S(t, α) = [0, t] × [0, α].
(The precise value of α will be specified later; it will be different to the one given in the previous subsection.) Given α and r, we write ρ = r -α 2 r . Then, for every point z = (x z , y z ) ∈ S, we define the rectangle

S z = S z (α) := [x z , x z + ρ] × [0, α] .
We need the following auxiliary lemma. Lemma 2.3. Let t > 0 and 0 < α ≤ r. Then, for every pair of points z ∈ S(t, α) and z ∈ S z (α), we have d E (z, z ) ≤ r (see Figure 3).

Proof. It is enough to show that the upper-left corner z 1 = (x z , α) and the bottomright corner z 2 = (x z + ρ, 0) of S z (α) satisfy d E (z 1 , z 2 ) ≤ r. Then all the points inside S z (α) lie at distance at most r, and in particular d E (z, z ) ≤ r.

We have

(d E (z 1 , z 2 )) 2 = ρ 2 + α 2 = r 2 -α 2 1 -(α/r) 2 ≤ r 2 ,
and the lemma follows. Our next task is to bound the graph distance between u and v in G S(t,α),u,v (r, t) by finding a path of length at most t r 1 + δr -4/3 from u to v, for some δ that will be made precise in the following proposition. Proposition 2.3. Let F > 0 and J > 3(F + 1) be constants and define g(x) =

x -log(1 + x). For every J ≤ δ ≤ F r 4/3 , there exists an α such that the following holds: fix t ≥ 0 and consider G S,u,v (r, t) to be a random geometric graph with u = (0, 0), v = (t, 0) in the rectangle S = S(t, α). Then we have

Pr d G (u, v) > t r 1 + δr -4/3 ≤ n exp - (F + 1)δ 1/2 r 4/3 2J 3/2 +exp -g (δ/J) 3/2 t r .
Proof. Let us first define some parameters that will be useful in our analysis. Set C = J -3/2 and let B be an arbitrary positive constant satisfying

B 2 + 2C/B < 1/(F + 1). (2.8)
Some elementary analysis shows that such B must exist. In fact, the equation B 2 + 2C/B = 1/(F + 1) has exactly two positive solutions B 1 and B 2 for any 0 < C =

J -3/2 < 1 (3(F +1)) 3/2
, and any 0

< B 1 < B < B 2 < 1/ √ F + 1 satisfies (2.8).
Fix some δ with J ≤ δ ≤ F r 4/3 , and set

α = Bδ 1/2 r 1/3 . (2.9)
In order to use Lemma 2.3, let us first show that α ≤ r. Since δ ≤ F r 4/3 by hypothesis of the proposition, we have

α ≤ (B √ F )r, (2.10)
and B √ F < 1, since B < 1/ √ F + 1. Moreover, we have

ρ = r -α 2 /r ≥ (1 -B 2 F )r. (2.11)
Let us consider the integer k = t r (1 + δr -4/3 ) and let A 1 be the event that d G (u, v) ≤ k; that is, there exists a path P = (u = z 0 , z 1 , . . . , z , v) in G S,u,v (r, t) from u to v of length at most k. Such a path will only use vertices inside S = S(t, α), but due to some technical considerations in the argument, we extend the Poisson point process of our probability space to the semi-infinite strip S(α) = [0, ∞) × [0, α]. Our goal is to show that the probability of A 1 is large.

As we did in the proof of Proposition 2.1, now we define random variables x i and a i for every i ≥ 1. Set x 0 = 0 and a 0 = 0. For each i ≥ 1, consider the rectangle

R i = R i (α) := (x i-1 + ρ/2, x i-1 + ρ] × [0, α]. If R i
contains at least a vertex, let z i be the vertex with largest x-coordinate inside R i . In such a case, define x i to be the x-coordinate of z i and a i = x i-1 + ρ -x i . Otherwise, we stop the process.

Let τ = min{i ≥ 1 : R i contains no points} be the stopping time of the process.

Claim: Conditional on τ ≥ k, if x k-1 + ρ ≥ t, then A 1 holds.
Proof of the claim. Assume that τ ≥ k and that x k-1 + ρ ≥ t. Observe that for every i < k, we have 0 ≤ a i ≤ ρ/2. Moreover, by construction of the process, for every 1 ≤ i < k, we have z i ∈ R i ⊆ S zi-1 and, since α ≤ r, Lemma 2.3 implies that z i is adjacent to z i-1 . Thus, the vertices z 0 , z 1 , . . . , z k-1 form a path. In particular,

x 1 ≥ ρ/2.
(2.12) Since x k-1 + ρ ≥ t, we know that there exists a value ≤ k -1 such that x -1 + ρ ≥ t and also x -1 ≤ t, and thus, by Lemma 2.3, z and v are connected by an edge. The path P = (u = z 0 , z 1 . . . , z , v) has length + 1 ≤ k, connects u and v and is fully contained in S. Therefore, A 1 is satisfied, which completes the proof of the claim.

It suffices to show that we have with high probability τ ≥ k, and that conditional on it, with high probability x k-1 + ρ ≥ t.

For every 0 ≤ i < τ , let

A (i)
2 be the event that a j ≤ ρ/2 for every 1 ≤ j ≤ i and let A 2 = A (k-1) 2 be the event that τ ≥ k. Conditional on A 2 , we have that the regions R 1 , . . . , R k-1 are disjoint. Hence, we deduce that conditional on A 2 , the joint distribution of a 1 , . . . , a k-1 is the same as the joint distribution of ã1 , . . . , ãk-1 , with ã1 , . . . , ãk-1 being k -1 independent exponentially distributed random variables with parameter α. In particular, conditional only on A (i-1) 2

, we also have that a i is stochastically dominated by ãi , and hence,

Pr(A 2 ) = k-1 i=1 Pr a i ≥ ρ/2 | A (i-1) 2 ≤ k-1 i=1 Pr(ã i ≥ ρ/2) . Since αρ/2 ≥ (1 -B 2 F )αr/2 = (1 -B 2 F )Bδ 1/2 r 4/3 /2, we have Pr(ã i ≥ ρ/2) = e -αρ/2 ≤ e -(1-B 2 F )Bδ 1/2 r 4/3 /2 . and that Pr(τ < k) = Pr(A 2 ) ≤ ne -(1-B 2 F )Bδ 1/2 r 4/3 /2 .
(2.13) Also, if we let a = (2.14)

Observe that now E(ã) = k-1 α . Let A 3 be the event that ã ≤ (1 + Cδ 3/2 ) k-1 α . We first show that A 3 implies the event {kρ -ã > t}. Conditional on A 3 , using the definition of α, the fact that δ -3/2 ≤ C and that δ ≤ F r 4/3 , we have

kρ -ã > kρ - 1 + Cδ 3/2 (k -1) α ≥ kr 1 - α 2 r 2 - (1 + Cδ 3/2 ) αr ≥ t(1 + δr -4/3 ) 1 -δr -4/3 B 2 + (δ -3/2 + C) B ≥ t(1 + δr -4/3 ) 1 -δr -4/3 B 2 + 2C B = t 1 + δr -4/3 1 -δr -4/3 + 1 B 2 + 2C B ≥ t 1 + δr -4/3 1 -(F + 1) B 2 + 2C B > t. (2.15)
Now, we can use (2.14) and the upper-tail bound in Lemma 2.1 to prove

Pr(A 3 ) = Pr ã ≥ (1 + Cδ 3/2 ) k -1 α ≤ e -g((δ/J) 3/2 )(k-1
) ≤ e -g((δ/J) 3/2 )( t/r -1) .

(2.16)

By expanding the definition of x k-1 , we can write x k-1 = (k -1)ρ -a. Thus, using (2.13), (2.14), (2.15) and (2.16) we obtain 1) .

Pr({τ < k} ∪ {x k-1 + ρ ≤ t}) = Pr(A 2 ) + Pr(x k-1 + ρ ≤ t, A 2 ) ≤ ne -(1-B 2 F )Bδ 1/2 r 4/3 /2 + Pr(kρ -a ≤ t, A 2 ) ≤ ne -(1-B 2 F )Bδ 1/2 r 4/3 /2 + Pr(kρ -ã ≤ t) ≤ ne -(1-B 2 F )Bδ 1/2 r 4/3 /2 + Pr(A 3 ) + Pr(kρ -ã ≤ t, A 3 ) ≤ ne -(1-B 2 F )Bδ 1/2 r 4/3 /2 + e -g((δ/J) 3/2 )( t/r -
(2.17) Moreover, by the properties of B and the definition of C, we have (1

-B 2 F )B > (1 -B 2 (F + 1))B > 2C(F + 1) = 2(F + 1)J -3/2 . Thus, Pr(x k-1 + ρ ≤ t) ≤ ne - (F +1)δ 1/2 r 4/3 J 3/2
+ e -g((δ/J) 3/2 )( t/r -1) , concluding the proof of the proposition.

Remark 2.1. Observe the trade-off between δ and the success probability in the proof of Proposition 2.3: for a given value of δ, we set α = Θ( √ δr 1/3 ). That is, for a given radius r, the smaller δ, the smaller α. Proposition 2.3 computes the probability that a path using vertices only within a strip of width α can be found. Clearly, the smaller δ, the straighter a path has to be, and the smaller the rectangle in which we have to find a path has to be, therefore making also α smaller. On the other hand, for smaller α, the probability of indeed finding a path in such a small strip also gets smaller.

Proposition 2.4. Given t > 0 and the vertices u = (0, 0) and v = (t, 0), let γ = γ(t) be defined as in the statement of Theorem where we have used that g(x) ≥ x/2 if x ≥ 3, and that C = J -3/2 . Therefore, Pr(d G (u, v) > k) ≤ n -5.0009/2 + n -5.004/2 = o(n -5/2 ).

(u, v) > k) = o(n -5/2 ). On the one hand, δ ≥ D log 2 n r 8/3 implies (1 -B 2 F )Bδ 1/2 r 4/3 2 -log n ≥ (1 -B 2 F )BD
Corollary 2.1. Statement (ii) in Theorem 1.1 is true.

Proof. We will use an argument similar to that at the end of Subsection 2.1 to relate the models G S(t,α),u,v (r, t) and G u,v (n, r). However, such endeavour entails extra difficulties. Given two vertices u, v ∈ S n = [-√ n/2, √ n/2] 2 at Euclidean distance t > 0, there are exactly two isometries that map them to (0, 0) and (t, 0), denoted by π + and π -. Unfortunately, the preimage of the rectangle S(t, α) under such isometries may not be entirely contained in the square S n . In order to overcome this obstacle, we just need to show that the internal vertices of the path from (0, 0) to (t, 0) that we built in the proof of Proposition 2.3 are contained in a smaller rectangle whose preimage under either π + or π -is contained in S n . This will be enough for us to conclude the existence (with sufficiently high probability) of a path in G u,v (n, r) between u and v of the desired length.

Recall the definition of α given in (2.9). Observe that from (2.12) together with (2.18),

x 1 ≥ ρ/2 > 4α/3 with probability at least 1-o(n -5/2 ). In particular, this event implies

Figure 1 :

 1 Figure 1: Graph distance vs. Euclidean distance between two points u and v in V

  proposed a new technique to study several problems on random geometric graphsthe so called Bin-Covering technique -which tries to cover the endpoints of a path by bins. They consider, among others, the problem of determining D G (u, v), which is the length of the shortest Euclidean path connecting u and v. Recently, Mehrabian and Wormald[START_REF] Merhabian | On the stretch factor of randomly embedded random graphs[END_REF] studied a similar problem to the one in[START_REF] Muthukrishnan | The bin-covering technique for thresholding random geometric graph properties[END_REF]. They deploy n points uniformly in [0, 1] 2 , and connect any pair of points with probability p = p(n), independently of their distance. In this model, they determine the ratio of D G (u, v) and d E (u, v) as a function of p.

  where the position of u and v is selected independently and uniformly at random in S n , and where V is a set obtained from a homogeneous Poisson point process of intensity 1 in the square S n of area n. Observe that V consists of N points in the square S n chosen independently and uniformly at random, where N is a Poisson random variable of mean n. Exactly as we did for the model G (n, r), we connect u 1 , u 2 ∈ V by an edge if and only if

Figure 2 :

 2 Figure 2: Example of some values of xi and their corresponding ai (i ∈ {1, 2, 3}).

. 2 )

 2 Thus, we can use Lemma 2.2 to show that the path P is contained in the strip R × [-α, α]. Moreover, we must have x1 -x 0 ≤ r (since u = z 0 and z 1 are adjacent vertices).

  r, t) to any arbitrary square S n of area n containing u = (0, 0) and v = (t, 0) (i.e. we consider the subgraph induced by the vertices lying inside of that square), since the graph distance between u and v can only increase when doing so. Moreover, by rotating and mapping an appropriate square S n to S n = [-√ n/2, √ n/2] 2 , we conclude that statement (i) in Theorem 1.1 holds in G u,v (n, r) with probability 1 -o(n -5/2 ). Hence, in view of the de-Poissonization argument described in the beginning of Section 2, this same property holds in G (n, r) with probability 1 -o(n -2 ), for a given pair of vertices u, v. The statement follows by taking a union bound over all at most n 2 pairs of vertices.

Figure 3 :

 3 Figure 3: The rectangle Sz

k- 1 i=1 a i and ã = k- 1 i=1

 11 ãi , conditional on A 2 (or in other words, on τ ≥ k), by the same argument, for every β ≥ 0, we have Pr(a ≥ β, A 2 ) = Pr(ã ≥ β, A 2 ) .

that z 1

 1 is outside of the square [0, 1.01α] × [0, α]. If z (the last internal vertex of the path P found) is outside [t -1.01α, t] × [0, α], we obtain a path connecting u and v of length +1 ≤ k with all its internal vertices in R:= [1.01α, t-1.01α]×[0, α]. Otherwise, suppose that z lies in [t -1.01α, t] × [0, α].Then, also with probability 1 -o(n -5/2 ), we can find some point ẑ in [t -1.01α -r/2, t -1.01α) × [0, α]: indeed, since ρ ≤ r, the region in which we want ẑ is bigger than the regions S i in the proof of Proposition 2.3 and Proposition 2.4. We now can use Lemma 2.3 to show that z can be replace by ẑ

Figure 4 :

 4 Figure 4: The black area represents the copy of R contained in Sn possible value of β. Notice that the minimal angle β under our assumptions on y u and y v is obtained for u = (-√ n/2, -√ n/2 + α) and v = ( √ n/2, √ n/2 -α), and thus β ≥ arctan √ n-2α √ n , or equivalently tan ( π 2 -β) ≤ √ n √ n-2α . In this case,

  Therefore, it suffices to apply Proposition 2.3 with δ = γ . It is straightforward to check that the restrictions (2.8) and J > 3(F + 1), required in Proposition 2.3, hold. We also need to show that J ≤ γ ≤ F r 4/3 . Notice that D log 2 n

	r ≥ 224	√	log n. Then, we have
				d G (u, v) >	t r	1 + γr -4/3	,
			γ = max E	log n t/r -1	2/3	, D	log 2 n r 8/3 , 3 2/3 J .
	≤ F r 4/3 , since ≤ F r 4/3 , since t/r -1 ≥ 1, and r 8/3 since r 2 ≥ E 3/2 log n, which follows from our assumption of r ≥ 224 r ≥ 224 √ log n ≥ D 1/4 √ log n; also E log n 2/3 t/r -1 √ log n; and finally
	3 2/3 J ≤ F r 4/3 since r = Ω( √	log n). Moreover, γ ≥ 3 2/3 J ≥ J.
	Note that this choice of constants combined with (2.10) and (2.11) implies
			α ≤	0.01 2.02 √	2	r ≤ r/3	and	ρ ≥ 8r/9 ≥ 8α/3.	(2.18)

1.1. Let G S,u,v (r, t) be a random geometric graph in the rectangle S = S(t, α), with additional vertices u and v. Suppose that with probability at most o(n -5/2 ).

Proof. First, observe that, if t ≤ r, then d G (u, v) = 1 with probability 1, and the statement of the proposition holds trivially. Thus, we assume henceforth that t > r.

Set B = 0.01/(2.02 √ 2), C = 10 -4 , F = 1, D = 4 • 10 6 , E = 1358 and J = 10 8/3 . Set Note that γ ≤ γ for γ as given in Theorem 1.1: indeed, the second and the third term are equal, and for the first term, for t > r, we have that (3/(1 + t/r)) > (1/( t/r -1))

holds.

We will now apply (2.17) in the proof of Proposition 2.3 with this given δ, in order to show that Pr(d G

  On the other hand, δ ≥ 3 2/3 J and δ ≥ E

						2	1/2 log n	-log n
		>	7.0009 2	log n -log n =	5.0009 2	log n .
				log n t/r -1	2/3	imply
	g (δ/J) 3/2 ( t/r -1) >	(δ/J) 3/2 2	( t/r -1) ≥	1 2	CE 3/2 log n >	5.004 2	log n,

Acknowledgements

The authors are grateful to the anonymous referees for their constructive comments.

in P . Observe that ẑ is connected to v, since 1.01α + r/2 ≤ ρ, and also ẑ is connected to z -1 , since its x-coordinate x satisfies |x -x -1 | ≤ max{ρ, r/2 -ρ/2} ≤ ρ. Thus, we can replace z with ẑ , and obtain a new path connecting u and v of length +1 ≤ k with all its internal vertices in R. We will show that either π + (R) or π -(R) is always contained in S n . We first introduce some definitions.

Consider two points u = (x u , y u ) and v = (x v , y v ) in S n . By symmetry we may assume that x u < x v and y u ≤ y v . Let β be the angle of the vector uv with respect to the horizontal axis. Again by symmetry, we may consider β ∈ [0, π/4].

We consider now two rectangles of dimensions α × t placed on each side of the segment uv. Let R + be the rectangle to the left of uv (that is, R + = π + (R)), and let R -be the rectangle to the right of uv (also, R -= π -(R)). We will show that at least one of these rectangles contains a copy of R fully contained in S n . This choice will determine which of the isometries, π + or π -, map R inside S n .

Notice that the intersection of R + and R -with each of the halfplanes x ≤ x u ,

x ≥ x v , y ≤ y u and y ≥ y v gives 4 triangles. We call them T + u , T - v , T - u and T + v respectively. All these triangles are right-angled, and denote by t + u , t - v , t - u and t + v the side of the corresponding triangle that it is parallel to the segment uv. Notice that

Call a triangle T * w , with w ∈ {u, v} and * ∈ {+, -}, safe if |t * w | ≤ 1.01α. Note that if T + u and T + v are safe or fully contained in the square, then R + contains the desired rectangle R, and analogously for R -.

Since we assumed that β ≤ π/4, we have |t

it is clear that either R + or R -contain the desired copy of R. Thus, we may assume that β > 0.

We can also assume that both u and v are on the boundary of S n , as otherwise we extend the line segment uv to the boundary of the square, and the original rectangles are contained in the new ones.

Recall that T + u and T - v are safe. If y v < √ n/2 -α, then T + v is completely contained in the square, and hence R + satisfies the conditions. Similarly, if y u > -√ n/2 + α, R -satisfies the conditions. Thus, assume that y v ≥ √ n/2 -α and y u ≤ -√ n/2 + α.

We will show that R -contains the desired copy of R: as before, T - v is safe, so it remains to consider T - u . We have |t - u | = α tan ( π 2 -β). For 0 < β ≤ π/4, tan ( π 2 -β) is decreasing in β, and it therefore suffices to show that T - u is safe for the smallest


Franck Corbier

Accelerate the Development of Certified Software for Train Control & Monitoring Systems

Keywords: Transportation, Railway, TCMS -

The rail transportation industry is highly dynamic and driven by the need to meet with mandatory safety criteria. Modern railway transportation systems are ever more sophisticated and rely heavily on embedded systems and communication networks. Train builders and their suppliers have to develop and deliver these increasingly sophisticated embedded systems while meeting increasingly stringent quality, safety and certification constraints.

Driven by the EN-50128 standard 'Software For Railway Control And Protection Systems' [3], which proposes needed methods to be used in order to provide software for control systems, train builders have applied processes and adopted tools that help to develop and validate software in compliance with safety constraints. This paper will illustrate how model based system engineering approaches and frameworks provide full business process support to:  Validate the functional requirements of the systems (requirement management);  Design control and safety functions at the train, vehicle and equipment levels (functional design and software architecture);  Design and develop the software applications of electronic controllers;  Integrate, validate and qualify electronic systems in a progressive integration process;  Manage all changes during the lifecycle of the product (change management but also option and variant management);  Demonstrate the compliance with the required safety integrated level,  Train and educate drivers, rail network operators and maintenance technicians. This kind of development and validation process and the associated tools are now successfully used by many of today's train manufacturers and equipment providers.

Control engineers in railway commonly used development tools providing the IEC61131-3 languages 'Programmable Controllers -Programming Languages' [1] that are deployed in industrial automation for a long time. In this paper, we will give an overview of the languages properties and their integration in the development of software for railway.

Since the latest version of EN-50128 introduced in 2011, train control system providers are facing increased constraints on the development processes but also on the verification and validation activities. The certification of a railway control system demands ever more documentations and demonstrations that the software complies with the safety rules and criteria. All these additional tasks increase the development and validation load but are nevertheless legitimate since they are designed to ensure the safety of the product. As the railway operators ask for SIL2 product, this paper will present a new technologies helping railway OEMs and suppliers to increase the Safety Integrated Level of control software from SSIL0 to SSIL2 and to master the cost of the product.

State of The Art

In railway industry, the functions (i.e. doors, traction, brake, lighting, fire detection, air conditioning) done with conventional control system (hardware based contactors, relays, circuit breakers, etc.) are increasingly being managed by intelligent control units. Those controllers are coordinated or synchronized by a centralized and networked control system named TCMS (Train Control and Monitoring System). Depending of the safety level of the functions, the software is developed using various tools:  softPLC providing IEC61131-3 languages  Model-Based System Engineering frameworks (providing previous languages or data flow, state machine, etc.) and dedicated code generators;  Hand coded in C or ADA for legacy reasons. Some of these tools are coming from automation and control engineering (mainly for SSIL0-2 software applications) and from the avionic industry (EN50128-SSIL4 is roughly the same asDO178C-DALA) even if a product qualified for a standard is not automatically qualified for another standard).

With the growth of embedded controllers and networks in the rolling stock, train manufacturers and their sub-system suppliers need development methods and tools that are able to support the design and help stakeholder to validate requirements and functional specifications as early as possible in the project. In fact, major issues in developing these systems is the development process itself (focusing on the delivery of safety documentation) and the decreasing time and budget available for the test activities including physical tests. Like in other industries, the railway control engineers are looking for virtual test, progressive integration and qualification of the electronic control units and electro mechanical equipment.

In parallel, rolling stock and sub-system providers need to assess the overall safety of the systems, for both passenger and equipment, and to demonstrate compliance with the relevant railway rules & standards (EN-50126, EN-50128 and EN-50129 for hardware, software and product). Traditionally, the design process for embedded electronic boards is built upon a number of elementary practices and follows a V Cycle as show in figure 1. The EN50128 recommend the implementation of this kind of development process including verifications and validation tests.

Figure 1: V-Cycle development process with documentation delivery

Moving from a paper based development cycle…

Functional and design specifications are exclusively supported by paper based communications. The project goes slowly because the verification of textual documents requires many review iterations and exchanges between business units, central office, sub-contractors and customers. Too much time is spent in "meeting" clarification and explanation making the sentences of the specifications clear, complete and understandable by all stakeholders involved in the project. Today, too much effort is spent on "paperwork" rather than on systems design, integration and physical tests. The target is not to reduce the number and the content of paper (mandatory) but to ease this production.

The deployment of a V-Cycle methodology also introduces problem. As we can see in figure 1, the left side of the V-Cycle supports creation, definition, and refinement activities while the right side covers the testing activities. It means that if coding errors can be found immediately during the programming and software module testing phases, the specification errors can be detected only at the late phases of the V-Cycle during integration and validation phases. Consequently, projects frequently miss their delivery deadlines and are difficult to keep within budget. This legacy process raises important managerial issues that need to be addressed.

… to a Model Based System Engineering Approach

The first proposal was to move from a paper based approach to a model based system engineering approach that provides executable specifications in order to speed up the development of the control software and meet project deadlines. This process was first launched with Electronic Control Units (ECUs) for Train Control & Monitoring System projects and is today deployed on other sub systems like brake, doors or passenger information systems. For this purpose, Dassault Systèmes provides a development framework that enables the achievement of new productivity related objectives.

Model a concept, a control function, a software module

This development framework, named ControlBuild, provides a MIL (Model In the Loop) approach allowing railway engineers to specify and design a model of each control function prior to any implementation. This approach is made possible through the use of known, open and standardized languages.

In the railway industry, control engineers develop the embedded software using the languages defined within the IEC61131-3 standard for many years (see Figure 2). Both control engineers (including software and electrical engineers) and maintenance engineers (end users) daily use Sequential Function Chart, Structured Text, Ladder and Function Bloc Diagram for software development and modifications i.e. applying system corrections and improvements. Currently, dozens of PLC and CPU manufacturers provide coding tools based on the IEC61131-3 standard, taking the main place in the coding activities. Therefore, Control Engineers naturally improve their efficiency focusing more on design activity than on the coding phase thanks to these languages. The only requisites were that the tool has to provide:  Code generators to ease the production of the source code of the software (in compliance with the required coding rules and properties);  Tools and advanced features to verify, to test and to assess the models (static metrics, quality level of the model, assertions, automatic tests, code coverage, dead code identification, traceability with the requirements, etc.) as presented in the figure 3.

The functional design can be manually validated using scope, button, virtual panels and synoptic. Tests of the specification and design models can also be automatically executed using test procedures and functional tests features. This incremental validation approach enables time to be saved on software design and integration testing activities. The detection and the correction of an error during the system validation phases is close to 100 times more expensive than when the tests are made on the model during the design activity.

Map the functional model on the hardware architecture

The model of the system is independent of the hardware (manufacturer, brand, operating system, etc.) and the architecture (i.e. networks and protocols). It means that we have designed the model of the software which can be deployed on one controller or distributed on many controllers. At this step of the design phase, we have to map parts (called POU -Program Organisation Unit) of the model on a defined hardware topology of the train system or sub-system. The use of the FBD (Function Block Diagram) language helps the software architect to allocate the POU on the virtual controllers (refer to figure 4).

Figure 4: System Mapping -Models to Controllers

Certification of a Software Application

For railway control and protection systems, a certificate ensures that a product is compliant with the safety standards. To get a certificate, the product must be certifiable. It means that the product has been developed in such a way that it can be certified. The EN 50126 standard addresses system issues on the widest scale [START_REF]EN50126, Railway applications -The specification and demonstration of Reliability, Availability, Maintainability and Safety (RAMS)[END_REF] and the EN 50129 standards addresses the approval process for individual systems [START_REF]Railway applications -Communication, signalling and processing systems -Safety related electronic systems for signalling[END_REF]. Railway products embed software applications that have to be certified regarding the EN50128 standard. As in other industries, this railway standard define the notion of safety level to be achieved (SILgeneric standard, DAL for avionic, ASIL for automotive and SSIL for Railway) [START_REF] Boulanger | Formal Methods -Industrial Use from Model to Code[END_REF]. In avionic, the safety level determines the goal to achieve while in railway it is rather a means to implement.

SSIL -Software Safety Integrated Level

Even if the standard define 5 SSIL from 0 to 4 (4 = high level of criticality), the railway industry commonly proposes to use only 3 safety levels (technics to achieve the level 1 are closer to level 2, and level 3 is also not so far from level 4):  SSIL0: comfort, passenger information, energy …  SSIL2: traction/brake, doors  SSIL4: automatic train control To get the certification that a software application is compliant with a given SSIL, the software development team has to demonstrate to the evaluation engineers / certification authorities that:  A well-defined software development process is used (quality assurance like ISO9001) including verification and validation tests,  The software application covers every and only the customer requirements (using traceability requirement tools like Reqtify),  The software application complies with the safety coding rules (like MISRA rules -Motor Industry Software Reliability Association),  The software application is covered by test procedures,  All associated document specification and test reports are available.

Be careful, SSIL0 doesn't mean 'no safety'; it is just the first level of safety and requires the same quality assurance process than for other SSIL. The safety case is just lighter.

Challenges and issues

Train operators are demanding even more SIL2 products for their trams, metros, inter-cities and high speed trains. It means that the Software Safety Integrated Level of the software applications of these new products have also to be compliant with SSIL2. In parallel, the latest version of the EN-50128 standard introduces more stringent activities and constraints on the development of software systems. Event for SSIL0 software application, number of companies improved their development process replacing manual programming (mainly C code) by model based design and automatic code generators.

Development tools supporting the IEC61131-3 standard are the most used in the railway community of control engineers for SSIL0-2 software applications. A well-known tool used in the avionic industry, is rather used for the development of SSIL4 software for signalling -automatic train control systemsthat concern less than 5% of the embedded software in a train).

An IEC61131-3 tool usually provides two kinds of code generators:  Generation of a sequence of instructions to be executed by an interpreter. The certification is difficult to achieve because the execution machine / runtime must be certified by the tool provider. Some softPLC providers got a SIL certificate (linked to the IEC61508 generic standard) covering a limited set of instructions (AND, OR, NOT, etc.) and reserved for the execution of Boolean expressions;  Generation of a source code to be compiled and executed on a real time platform.

To reduce the dependency of the source code to the compiler, the EN50128 standard propose to define a subset of the language allowing a controllable execution (no jump, limited use of pointer, no dynamic memory access, no if without else, a default for switch case, only one return inside a function…).

ERTS 2 -2016 There are many issues automatically solved because the IEC61131-3 languages are already limited (no pointer in structured text, no direct access to the operating system, etc.). Limited doesn't mean simple but sufficient to provide the code for the controller of the traction/brake systems or for the other converters, doors, air conditioning, pantograph, etc. More than that, the modelling tools provide features and metrics allowing the developer to measure the complexity of the structure of the models (depth of the hierarchy, number of interfaces, number of lines, number of blocks, and number of sub-loops for example). As the IEC61131-3 standard includes textual language (ST: Structured Text), the model described in ST must also be compliant with rules such as those defined by the MISRA group. However, at the end of the development process, the generated software application has to be certified like a manually programmed certifiable software application.

Certified SSIL2 Software generated by a Certified Code Generator

The aim of a software code generator is to help developers to remove the coding activity and to ease some verification and safety demonstration. It is why most of the development tools provide an automatic translator that takes into account a systems model and delivers source code for implementation on electronic control units. At this step, the software engineers have to provide deliverables and demonstrations for the evaluation of the static and dynamic criteria defined in the safety standard:  First, from a quality point of view, the generated software has to be readable, testable, verifiable and, maintainable. This means that coding rules have to be assessed as if the software were written by hand (depth of the function calls, number of overlapping 'IF' statements, number of inputs, outputs & parameters, the justification of use of pointer or global variables, dead code, out of boundary protection, defensive code, etc.).  Another activity is to demonstrate that the proposed test cases totally cover the test requirement specification document, and then the customer requirements by using requirement traceability tools.  Then, the software engineers have to demonstrate that the generated software has the same behaviour as the IEC61131-3 models if validation activities were already done at the design level. If not, all unit tests have to be executed and reported as if the software code were manually written.

The aim of a certified code generator is to insure that the generated code satisfy a certain level of confidence regarding the accuracy of the outputs relative to inputs. A certified code generator helps developers to remove unit tests and number of manual demonstrations.

Prerequisite: use of a well-defined input format

Since the first release of the IEC61131-3 programming standard, users want to be able to exchange their programs, libraries and projects between development environments. in fact, there are more than one hundred tools providing the IEC61131-3 textual and graphical languages but saving the description using their own specific format. Although this was not the intent of the standard itself, it was a task that the independent organization PLCopen [START_REF]PLCopen for efficiency in automation[END_REF] committed itself to. IEC61131-3 is focused on the software development environment.

This resulted in a workgroup named TC6 for XML. This committee defined an open interface between all different kinds of IEC61131-3 software tools. It provides the ability to transfer the information that is graphically presented on the windows of one coding tool to other coding tool. This format supports textual description (i.e. for structured text, SFC activity and receptivity) but also graphical information, like position and size of the blocks, connections (including route) between the objects of the language.

The SSIL2 Certified Code Generator

With aims of sustainability, we decided that our SSIL2 Certified Code Generator would take XML files defined by the PLCopen organization as the input and provide SSIL2 C code as the output. As in avionic, it is mandatory that the SIL level of the code generator has the same SIL level as the code it produce.

The figure 6 shows that the SSIL2 Certified Code Generator is composed of 4 modules (light grey):  The PLCopen XML comparator: this module verifies that the XML input file is compliant with the version of the tool.  The IEC61131-3 library: the C functions associated to each IEC61131-3 function block have been certified by static analysis and unit tests.  The XML to C code generator.  The report generator of the code generation execution These 4 certifiable software modules have been developed following a SSIL2 development, verification and validation process. The evaluation process has been executed by the certification authority named CERTIFER [START_REF] Certifer | [END_REF] which deliver the 8270/0157 certificate: SILCoder version 3.00 meets the 551L2 requirements of the standard EN50128:2011. The certificate is associated to an annex which gives the list of the definition documents (requirements, plans, specifications, design, tests and reports) with the related version used for the evaluation, the results of the evaluation and also the perimeter of the certification.

Some tangibles and intangibles results

Contribution of Model based design and model based testing methodologies

A model-based development process must be used closely to the requirement management process. Modelling frameworks allows stakeholders to simulate the requirements, the functional specification, the design specification and then the software and the hardware (engineers from the automotive fully play with Model In the Loop, Software In the Loop, Hardware In the Loop and today Driver or Human In the Loop). The customer requirements and needs are immediately understood by the system provider. Functional specification can be dynamically explained to the customer, same for design, software and hardware. For example, OEM provides us return of investment showing that only 3 reviews compare to 10 before are necessary to get the design acceptance from the customer.

Executable specifications allow easy verification and validation with both internal and external non experts groups. End users like drivers or maintenance engineers are now involved in the specification reviews and design reviews rather than at the end of the project. Verification and validation at every stage of the V-Cycle significantly decreases system failure risks by finding errors at the earliest.

Value of a Certified Code Generator

The SSIL2 Certified Code Generator delivers SSIL2 Certified C Code. As the software application is not just certifiable but automatically certified, software engineers can remove some demonstrations activities from their development and verification processes. The certificate gives the commitment and ensures that the software application is compliant with the SSIL2 readability and testability criteria and the recommended programming quality rules.

The software to embed in a real time controller is composed by two parts (see figure 6): the software application that is specific to one project and the "basic software" which makes the link between the software application and the operating system or directly the hardware if the target has not an OS.

As the specific application part is automatically certified using the Certified SSIL Code Generator, the product provider has only to certify the software dependent part of the generic target using his corporate SSIL2 Development and Verification & Validation process. In fact the certification of the basic software has to be done one time independently of the projects by the team managing the execution platform as a product: Each train maker defines and certifies a platform (hardware + basic software) that will be used during 5-10 years for multiple railway projects).

The final product is based on certified generic basic software (from OEM) and a certified specific software application (from our SILCoder). Certificates of the platform and of the code generator replace lot of tests and demonstration and save time and money. The Control System provider can focus directly on the software integration activities and software/hardware integration activities.

What we test is what we embed

How to demonstrate that the behaviour of the model has the same behaviour as the software? This big issue is solved using our IEC61131-3 development framework (ControlBuild) and our certified code generator (SILCoder). Each IEC61131-3 model is saved using the PLCopen format (native format for ControlBuild) and a unique C code is automatically generated from the XML PLCopen description. The C code of each model is generated one time and is linked with the right target platforms for simulation, code generation, debug, and deployment or test bench objectives as we can show in the figure 8. It means that we reuse the same generated C code (from the models) to provide software applications but also real time simulators (actuator, sensors and virtual controllers) for the HIL platforms (Hardware-inthe-Loop, Iron Bird in avionic). The test engineers focus on the virtual world making number of verifications. Like with Virtual Iron Bird, it is easier to execute test procedures and insert failures on the virtual environment (actuators/sensor) of the virtual product under test. Objective of virtual bench is to find and remove the last errors while the validation on HIL is to make the final demonstrations and to provide the deliverables for certification authorities.

Figure 2 :

 2 Figure 2: Simple models using IEC61131-3 languages

Figure 3 :

 3 Figure 3: Test functionalities and features

Figure 5 :

 5 Figure 5: Simple example of C code generation from LD or FBD model

Figure 6 :

 6 Figure 6: Simple example of C code generation from SFC

Figure 6 :

 6 Figure 6: Certified Code Generator

Figure 7 :

 7 Figure 7: Improvement of the V-Cycle using Model based Design and Certified Code Generator

Figure 8 :

 8 Figure 8: What is simulated is what is embedded

Functional Requirements Quality Metrics Code Coverage Time Stamping Events Automatic Test Vectors Test Procedures Boundaries Assertions Contracts Test execution Reports Replay Scenario

ERTS 2 -2016

Value for the railway industry

The IEC61131-3 languages, mainly used by the software engineers, will continue to take an ever more important place in development and delivery of software applications for railway. As the Certified Code Generator delivers the right level of safety for the control software, the gap to change the level of the certification of the software applications from SSIL0 to SSIL2 is significantly reduced. But for sure, it also means that the system provider must still upgrade his development and verification & validation process and integrate the use of modelling and code generation tools.

As in avionic with tools like SCADE®, the value of such a model based design and testing approach using ControlBuild coupled with a Certified Code Generator is promising in terms of productivity benefits and organizational effectiveness for the train makers. For the whole process, the time required to design and implement SIL2 Certified Train Control & Monitoring Systems is now halved with the ability to capitalize and reuse systems knowledge and certified assets on others projects through the standardization of train functions (and their certified products and software).

With the other IEC61131-3 tools providers, we are engaged in the official standardization of the XML PLCopen format. The name of this future standard will be IEC61131-10. This will open the use of the Certified SSIL2 Code Generator for all the developer's community. More and more products will be compliant with the SIL2 level for the safety of passengers, operators and systems.

Terminology


Anthony Bonato
email: abonato@ryerson.ca

David F Gleich
email: dgleich@purdue.edu.

Myunghwan Kim

Dieter Mitsche

Pawel Pralat

Amanda Tian

Stephen J Young

Pawe L Pra Lat

Dimensionality of Social Networks Using Motifs and Eigenvalues

come   L'archive ouverte pluridisciplinaire

Introduction

Empirical studies of on-line social networks as undirected graphs suggest these graphs have several intrinsic properties: highly skewed or even power-law degree distributions [START_REF] Barabási | Emergence of scaling in random networks[END_REF][START_REF] Faloutsos | On power-law relationships of the internet topology[END_REF], large local clustering [START_REF] Watts | Collective dynamics of "small-world" networks[END_REF], constant [START_REF] Watts | Collective dynamics of "small-world" networks[END_REF] or even shrinking diameter with network size [START_REF] Leskovec | Graph evolution: Densification and shrinking diameters[END_REF], densification [START_REF] Leskovec | Graph evolution: Densification and shrinking diameters[END_REF], and localized information flow bottlenecks [START_REF] Estrada | Spectral scaling and good expansion properties in complex networks[END_REF][START_REF] Leskovec | Community structure in large networks: Natural cluster sizes and the absence of large well-defined clusters[END_REF]. These are challenging properties to capture in concise models of social network connections and growth [START_REF] Kim | Multiplicative attribute graph model of real-world networks[END_REF][START_REF] Kolda | A scalable generative graph model with community structure[END_REF][START_REF] Gleich | Moment based estimation of stochastic Kronecker graph parameters[END_REF], and many models only possess them in certain parameter regimes. One model that captures these properties asymptotically is the geometric protean model (GEO-P) [START_REF] Bonato | Geometric protean graphs[END_REF]. It differs from other network models [START_REF] Barabási | Emergence of scaling in random networks[END_REF][START_REF] Leskovec | Graph evolution: Densification and shrinking diameters[END_REF][START_REF] Kumar | Stochastic models for the web graph[END_REF][START_REF] Leskovec | Kronecker graphs: An approach to modeling networks[END_REF] because all links in geometric protean networks arise based on an underlying metric space. This metric space mirrors a construction in the social sciences called Blau space [START_REF] Mcpherson | Evolution on a dancing landscape: Organizations and networks in dynamic Blau space[END_REF]. In Blau space, agents in the social network correspond to points in a metric space, and the relative position of nodes follows the principle of homophily [START_REF] Mcpherson | Birds of a feather: Homophily in social networks[END_REF]: nodes with similar socio-demographics are closer together in the space.

In order to accurately capture the observed properties of social networks-in particular, constant or shrinking diameters-the dimension of the underlying metric space in the GEO-P model must grow logarithmically with the number of nodes. The logarithmically scaled dimension is a property that occurs frequently with network models that incorporate geometry, such as in multiplicative attribute graphs [START_REF] Kim | Multiplicative attribute graph model of real-world networks[END_REF] and random Apollonian networks [START_REF] Zhang | High-dimensional Apollonian networks[END_REF]. Because of its prevalence in these models, the logarithmic relationship between the dimension of the metric space and the number of nodes has been called the logarithmic dimension hypothesis [START_REF] Bonato | Geometric protean graphs[END_REF]. This hypothesis generalizes previous analysis which shows that individuals in a social network can be identified with relatively little information. For instance, Sweeney found that 87% of the U.S. population had reported attributes that likely made them unique using only zip code, gender and date of birth, and concluded that few attributes were needed to uniquely identity a person in the U.S. population [START_REF] Sweeney | Uniqueness of simple demographics in the U.S. population[END_REF]. Here, we find evidence of the log-dimension property in real world social networks.

We emphasize that the present paper is the first study that we are aware of which attempts to quantify the dimensionality of social networks and Blau space. While we do not claim to prove conclusively the logarithmic dimension hypothesis for such networks, our experiments, such as those of [START_REF] Sweeney | Uniqueness of simple demographics in the U.S. population[END_REF], suggest a much smaller dimension in contrast to the overall size of the networks. Interestingly, speculation on the low dimensionality of social networks arose independently from theoretical analysis of mathematical models of social networks in [START_REF] Kim | Multiplicative attribute graph model of real-world networks[END_REF][START_REF] Bonato | Geometric protean graphs[END_REF][START_REF] Zhang | High-dimensional Apollonian networks[END_REF].

Our findings provide evidence for dimensional properties underlying social networks that have a number of potential applications in future studies. First, the dimensional properties could be used for further classification and characterization of different types of networks. Second, many NP-hard optimization problems related to graph properties and community detection are polynomial time solvable in a dimensional metric space, and thus, our findings suggest new techniques to explore for understanding why we may expect to solve these problems in social networks. Finally, if techniques to find these dimensions emerge, we should be able to create powerful new methods to harness the insight they offer into the network structure.

MGEO-P

The particular network model we study is a simple variation on the GEO-P model that we name the memoryless geometric protean model (MGEO-P), since it enables us to approximate a GEO-P network without using a costly sampling procedure. Both GEO-P and the MGEO-P model depends on five parameters described in Table 1.

The nodes and edges of the network arise from the following process. Initially the network is empty. At each of n steps, a new node v arrives and is assigned both a random position q v in R m within the unit-hypercube [0, 1] m and a random rank r v from those unused ranks remaining in the set 1 to n. The influence radius of any node is computed based on the formula:

I(r) = 1 2 r -α n -β 1/m .
With probability p, the node v forms an undirected connection to any preexisting node u where D(v, u) ≤ I(r v ), where the distances are computed with respect to the following metric:

D(v, u) = min { q v -q u -z ∞ : z ∈ {-1, 0, 1} m } ,
and where • ∞ is the infinity-norm. We note that this implies that the geometric space is symmetric in any point as the metric "wraps" around like on a torus. The volume of space influenced by the node is r -α v n -β . Then the next node arrives and repeats the process until all n nodes have been placed. In the MGEO-P model, the process ends here, whereas in the GEO-P model, the network then removes the least-recently added node, and inserts a new node following the same procedure. This iterative replacement process continues until it reaches it reaches a random point.

Figure 1 illustrates two features of the model. First, after a few steps, only a few nodes exist and even a large influence region will only produce a few links. Second, when the number of steps approaches n, a large influence region will produce many links. The idea behind the model is a simple abstraction of the growth of an on-line social network. When the network is first growing (few steps), even influential members will only know a few other members who have also joined. But after the network has been around for a while (many steps), influential members will begin with many friends.

We formally prove that the MGEO-P model has the following properties. Let α ∈ (0, 1), β ∈ (0, 1α), p ∈ (0, 1] and m be positive integer. The following statements hold with probability tending to 1 as n tends to ∞. See the MGEO-P section of File S1 for the proofs. We actually show these results hold with extremely high probability, which is a stronger notion that implies probability tending to 1.

1. Let v be a node of MGEO-P(n, m, α, β, p) with rank R that arrived at step t. Then

deg(v) = i -1 n -1 p 1 -α n 1-α-β + (n -i)pR -α n -β •   1 + O   log 2 (n) n 1-α-β     .
This result implies that the degree distribution follows a powerlaw with exponent η = 1 + 1 α . 2. The average degree of node of MGEO-P(n, m, α, β, p) is

ρ = p 1 -α n 1-α-β   1 + O   log 2 (n) n 1-α-β     . 3. The diameter of MGEO-P(n, m, α, β, p) is n Θ( 1 m ) .
This last property suggests that, ignoring constants, for a network with n nodes and diameter D, the expected dimension based on the MGEO-P model is

m ≈ log n log D .
Thus, like some network models that incorporate geometry [START_REF] Kim | Multiplicative attribute graph model of real-world networks[END_REF][START_REF] Zhang | High-dimensional Apollonian networks[END_REF], in the MGEO-P model, the dimension m must scale logarithmically in order for the diameter to remain constant as n increases.

Experimental Design and Graph Summaries

Both graph motifs and spectral densities are numeric summaries of a graph that abstract the details of a network into a small set of values that are independent of the particular nodes of a network. These summaries have the property that isomorphic graphs have the same values, and we will use these summaries to determine the dimension of the metric space that best matches Facebook and LinkedIn networks as illustrated in Figure 2. Graph motifs, graphlets, or graph moments are the frequency or abundance of specific small subgraphs in a large network. We study undirected, connected subgraphs up to four nodes as our graph motifs (with the exception of the number of edges, or two-node motifs, as the networks are created to preserve this count). This is a set of 8 graphs shown in at the bottom of Figure 2 along with the single two node graph of an edge. The spectral density of a graph is the statistical distribution of eigenvalues of the normalized Laplacian matrix as indicated in the upper right of that figure.

These eigenvalues indicate and summarize many network properties including the behavior of a uniform random walk, the number of connected components, an approximate connectivity measure, and many other features [START_REF] Chung | Spectral Graph Theory[END_REF][START_REF] Banerjee | Graph spectra as a systematic tool in computational biology[END_REF]. Thus, the spectral density of the normalized Laplacian is a particularly helpful characterization that captures many such separate network properties.

We study dimensional scaling in social networks by comparing samples of the MGEO-P networks of varying dimensions with samples of social network data from Facebook and LinkedIn. We pay particular attention to the relationship between the number of nodes n of the network and the dimension m of the best fit MGEO-P network. In order to determine what underlying dimension for MGEO-P best fits a given graph, we employ two distinct methods. For one experiment, we use features known as graph motifs, graphlets, or graph moments in concert with a support vector machine (SVM) classifier. This approach has been used successfully to determine the best generative mechanism of a network [START_REF] Memišević | An integrative approach to modeling biological networks[END_REF] and to select parameters of a complicated network models to fit real-world data [START_REF] Gleich | Moment based estimation of stochastic Kronecker graph parameters[END_REF][START_REF] Moreno | Learning mixed Kronecker product graph models with simulated method of moments[END_REF]. In a second experiment, we use spectral densities of the normalized Laplacian matrix of a graph and a Kullback-Leibler divergence (KL divergence) similarity measurement, which has been used to match protein networks between species [START_REF] Patro | Global network alignment using multiscale spectral signatures[END_REF][START_REF] Banerjee | Structural distance and evolutionary relationship of networks[END_REF]. We find evidence of the logarithmic dimension hypothesis in both cases.

The data

Facebook distributed 100 samples of social networks from universities within the United States measured as of September 2005 [START_REF] Traud | Social structure of Facebook networks[END_REF], which range in size from 700 nodes to 42,000 nodes. We call these networks the Facebook samples. The LinkedIn samples were created from the LinkedIn connection network together with the creation time of each connection from May 2003 to October 2006. To perform our experiments on networks of different size, we build 71 snapshots of the LinkedIn network at various timestamps. We then extracted a dense subset of their graph at various time points that is representative of active users; we used the 5-core of the network for this purpose [START_REF] Seidman | Network structure and minimum degree[END_REF]. The k-core of a network is a maximum size subset of vertices such that all vertices have degree k. See Figure 3 and the full statistics tables of File S1 for additional properties of these networks. In both networks, the number of edges per node grows at essentially the same rate.

Results

The results of our dimensional fitting for graphlets are shown in Figure 4 and the results of the fitting using spectral densities are in Figure 5. For both datasets and both types of statistics, the best-fit dimension scales logarithmically with the number of nodes and closely tracks a simple model prediction based on the diameter D of the network (the model curve plots m = log(n)/ log(D)). These experiments corroborate the logarithmic dimension hypothesis; although the precise fits differ as shown in Table 2.

The most important feature of these results is that both methodologies show similar scaling in how the dimensionality scales with network size. There are minor differences between the precise predicted dimensions-for instance, the spectral density approach predicts slightly higher dimensions for Facebook than does the graphlet approach-but the results agree to a reasonable degree with the dimension predicted by the model: log(n)/ log(D). Also, the confidence bounds are small around the chosen dimension.

Sensitivity and robustness

We investigate the sensitivity of the graphlet results in two settings. If we reduce the training set size of the SVM classifier by using a random subset of 20% of the input training data and then rerun the training and classification procedure 50 times, then we find a distribution over dimensions that we report as a box-plot, shown in Figure 4. In File S1, we further study perturbation results that argue against these results occurring due to chance. In particular, we find that these dimensions are robust to moderate changes to the network structure (File S1, Figure S2) and we find that our methodology does not predict useful dimensions of Erdős-Rényi random graphs or random graphs with the same degree distribution (File S1, Figure S1). We do not report a precise p-value as there are no widely accepted null-models for network data. We study the sensitivity of the spectral densities that look for matches that are within 105% of the true minimum divergence. This defines a dimension interval around each match that is small for all of our examples.

Discussion

There is a growing body of evidence that argues for some type of geometric structure in social and information networks. An important study in this direction views networks as samples of geometric graphs within a hyperbolic space [START_REF] Krioukov | Hyperbolic geometry of complex networks[END_REF][START_REF] Krioukov | Network cosmology[END_REF][START_REF] Krioukov | Duality between equilibrium and growing networks[END_REF]. Recent work has further shown that hyperbolic embeddings reproduce shortest path metrics in real-world networks [START_REF] Zhao | Efficient shortest paths on massive social graphs[END_REF]. In both MGEO-P and hyperbolic random geometric networks, highly skewed or power-law degree distributions are imposed-either directly as in MGEO-P, or implicitly as in the hyperbolic space scaling. These results further support hidden metric structures in networks by empirically confirming a prediction about the dimension of the metric space made by one particular model. The importance of this finding is that it provides new insight into how the metric space must behave as the network grows. Previous studies assume a fixed dimension metric structure and our results indicate that a variable dimension may be more appropriate. In practice, estimating the dimensions of these networks could be useful for anomaly detection in the network and characterizing new types of network data.

Note that these results do not conclusively argue that MGEO-P is a perfectly accurate model for social networks; there are meaningful differences between the spectral histograms from MGEO-P and real social networks, see Figure 6. There are also similar differences in the graphlet counts. Our results support a different hypothesis. The closest MGEO-P network to a given social network has a metric space whose dimension scales logarithmically with the number of nodes. In File S1 (Sensitivity studies section), we have determined that this property is not due to either the edge density or the degree distribution; thus, our findings appears to reflect a new intrinsic property of social networks.

This finding suggests a number of opportunities for designing social network models with metric spaces that evolve in time. We believe that such models offer the opportunity to identify new properties of social network based on emergent properties of the models. One question to address is how the metric space and connection radius change, if at all, as the network grows. Answering this question would provide insight into the value of additional users of a network. Additionally, our results suggest that many network models that assume a fixed dimension should be reevaluated.

Materials and Methods

Powerlaw fitting

To determine the powerlaw exponent η, we use the Clauset-Shalizi-Newman power-law exponent estimator [START_REF] Clauset | Power-law distributions in empirical data[END_REF] as implemented by Tamás Nepusz [START_REF] Nepusz | plfit software[END_REF].

Diameters

The MGEO-P model of a network predicts that the dimension m should approximate log(n)/ log(D), where D is the diamater. However, as D is sensitive to outliers we use the 99% effective diameter computed via an asymptotically accurate approximation scheme [START_REF] Palmer | ANF: A fast and scalable tool for data mining in massive graphs[END_REF] as implemented in the SNAP library on 2011-12-31. The effective diameter of all Facebook networks ranges between 3.5 and 4.6, with a mean of 4.1. For the LinkedIn data, the effective diameter ranges between 4.3 and 5.9, with a mean of 5.4. In both networks, larger graphs have bigger effective diameters, although the differences are slight and the full data is available in the File S1, Full statistics tables.

Graphlets

To compute graphlets, we employ the rand-esu sampling algorithm [START_REF] Wernicke | Efficient detection of network motifs[END_REF] as implemented in the igraph library [START_REF] Csardi | The igraph software package for complex network research[END_REF]. This algorithm approximates the count of each subgraph via a stochastic search, which then depends on the probability of continuing to search. Thus, if the probability is near 1 then the scores are nearly exact, but very expensive to compute, and small probabilities truncate the search early to produces fast estimates. The value we use is 10/n. We use log-transformed output from this procedure in order to capture the dynamic range of the resulting values.

Spectral densities

We approximate the spectral density via a 201-bin histogram of the eigenvalues of the normalized Laplacian, which all fall between 0 and 2. (The choice of 201 was based on prior experiences with the spectral histograms of networks.) To compute eigenvalues of a network, we employ the recently developed ScaLA-PACK routine using the MRRR algorithm [START_REF] Dhillon | A new O(n 2 ) algorithm for the symmetric tridiagonal eigenvalue/eigenvector problem[END_REF][START_REF] Dhillon | The design and implementation of the MRRR algorithm[END_REF][START_REF] Vömel | ScaLAPACK's MRRR algorithm[END_REF].

SVM

We used a multi-class support-vector machine (SVM) based classification tool from Weka [START_REF] Witten | Data Mining: Practical machine learning tools and techniques[END_REF] to predict the relationship between the graphlets and the dimension. We considered alternatives, such as alternating decision trees and logistic regression; however, we settled on the SVM approach as it has the most flexible classification boundary to fit the highly nonlinear relationships between graphlet counts and dimensions.

Setting MGEO-P Parameters

Consider a graph G = (V, E) that we wish to compare to an MGEO-P sample. The MGEO-P model depends on four parameters: n, m, α, and β. The choice of n is straightforward as we use the number of nodes of the original graph. Both α and β can be chosen independently of the dimension m. Specifically, both α and β determine the average degree of the network and the exponent of the power law in the degree distribution, up to lower-order terms, as shown by property 1 and property 2. By computing just these two simple statistics of a network-the exponent of the power law and the average degree-we can invert these relationships and choose these parameters. Let η be the power-law exponent and ρ be the average degree. Then:

α + β = 1 -log(ρ)/ log(n) and α = 1 η-1 . In order to derive this simple expression, we make the simplifying assumption that p does not go to zero too quickly, for example p = n -o (1) , in which case: log(rho) = (1 -α -β) log(n) + o(1) follows from the expression for the average degree of a MGEO-P network. We use the following treatment of the probability p in order to maximize the clustering coefficient of the network. We first generate an MGEO-P network with p = 1. Then suppose that the original network had E = nρ/2 edges, we continue by randomly deleting edges until the output has exactly the same number of edges as the input network. This step can be interpreted as using the value of p necessary to get the same edge count as the original graph. In the case where there are insufficient edges, we leave the output from the MGEO-P generator untouched. This process effectively chooses p as large as possible, which gives us the largest local clustering. An example describing the MGEO-P process on a graph with 250 nodes in the unit square with torus metric, where α = 0.9 and β = 0.04 and p = 1. Each figure shows the graph "replicated" in grey on all sides in order to illustrate the torus metric. Links are drawn to the closest replicated neighbor. The blue square indicates the region [0, 1] 2 . Top row (left to right) The MGEO-P process begins with relatively few nodes, and thus, nodes must have large influence radii (red squares) to link anywhere. As more nodes arrive, large radii result in many connections, modeling influential users, and small radii result in a few connections, modeling standard users. Bottom row Illustrates the final constructed graph. We see similar scaling for both types of networks, but with slightly different offsets. For Facebook, log 10 (edges) = 1.06 log 10 (nodes) + 1.35 with R 2 = 0.945; for LinkedIn log 10 (edges) = 1.07 log 10 (nodes) + 0.56 with R 2 > 0.999. The regularity in the LinkedIn sizes is due to our construction of those networks. Figure 4. Facebook dimension at left, LinkedIn dimension at right. Each red dot (SVM) is the predicted dimension computed via graphlet features and a support vector machine classifier. For the Facebook data, we find that m = 2.06 log(n)/ log(10) -3.00. For the LinkedIn data, we find that m = 0.7333 log(n)/ log(10) + 1. And these are plotted as the red linear fit line. Our theoretical model predicts a dimension of log(n)/log(D) and we plot this as the dashed line. In each figure, we show the variance in the fitted dimension as a box-plot. We estimate the variance by using only 20% of the original training data and repeating over 50 trials. There are only a few outliers for small dimensions. We also show any other other dimensions within 5% of this divergence value. The dimensions shift modestly higher for Facebook and remain almost unchanged for LinkedIn. Both still are closely correlated with the theoretical prediction based on the model based on log(n)/ log(D) (dashed line). The linear fits to the predicted dimensions is plotted as the red linear fit line. Figure 6. For three of the Facebook networks, we show the eigenvalue histogram in red, the eigenvalue histogram from the best fit MGEO-P network in blue, and the eigenvalue histograms for samples from the other dimensions in grey. The MGEO-P model correctly captures the peak of the distribution around 1, but fails to completely capture the tail between 1 and 2. Thus, we see meaningful difference between these profiles and hence, do not suggest that MGEO-P captures all of the properties of real-world social networks.

Tables Table 1: The parameters of the MGEO-P model n the total number of nodes m the dimension of the metric space 0 < α < 1 the attachment strength parameter 0 < β < 1 -α the density parameter 0 < p ≤ 1 the connection probability The specific dimensional scaling lines fit to the data in Figures 4 and5 illustrate the growth of the network is logarithmic in the number of nodes.

Figure 1 .

 1 Figure 1.An example describing the MGEO-P process on a graph with 250 nodes in the unit square with torus metric, where α = 0.9 and β = 0.04 and p = 1. Each figure shows the graph "replicated" in grey on all sides in order to illustrate the torus metric. Links are drawn to the closest replicated neighbor. The blue square indicates the region [0, 1] 2 . Top row (left to right) The MGEO-P process begins with relatively few nodes, and thus, nodes must have large influence radii (red squares) to link anywhere. As more nodes arrive, large radii result in many connections, modeling influential users, and small radii result in a few connections, modeling standard users. Bottom row Illustrates the final constructed graph.

Figure 2 .

 2 Figure2. At left and center, we have the steps involved in fitting via graphlets; at right and center, we have the steps involved in fitting via spectral histogram. Throughout, red lines denote the flow of features for the MGEO-P networks whereas blue lines denote flow of features for the original networks. At the bottom, we show an enlarged representation of the 8 graphlets we use.

Figure 3 .

 3 Figure 3. The scale of the network data involved in our study varies over three orders of magnitude. We see similar scaling for both types of networks, but with slightly different offsets. For Facebook, log 10 (edges) = 1.06 log 10 (nodes) + 1.35 with R 2 = 0.945; for LinkedIn log 10 (edges) = 1.07 log 10 (nodes) + 0.56 with R 2 > 0.999. The regularity in the LinkedIn sizes is due to our construction of those networks.

Figure 5 .

 5 Figure 5. Facebook data at left, LinkedIn data at right. Each blue point (Eigen) is the dimension of the MGEO-P sample with the minimum KL-divergence between the graph and the MGEO-P sample. We also show any other other dimensions within 5% of this divergence value. The dimensions shift modestly higher for Facebook and remain almost unchanged for LinkedIn. Both still are closely correlated with the theoretical prediction based on the model based on log(n)/ log(D) (dashed line). The linear fits to the predicted dimensions is plotted as the red linear fit line.

Table 2 :

 2 Dimension scaling for Facebook and LinkedIn

	Data	Dimension fit	Coefficients		95% Confidence
		a log(n)/ log(10) + b a	b	a	b
	Facebook Graphlet	2.06	-3.00 (1.851, 2.264)	(-3.821, -2.182)
		Spectral density	1.21	1.65	(0.9782, 1.446)	(0.7272, 2.578)
	LinkedIn Graphlet	0.98	1.01	(0.786, 1.178)	(0.1591, 1.87)
		Spectral density	0.77	1.1	(0.56, 0.99)	(0.23, 1.95)

Acknowledgments

We acknowledge Jure Leskovec for allowing us to access the LinkedIn dataset.

Supporting Information Legends

File S1: This supporting document contains the following components of our analysis. (i) Formal proofs of the MGEO-P properties. (ii) Full statistical information about each of the Facebook and LinkedIn networks including the graphlet counts. (iii) Figure S1: Predicted dimensions of random graphs with the same degree distribution. (iv) Figure S2: The change in predicted dimension found by perturbing the graph structure. (v) A discussion of the sensitivity results about the predicted dimension.