Xavier Descombes

Grégoire Malandain

Caroline Fonta

Lazlo Negyessy

Rajmund Mosko

Lázló Négyessy

Automatic dendrite spines detection from X-ray tomography volumes

Keywords: Dendritic spine, X-Ray micro-tomography, 3D Hough transform, Point Process

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

INTRODUCTION

Dendritic spines are small bulbous structures connected to the stem dendrite through a tiny neck. They have a huge biological importance as they represent the majority of the synaptic connections of the cerebral cortex providing the backbone of the cortical connectome [START_REF] Yuste | Dendritic spines and distributed circuits[END_REF]. Studying their density and repartition along dendrites provide essential information to explain neural network functioning, and pathology of diverse disorders [START_REF] Fiala | Dendritic spine pathology: cause or consequence of neurological disorders?[END_REF][START_REF] Irwin | Dendritic spine structural anomalies in fragile-x mental retardation syndrome[END_REF][START_REF] Tackenberg | Thin, stubby or mushroom: spine pathology in alzheimer's disease[END_REF][START_REF] Lee | Cocaine-induced dendritic spine formation in d1 and d2 dopamine receptor-containing medium spiny neurons in nucleus accumbens[END_REF]. In addition, their repartition according to a morphologic type (mushroom, stubby or thin) provides more detailed quantative biomarkers.

The first step to study spine repartition consists in detecting and segmenting the dendritic spines in neurons images. In the literature, imaging protocols generally rely on fluorescence microscopy techniques, either confocal scanning microscopy [START_REF] Rodriguez | Automated three-dimensional detection and shape classification of dendritic spines from fluorescence microscopy images[END_REF], 2-photons scanning microscopy [START_REF] Cheng | A novel computational approach for automatic dendrite spines detection in twophoton laser scan microscopy[END_REF][START_REF] Zhang | Dendritic spine detection using curvilinear structure detector and lda classifier[END_REF], or electron microscopic reconstruction. Such protocols provide contrasted images at high resolution that facilitate the detection process. However, the obtained images are usually anisotropic, with a slice thickness that may be much larger than the pixel size. As a consequence, some spines may be missed. In addition, the imaged field of view of the scanning microscopy technique is quite small (typically one hundred of 512x512 slices with a pixel size of about 0.05 µm): a larger number of samples is required for statistical analysis.

In this paper, we consider an alternative imaging protocol, namely the X-rays micro-tomography (XRMT). It overcomes the two main drawbacks of the scanning microscopy protocols, since the obtained 3D images are isotropic (the voxel size is identical along the 3 dimensions), and it allows to image larger samples (typically the obtained images are of size 1000 3 with a voxel size of 0.308µm 3 ). Moreover, in contrast to scanning microscopy, no occlusion or masking by overlapping structures occur in the deeper slices (structures closer to the objective lens occlude deeper ones beneath).

One advantage of fluorescence microscopy is due to the potential use of dedicated labels, which are cellular specific. However, they require either intracellular staining or use of transgenic animals. Moreover, photo-bleaching is clearly a limitation to such staining. Our staining technique (Golgi [START_REF] Fairen | A new procedure for examining golgi impregnated neurons by light and electron microscopy[END_REF]) is more versatile, and is clearly sufficient to describe the fundamental brain circuits with a high specificity (Ramon y Cajal describes the different morphologic spine types with such staining at the end of the 19th century).

XMRT also exhibits a higher level of noise than laser scanning microscopy. This implies that popular existing approaches (e.g. SpineLab [START_REF] Jungblut | Spinelab: tool for three-dimensional reconstruction of neuronal cell morphology[END_REF] or Neuronstudio [START_REF] Rodriguez | Automated three-dimensional detection and shape classification of dendritic spines from fluorescence microscopy images[END_REF]) are likely to fail. Hence, robust algorithms are mandatory. Moreover, some of the existing methods requires user interaction, which is not compatible with the large volumes produced by XMRT. Therefore, a fully automated approach has to be designed.

Before defining a strategy to detect dendritic spines, let us notice some difficulties linked with the problem: (i) spine size is close to the data resolution, (ii) spine may appear not connected with the dendrite. Therefore we have to avoid (i) denoising (as proposed in [START_REF] Cheng | A novel computational approach for automatic dendrite spines detection in twophoton laser scan microscopy[END_REF]) that may remove some spine, (ii) basing the detection on the dendrite tree extraction (as proposed in [START_REF] Rodriguez | Automated three-dimensional detection and shape classification of dendritic spines from fluorescence microscopy images[END_REF]) that may omit some spines and (iii) detecting spines on the whole dataset that may introduce many false positives (as proposed in [START_REF] Cheng | A novel computational approach for automatic dendrite spines detection in twophoton laser scan microscopy[END_REF]). Besides, we expect to make the detection in 3D but not on the Maximum Intensity Projection as in [START_REF] Choy | Multi scale and slice-based approach for automatic spine detection[END_REF]. To circumvent these difficulties we propose to use shape information and localization with respect to dendrites to detect spines. We first extract dendrites as tubular structures. These tubular features are used to define a spatial mask where spine may appear. A spine detector embedding shape information is finally applied on this mask.

We detail the method in section 2 and show some result in section 3 before concluding in section 4.

METHODS

Extracting dendrite medial axis

First, the medial axis of the dendrites are detected thanks to a multi-scale Hessian-based method [START_REF] Krissian | Model-based detection of tubular structures in 3d images[END_REF]. Basically, at each scale σ the Hessian matrix (i.e. the matrix of the image second derivatives) is calculated by a convolution of the image I with the second derivatives of a 3D Gaussian kernel of standard deviation σ. At the center of a tubular structure, the Hessian matrix exhibits two high eigenvalues, whose associated eigenvectors v 1 and v 2 are orthogonal to the direction of the tubular structure, and one low eigenvalue, whose associated eigenvector obviously gives the direction of the tubular structure. Once the likely direction of the dendrite is identified at point M , a function response is calculated that integrates the radial derivative along a circle of radius σ:

R σ (M ) = 1/2π 2π 0 v(θ).∇ σ I(M + σv(θ))dθ (1) with v(θ) = cos θv 1 + sin θv 2 .
The responses at each scale are then combined into a multiscale response, thanks to normalized γ-derivatives [START_REF] Lindeberg | Scale-Space Theory in Computer Vision[END_REF], and last the local extrema of this response, with respect to the directions v 1 and v 2 are extracted.

Dendrite detection

To denoise the obtained dendrite medial axis image, we consider a segment approximation by computing the 3D Hough transform. To parametrize lines in the 3D space, we consider a sphere including the data volume. Any line crossing the sphere is then parametrized by the four angles corresponding to the two intersection points between the line and the sphere. Any voxel P above a threshold T , in the medial axis image will vote for the following points in the 4D Hough space (see figure 1):

{(θ 1 , φ 1 , θ 2 , φ 2 ), θ 1 , θ 2 ∈ [-π, π], φ 1 , φ 2 ∈ [0, π 2 ], P S 1 S 1 S 2 , } , (2)
where S 1 (resp, S 2 ) is the point on the sphere with angles

(θ 1 , φ 1 ) (resp.(θ 2 , φ 2
)) The candidate lines are obtained by extracting the local maxima in the Hough space. Finally, the dendrites are approximated from the candidate lines by considering the segments constituted of voxels whose gray value Fig. 1. The voting mechanism for the 3D Hough transform is higher than the threshold T and having a minimal length l. Once the segments describing the dendrite axis have been obtained, we reconstruct the dendrites with a conditional dilation M D = (S D K 1 ) ∩ Y T , where is the Minkowsky addition and Y T is the data thresholded at T . As dendrites are mainly vertically oriented, we consider a disk in the XY plane for the structuring element K 1 , the size of which is defined by the data resolution and the typical dendrites diameter (three voxels in our case as the typical dendrite diameter is 1µm [14]).

Spine detection

Once we have detected dendrites we define a mask where spines may lie. This mask is obtained using a dilation of the dendrite mask as follows:

M S = M D K 2 -M D , (3)
where M D is the dendrite mask, the Minkowsky addition. The size of the structuring element K 2 is defined by the resolution and the maximum distance between a spine and its dendrite (again, a XY disk of size three voxels is used as a structuring element corresponding to the typical size of spine stalk [START_REF] Jacobson | Dimensions of the dendritic spine in the sensorimotor cortex of the rat, cat, squirrel monkey and man[END_REF]). The spine mask M S allows a reduction of the computation time in the detection process but also avoids false positives due to irregularities on dendrites.

To detect spines we consider a point process in M S . A configuration of points is denoted x = {x 1 , ..., x n }, x i ∈ M S . This probabilistic framework embeds geometric constraints through the definition of a data term and structural constraints through the definition of a prior. It has been successfully applied for multiple objects detection problems such as trees, flamingos or cells detection. We define the point process by an unormalized density h with respect to the Poisson measure [START_REF] Descombes | Stochastic geometry for image analysis[END_REF] as follows:

h(x) = f (x)g(x|y), (4)
where f (x) is the prior and g(x|y) the data term depending on data y. The prior is a hard core process:

f (x) = (i,j),i =j δ(d(x i , x j ) > D), (5)
where δ(a) = 1 if a is true and 0 otherwise. Therefore the prior states that the minimal distance between two spines must be greater than D (D = 3 voxels in our case). To compute the data term we consider that a spine is included in a 3×3×3 cube (the typical spine bulb size being 1µm in mammals [START_REF] Jacobson | Dimensions of the dendritic spine in the sensorimotor cortex of the rat, cat, squirrel monkey and man[END_REF]) and consider the contrast between this cube and its neighborhood by computing the Bhattacharyya distance:

b(x i ) = 1 4 (µ in -µ out ) 2 σ 2 in + σ 2 out + 1 2 log σ 2 in + σ 2 out 2 σ 2 in σ 2 out δ( µ in µ out > 1), (6)
where µ in and σ 2 in (resp. µ out and σ 2 out ) are the mean and variance of the set {y s : d(i, s) < 2} (resp. {x s : 2 ≤ d(i, s) < 3}). Finally, the data is written as follows:

g(x|y) = exp - i u(x i ) , (7)
where:

u(x i ) = 1 -b(xi) b0 if b(x i ) < b o exp -(b(xi)-bo) 3b0 otherwise ( 8
)
b 0 is a contrast parameter that defines spine candidates when the corresponding Bhattacharyya distance is greater than b 0 .

To obtain the optimal solution with respect to h(.), we consider a fast approximation of the Multiple Birth and Death algorithm [START_REF] Descombes | Object extraction using a stochastic birth-and-death dynamics in continuum[END_REF]:

Birth: Consider as candidates the set C of all points in M S Sort: Sort the candidates x i ∈ C with respect to u(x i ) in the decreasing order (the worst first), to obtained x = {x 1 , ..., x n } Death: For i = 1 to n, if h(x) < h(x/{x i }) remove x i from x

RESULTS

To exemplify the proposed approach, a subvolume (220 × 180 × 100) has been extracted from a XRMT volume that is given on figure 2. As expected, the spines appear as small objects, whose size is close to the image resolution, along the tubular structures representing dendrites (see figure 3). Using the localization information to detect spine is essential to prevent false positives due to noise or to the deviation of dendrites from a cylinder model (see figure 4 for a potential false positive avoided thanks to its distance with respect to dendrites). The parameters linked with the geometry (size of the mask and structuring elements) are automatically fixed considering the image resolution and the dendrite spine size referred in the litterature [START_REF] Jacobson | Dimensions of the dendritic spine in the sensorimotor cortex of the rat, cat, squirrel monkey and man[END_REF]. The remaining parameters (threshold in the 3D Hough transform and in the detection step) are data dependant as they are linked with the image contrast. Next step will consists in estimating these parameters from the image histogram. Figure 5 shows the detected dendrite medial axis and the spine mask derived from the dendrite detection. The obtained spine detection is shown on figure 6. The obtained results are promising and correspond to a visual inspection of the data. We are currently performing a manual detection by several experts to provide a quantitative evaluation of these results.

CONCLUSION

We have proposed an automated algorithm for detecting dendritic spines from XRMT data. XRMT data allows imaging a large volume of tissue, and therefore a higher number of spines than laser scanning microscopy. We have shown that despite the lower image quality compared to microscopic data, we were able to extract dendritic spines. The main idea of the proposed approach is to define a mask for performing the spine detection without facing the false positive problem. We therefore first extract the dendrites themselves and then compute the spine mask based on prior knowledge on their localization with respect to dendrites. Results are visually promising, and a forthcoming validation study will allow to better assess the quality of the detection by providing a quantitative evaluation. In a next step, we will consider the problem of spine classication in the three classical classes: stubby, thin, mushroom.

Fig. 2 .

 2 Fig. 2. Original XMRT data. Left: axial slice (69µm×57µm) ; right: sagittal slice (69µm × 31µm) .

Fig. 3 .

 3 Fig. 3. Examples of spines. Left: constrasted spine; right: poorly contrasted spine.

Fig. 4 .

 4 Fig. 4. Noise similar to spine but too far from dendrite.

Fig. 5 .

 5 Fig. 5. Definition of a spine mask. Left: medial axis (perspective view from front side); right: dendrite in white and spine mask in gray (perspective view from top).

Fig. 6 .

 6 Fig. 6. Dendrites (in green) and spines (in red)


Luca Castelli Aleardi

Olivier Devillers
email: olivier.devillers@inria.fr

Eric Fusy
email: fusy@lix.polytechnique.fr

Luca Castelli

Éric Fusy

Canonical Ordering for Triangulations on the Cylinder, with Applications to Periodic Straight-line Drawings

come

Introduction

The problem of efficiently computing straight-line drawings of planar graphs has attracted a lot of attention over the last two decades. Two combinatorial concepts for planar triangulations turn out to be the basis of many classical straight-line drawing algorithms: the canonical ordering (a special ordering of the vertices obtained by a shelling procedure) and the closely related Schnyder wood (a partition of the inner edges of a triangulation into 3 spanning trees with specific incidence conditions). Algorithms based on canonical ordering [START_REF] De Fraysseix | How to draw a planar graph on a grid[END_REF][START_REF] Kant | Drawing planar graphs using the canonical ordering[END_REF] are typically incremental, adding vertices one by one while keeping the drawing planar. Algorithms based on Schnyder woods [START_REF] Schnyder | Embedding planar graphs on the grid[END_REF] are more global, the (barycentric) coordinates of each vertex have a clear combinatorial meaning (typically the number of faces in certain regions associated to the vertex). Algorithms of both types make it possible to draw in linear time a planar triangulation with n vertices on a grid of size O(n × n). They can also both be extended to obtain (weakly) convex drawings of 3-connected maps on a grid of size O(n × n). The problem of obtaining planar drawings of higher genus graphs has been addressed less frequently [START_REF] Kocay | Drawing graphs on the torus[END_REF][START_REF] Gortler | Thurston Discrete one-forms on meshes and applications to 3D mesh parameterization[END_REF][START_REF] Mohar | Straight-line representations of maps on the torus and other flat surfaces[END_REF][START_REF] Mohar | Tessellation and visibility representations of maps on the torus[END_REF][START_REF] Chambers | Drawing Graphs in the Plane with a Prescribed Outer Face and Polynomial Area[END_REF][START_REF] Duncan | Planar drawings of higher-genus graphs[END_REF][START_REF] Zitnik | Drawing graphs on surfaces[END_REF], from both the theoretical and algorithmic point of view. Recently some methods for the straight-line planar drawing of genus g graphs with polynomial grid area (of size O(n 3 ), in the worst case) have been described in [START_REF] Chambers | Drawing Graphs in the Plane with a Prescribed Outer Face and Polynomial Area[END_REF][START_REF] Duncan | Planar drawings of higher-genus graphs[END_REF] (to apply these methods the graph needs to be unfolded planarly along a cut-graph). However, it does not yield (at least easily) periodic representations: for example, in the case of a torus, the boundary vertices (on the boundary of the rectangular polygon) might not be aligned, so that the drawing does not give rise to a periodic tiling.

Our main contribution is to generalize the notion of canonical ordering and the incremental straight-line drawing algorithm of de Fraysseix, Pach and Pollack [START_REF] De Fraysseix | How to draw a planar graph on a grid[END_REF] (shortly called FPP algorithm thereafter) to triangulations on the cylinder. For any triangulation G of the cylinder, our algorithm yields in linear time a crossing-free straight-line drawing of G on a regular grid (on the flat cylinder) of the form Z/wZ × [0..h], with w ≤ 2n and h ≤ n(2d + 1), where n is the number of vertices of G and d is the (graph-) distance between the two boundaries of G. By a reduction to the cylindric case (the reduction is done with the help of a so-called tambourine [START_REF] Bonichon | Edge partition of toroidal graphs into forests in linear time[END_REF]), we also get a drawing algorithm on the torus. Precisely, for any toroidal triangulation G, we can obtain in linear time a crossing-free straight-line drawing of G on a regular grid (on the flat torus) Z/wZ × Z/hZ, with w ≤ 2n and h ≤ 1 + n(2c -1), where n is the number of vertices of G and c is the length of a shortest non-contractible cycle. Since c ≤ (2n) 1/2 as shown in [START_REF] Albertson | On the independence ratio of a graph[END_REF], we have h ≤ (2n) 3/2 , so that the grid area is O(n 5/2 ).

For the toroidal case we mention that a notion of canonical ordering has been introduced in [START_REF] Castelli-Aleardi | Schnyder woods for higher genus triangulated surfaces, with applications to encoding[END_REF] (this actually works in any genus and yields an efficient encoding procedure) but we do not use it here. We also mention that, independently, an elegant periodic straight-line drawing algorithm for toroidal triangulations has been very recently described in [START_REF] Gonçalves | Toroidal maps : Schnyder woods, orthogonal surfaces and straight-line representation[END_REF], based on so-called toroidal Schnyder woods and face-counting operations; in their case the area of the periodic grid is O(n 4 ).

Preliminaries

Graphs embedded on surfaces. A map of genus g is a connected graph G embedded on the compact orientable surface S of genus g, such that all components of S\G are topological disks, which are called the faces of the map. The map is called planar for g = 0 (embedding on the sphere) and toroidal for g = 1 (embedding on the torus). The dual of a map G is the map G * representing the adjacencies of the faces of G, i.e., there is a vertex v f of G * in each face f of G, and each edge e of G gives rise to an edge e * = {v f , v f } in G * , where f and f are the faces on each side of e. A cylindric map is a planar map with two marked faces B 1 and B 2 whose boundaries C(B 1 ) and C(B 2 ) are simple cycles (possibly C(B 1 ) and C(B 2 ) share vertices and edges). The faces B 1 and B 2 are called the boundary-faces. Boundary vertices and edges are those belonging to C(B 1 ) (black circles in Fig. 1) or C(B 2 ) (gray circles in Fig. 1); the other ones are called inner vertices (white circles in Fig. 1) and edges.

Periodic drawings of cylindric triangulations

A cylindric triangulation is a cylindric map with no loops nor multiple edges and such that all non-boundary faces are triangles (in Section 5 we will allow such triangulations to have 1-cycles or 2-cycles separating the two boundary-faces). We introduce at first a notion of canonical ordering for cylindric triangulations:

Definition 1.
Let G be a cylindric triangulation with boundary-faces B 1 and B 2 , and such that the cycle C(B 1 ) has no chords (i.e., there is no edge that is not on C(B 1 ) and has both ends on C(B 1 )). An ordering π = {v 1 , v 2 , . . . , v n } of the vertices of G\C(B 1 ) is called a (cylindric) canonical ordering if it satisfies:

• For each k ≥ 0 the map G k induced by C(B 1 ) and by the vertices {v 1 , . . . , v k } is a cylindric triangulation.

• The other boundary-face B of G k (the one different from B 1 ), whose contour is denoted C k , contains B 2 (to clarify the meaning of "B contains B 2 ", remember that a cylindric map is a special kind of map on the sphere; since G k is a submap of G, any face of G k is made of a set of faces of G).

• The vertex v k is on C k , and its neighbours in G k-1 are consecutive on C k-1 . The notion of canonical ordering makes it possible to construct a cylindric triangulation G incrementally, starting from G 0 = C(B 1 ) and adding one vertex (and incident edges) at each step. This is similar to canonical orderings for planar triangulations, as introduced by de Fraysseix, Pach and Pollack [START_REF] De Fraysseix | How to draw a planar graph on a grid[END_REF] (the main difference is that for a planar triangulation one starts with G 0 being an edge, called the base-edge, whereas here one starts with G 0 being a cycle, seen as a cylindric map without non-boundary faces). Shelling procedure. We now describe a shelling procedure to compute a canonical ordering of a cylindric triangulation G with boundary-faces B 1 , B 2 . At each step the graph formed by the remaining vertices is a cylindric triangulation, one boundary face remains B 1 all the way, while the other boundary-face (initially B 2 ) has its contour, denoted by

B 1 B 2 v C 6 C 7 v e v w 2 w 1 w 1 w 1 w 1 F F * edge of F
C k , getting closer to C(B 1 ). A vertex v ∈ C k is free if v is incident to no chord of C k and if v / ∈ C(B 1 ) (see Fig. 2 top left).
The shelling procedure goes as follows (n is the number of vertices in G\C(B 1 )): for k from n to 1, choose a free vertex v on C k , assign v k ← v, and then delete v together with all its incident edges. The existence of a free vertex at each step follows from the same argument as in the planar case [START_REF] Brehm | 3-orientations and Schnyder 3-Trees decompositions[END_REF]. Indeed, if there is no chord incident to C k , then any vertex v ∈ C k \C(B 1 ) is free, while if there is a chord e for C k , then the set of chords incident to C k forms a system of archs (relative to C k ). If we look at a chord e = {u, v} that is "innermost" (i.e., seen as an arch, no other arch is nested inside), then the path between u and v on C k contains at least one vertex, which has to be free (see Fig. (each one along a path in the dual forest) are inserted in order to make the slopes of e and e r smaller than 1 in absolute value. Then the new vertex and its edges connected to the upper boundary can be drawn in a planar way.

Underlying forest and dual forest. Given a cylindric triangulation G (with boundary faces B 1 and B 2 ) endowed with a canonical ordering π, define the underlying forest F for π as the oriented subgraph of G where each vertex v ∈ C(B 2 ) has outdegree 0, and where each v / ∈ C(B 2 ) has exactly one outgoing edge, which is connected to the adjacent vertex of v of largest label in π. The forest F can be computed on the fly during the shelling procedure: when treating a free vertex v k , for each neighbour v of v k such that v / ∈ C k , add the edge {v, v k } to F , and orient it from v to v k . Since the edges are oriented in increasing labels, F is an oriented forest; it spans all vertices of G\C(B 2 ) and has its roots on C(B 2 ). The augmented map G is obtained from G by adding a vertex w 1 inside B 1 , a vertex w 2 inside B 2 , and connecting all vertices around B 1 to w 1 and all vertices around B 2 to w 2 (thus triangulating the interiors of B 1 and B 2 , see Fig. 2 bottom middle). Define F as F plus all edges incident to w 1 and all edges incident to w 2 . Define the dual forest F * for π as the graph formed by the vertices of G * (the dual of G) and by the edges of G * that are dual to edges not in F . Since F is a spanning connected subgraph of G, F * is a spanning forest of G * . Precisely, each of the trees (connected components) of F * is rooted at a vertex "in front of" each edge of B 1 , and the edges of the tree can be oriented toward this root-vertex (see Fig. 2

bottom right).

Drawing algorithm. Given a cylindric triangulation G such that C(B 1 ) has no chord, we first compute a canonical ordering of G, and then draw G in an incremental way. We start with a cylinder of width 2|C(B 1 )| and height 0 (i.e., a circle of length 2|C(B 1 )|) and draw the vertices of C(B 1 ) equally spaced on the circle (space 2 between two consecutive vertices). Then the strategy for each k ≥ 1 is to compute the drawing of G k out of the drawing of G k-1 by first stretching the cylinder (increasing the width by 2) and then placing the vertex v k and its incident edges (in G k ) in a planar way. Define the x-span of an edge e in the cylindric drawing as the number of columns [i, i + 1] × [0, +∞] that meet the interior of e (we have no need for a more complicated definition since, in our drawings, a column will never meet an edge more than once).

Consider the dual forest F * for the canonical ordering restricted to G k-1 . Note that the set of vertices of C k-1 that are neighbours of v k form a path on C k-1 . Traversing this path ). The vertices are treated in increasing label (the canonical ordering is the one computed in Fig. 2).

γ with the boundary-face of G k-1 (the one different from B 1 ) to the left, let e be the first edge of γ and e r be the last edge of γ. Let also a k be the starting vertex and let b k be the ending vertex of γ. Let P (resp. P r ) be the path in F * from e * (resp. e * r ) to the root in its connected component (which is a vertex "in front of" an edge of B 1 )). We stretch the cylinder by inserting a vertical strip of length 1 along P and another along P r (see Fig. 3). This comes down to increasing by 1 the x-span of each edge of G k dual to an edge in P , and then increasing by 1 the x-span of each edge dual to an edge in P r (note that P and P r are not necessarily disjoint, in which case the x-span of an edge dual to an edge in P ∩ P r is increased by 2). After these stretching operations, whose effect is to make the slopes of e and e r strictly smaller than 1 in absolute value, we insert (as in the planar case) the vertex v k at the intersection of the ray of slope 1 starting from a k and the ray of slope -1 starting from b k , and we connect v k to all vertices of γ by segments 1 . These two rays actually intersects at a grid point since the Manhattan distance between any two vertices on C k-1 is even. Fig. 4 shows the execution of the algorithm on the example of Fig. 1.

Proposition 1. For each cylindric triangulation G with no chord incident to C(B 1 ), one can compute in linear time a crossing-free straight-line drawing of G on an x-periodic regular grid Z/wZ × [0..h] where -with n the number of vertices of G and d the (graph-)distance between the two boundaries-w = 2n and h ≤ n(2d + 1), such that: every v ∈ C(B 1 ) has 1 In the FPP algorithm for planar triangulations, the step to make the (absolute value of) slopes of e and e r smaller than 1 is formulated as a shift of certain subgraphs described in terms of the underlying forest F . The extension of this formulation to the cylinder would be quite cumbersome. We find the alternative formulation with strip insertions more convenient for the cylinder (it also gives rise to a very easy linear-time implementation). y(v) = 0 (so every edge in C(B 1 ) has slope 0), and every edge belonging to C(B 2 )\C(B 1 ) has slope ±1.

Proof. The fact that the drawing remains crossing-free relies on the slope-property for the upper boundary and on the following inductive property, which is easily shown to be maintained at each step k from 1 to n:

Pl: for each edge e on C k (the upper boundary of G k ), let P e be the path in F * from e * to the root, let E e be the set of edges dual to edges in P e , and let δ e be any nonnegative integer. Then the drawing remains planar when successively increasing by δ e the x-span of all edges of E e , for all e ∈ C k .

We now prove the bounds on the grid-size. If |C(B 1 )| = t then the initial cylinder is 2t × 0; and at each vertex insertion, the grid-width grows by 2. Hence w = 2n. In addition, due to the slope conditions (slopes of boundary-edges are at most 1 in absolute value), the y-span (vertical span) of every edge e is not larger than the current width at the time when e is inserted in the drawing. Hence, if we denote by v the vertex of C(B 2 ) that is closest (at distance d) to C(B 1 ), then the ordinate of v is at most d • (2n). And due to the slope conditions, the vertical span of C(B 2 ) is at most w/2 ≤ n. Hence the grid-height is at most n(2d + 1). The linear-time complexity is shown next.

Linear-time implementation. An important remark is that, instead of computing the xcoordinates and y-coordinates of vertices in the drawing, one can compute the y-coordinates of vertices and the x-span of edges (as well as the knowledge of which extremity is the leftend vertex and which extremity is the right-end vertex). In a first pass (for k from 1 to n) one computes the y-coordinate of vertices and the x-span r e of each edge e ∈ G at the time t = k when it appears on G k (as well one gets to know which extremity of e is the left-end vertex). Afterwards if e / ∈ F , the x-span of e might further increase due to insertion of new vertices. Precisely, let v j be a vertex inserted afterwards (i.e., j > t), and (with the notations a j , b j of the drawing algorithm description) let = {v j , a j } and r = {v j , b j }. Note that e is stretched due to the insertion of the strip along P iff * is in the subtree T e of F * formed by the edges descending from e * . Similarly e is stretched due to the insertion of the strip along P r iff * r is in T e . To state it more clearly, each edge in T e is responsible for an increase (by 1) of the x-span of e. Hence the total x-span of each edge e ∈ G is given by r e + s e , where s e = 0 if e ∈ F , and, if e / ∈ F , s e is the number of edges in T e . Since all quantities s e can easily be computed in linear time, this gives a linear-time implementation.

Allowing for chordal edges at C(B 1 ). We finally explain how to draw a cylindric triangulation when allowing for chords incident to C 0 = C(B 1 ); it is good to view B 2 as the top boundary-face and B 1 as the bottom-boundary face (and imagine a standing cylinder). For each chordal edge e for the cycle C 0 , the component under e is the face-connected part of G that lies below e; such a component is a quasi-triangulation (polygonal outer face, triangular inner faces) rooted at the edge e. A chordal edge e of C 0 is maximal if the component Q e under e is not strictly included in the component under another chordal edge for C 0 . The size of such an edge e is defined as |e| = 2|V (Q e )| -4. (the size |e| is actually the width of the drawing of Q e using the FPP algorithm). If we delete the component under each maximal chordal edge (i.e., delete everything from the component except for the chordal edge itself) we get a new bottom cycle C 0 that is chordless, so we can draw the reduced cylindric triangulation G using the algorithm of Proposition 1. As we have seen in Section 3 (implementation paragraph), for each edge e of C 0 , the initial x-stretch is r e = 2 and then the further increase s e of the x-stretch equals the number of edges descending from e * in the dual forest F * . Note that we have actually some freedom to choose the initial x-stretch r e of each edge e ∈ C 0 (just it has to be a positive even number since at each step of the incremental algorithm the vertices of the current upper boundary have to be at even Manhattan distance). If e ∈ C 0 is on C 0 we take r e = 2. If e ∈ C 0 is not on C 0 (i.e., e was a maximal chordal edge for C 0 ), we take for r e the minimal even positive number such that r e + s e ≥ |e|, i.e., r e = 2 • max(1, (|e| -s e )/2 ). Hence, at the end of the execution of the drawing of G , the length e = r e + s e of each maximal chord e satisfies e ≥ |e|. Then for each maximal chord e of C 0 , we draw the component Q e under e using the FPP algorithm. This drawing has width |e|, with e as horizontal bottom edge of length |e| and with the other outer edges of slopes ±1. We shift the left-extremity of e so that the drawing of Q e gets width (e), then we rotate the drawing of Q e by 180 degrees and plug it into the drawing of G (see Fig. 5). The overall drawing of G is clearly planar. We obtain: Theorem 1. For each cylindric triangulation G, one can compute in linear time a crossingfree straight-line drawing of G on an x-periodic regular grid Z/wZ × [0, h], where -with n the number of vertices and d the (graph-) distance between the two boundaries-w ≤ 2n and h ≤ n(2d + 1). The drawing is x-monotone (the intersection with any vertical line is an interval) and the slopes of boundary-edges are at most 1 in absolute value.

Periodic drawings of toroidal triangulations

Let G be a toroidal triangulation. A non-contractible cycle of G is a cycle of edges of G that does not delimit a topological disk. Note that a pair Γ 1 , Γ 2 of oriented noncontractible cycles that are parallel (i.e., are homotopic and are oriented in a parallel way) delimits a cylindric triangulation T , which is formed by the faces to the right of Γ 1 and to the left of Γ 2 . The pair Γ 1 , Γ 2 is called a tambourine if the edges dual to those of T \{Γ 1 , Γ 2 } form a non-contractible cycle Γ that is homotopic to Γ 1 and Γ 2 (see dark-blue faces and dashed edges in Fig. 6 top-middle). The edges of T \{Γ 1 , Γ 2 } are said to be inside the tambourine. It can be shown (see [START_REF] Bonichon | Edge partition of toroidal graphs into forests in linear time[END_REF], the master's thesis of the third author Arnaud Labourel, and the next paragraph) that for each non-contractible cycle Γ of G, there exists a tambourine whose two cycles are parallel to Γ. Deleting the edges that are strictly inside the tambourine, one obtains a cylindric triangulation G with Γ 1 and Γ 2 as the contours of the boundary-faces. Note also that the distance d between Γ 1 and Γ 2 is smaller than the length of a shortest non-contractible cycle not parallel to Γ. We now apply the algorithm of Theorem 1 to G . If we augment the height h of the drawing to h = h + w + 1, and then wrap the x-periodic grid Z/wZ × [0..h] into a periodic grid Z/wZ × Z/h Z, and finally insert the edges inside the tambourine as segments2 , then the slope properties (edges on Γ 1 and Γ 2 have slope at most 1 in absolute value while edges inside the tambourine have slopes greater than 1 in absolute value) ensure that the resulting drawing is crossing-free (see Fig. 6). Observe that we can choose Γ so that the graph-distance between the two boundaries Γ 1 and Γ 2 (in G ) is smaller than the length γ of a shortest non-contractible cycle in G; and this choice for Γ can be done without computing a shortest non-contractible cycle. Indeed, let {Γ a , Γ b } be a basis of non-contractible cycles (computable in linear time, using for instance a cut-graph).

Denoting by Γ min a shortest non-contractible cycle of G, for sure at least one of Γ a or Γ b is not parallel to Γ min . Hence, for Γ a or for Γ b , the distance between the boundary-cycles (after deleting edges of the parallel tambourine) is smaller than |Γ min |. In other words if we choose the one cycle (among {Γ a , Γ b }) that yields the smaller distance between the two boundaries of G , then this distance will be smaller than γ. We obtain:

Theorem 2. For each toroidal triangulation G, one can compute in linear time a crossingfree straight-line drawing of G on a periodic regular grid Z/wZ × Z/hZ, where (with n the number of vertices and γ the length of a shortest non-contractible cycle) w ≤ 2n and h ≤ 1 + n(2γ + 1). Since γ ≤ √ 2n (as shown in [START_REF] Albertson | On the independence ratio of a graph[END_REF]), the grid area is O(n 5/2 ).

Existence of a tambourine. For the sake of completeness we include a proof of existence of a tambourine, which slightly extends the proof given in the master's thesis of Arnaud Labourel. A toroidal map is called weakly 3-connected if its periodic representation in the plane is 3-connected. Let G be such a map and let Γ be a non-contractible cycle of G.

We are going to show that G has a tambourine parallel to Γ. Let G be the cylindric map obtained after cutting G along Γ; we take the annular representation of G , calling Γ 1 (resp. Γ 2 ) the copy of Γ that is the outer (resp. inner) boundary. Let Γ be the smallest (in terms of the enclosed area) cycle that strictly encloses Γ 2 (i.e., encloses Γ 2 and is vertex-disjoint from Γ 2 ). Let Γ be the largest (in terms of the enclosed area) cycle that is strictly enclosed in Γ (i.e., is enclosed by Γ and is vertex-disjoint from Γ ). Note that by minimality Γ has no chord inside, and by maximality Γ has no chord outside. Hence, if we can show that there is no vertex in the area A (strictly) between Γ and Γ , then we can conclude that, in G, Γ and Γ form a tambourine parallel to Γ. Assume there is a vertex v in A. Call vertex of attachment for Γ a vertex w ∈ Γ such that there is a path from v to w visiting only vertices of A before reaching w. Again by minimality of Γ it is easy to see that there is a unique vertex of attachment v for Γ . Similarly (by maximality of Γ ) there is a unique vertex of attachment v for Γ . Then (again by minimality of Γ and maximality of Γ ) there is a closed curve γ that meets G only at the vertices v , v , and such that the interior of γ contains v but does not contain any of the two boundary-faces. Such a curve γ yields a 2-separator in the periodic representation of G, a contradiction.

Allowing for non-contractible 1-and 2-cycles

For a cylindric map, a non-contractible cycle is a cycle that separates the two boundary-faces (i.e., there is a boundary-face on one side and a boundary-face on the other side of the cycle).

Other cycles are called contractible. A weakly simple cylindric triangulation is a cylindric map whose non-boundary faces are triangles, whose contractible cycles have length at least 3, and such that each vertex has at most one incident non-contractible loop. A weakly-simple toroidal triangulation is a toroidal map with triangular faces, with all contractible cycles of length at least 3, and such that any two (non-contractible) loops incident to the same vertex can not be homotopic. These are the necessary and sufficient conditions for the triangulation (whether on the cylinder or on the torus) to have no loops nor multiple edges in the periodic representation; hence these are the conditions under which one can aim at a periodic crossing-free straight-line drawing. For weakly simple cylindric triangulations without loops (nor chords incident to C(B 1 )), exactly the same shelling procedure and iterative drawing algorithm can be taken as for simple cylindric triangulations. In case there are loops we have to explain how to deal with them (see Fig. 7 for an example). For the shelling procedure, if the current upper boundary C k is a loop -call v the incident vertex-then one deletes the loop and immediately takes v as the next free vertex (the fact that v is free is due to the fact that there is no other loop at v). In the drawing procedure (how to insert v and its incident loop into the drawing), one first adds v without its loop (by a classical one-step iteration of the drawing algorithm, involving two strip insertions), and then one draws the loop at v as an horizontal segment spreading over the whole width of the current periodic drawing (the loop is added without inserting any vertical strip). Finally one can deal with chords incident to C(B 1 ) in the same way as for simple cylindric triangulations. About weakly simple toroidal triangulated maps, the procedure is also the same as for simple toroidal triangulations since the above proof of existence of a tambourine holds in that case. And the grid bounds (for the cylindre and for the torus) are the same as for simple triangulations.

Figure 1 :

 1 Figure 1: A cylindric triangulation with boundary faces B 1 = {a, b, c, d} and B 2 = {h, i, j, k}. Left: annular representation. Right: x-periodic representation.

Figure 2 :

 2 Figure 2: Shelling procedure to compute a canonical ordering of a given cylindric triangulation. The underlying forest is computed on the fly; the last drawing shows the underlying forest superimposed with the dual forest. The graph is the one of Fig. 1.

  2 top right).

Figure

  Figure One step of the incremental drawing algorithm. Two vertical strips of width 1 (each one along a path in the dual forest) are inserted in order to make the slopes of e and e r smaller than 1 in absolute value. Then the new vertex and its edges connected to the upper boundary can be drawn in a planar way.

Figure 4 :

 4 Figure4: Algorithm (Prop. 1) to compute an x-periodic drawing of a cylindric triangulation (no chordal edges incident to B 1 ). The vertices are treated in increasing label (the canonical ordering is the one computed in Fig.2).

e 2 e 1 Figure

 1 Figure Drawing a cylindric triangulation with chords at B 1 . To make enough space to place the component under e 2 , one takes 8 (instead of 2) as the initial x-span of e 2 .

1 Γ 2 ΓFigure 6 :

 126 Figure 6: The main steps for drawing of a toroidal triangulation: 1) remove the edges inside a tambourine, (2) draw the obtained cylindric triangulation, 3) insert the edges of the tambourine back into the drawing.

2 2 1 Figure 7 :

 217 Figure 7: Top line: canonical ordering of a weakly simple cylindric triangulation (with some non-contractible 1-and 2-cycles). Middle line: iterative drawing algorithm. Bottom line: dealing with chordal edges.

  the flat cylinder of width w and height h is the rectangle [0, w] × [0, h] where the vertical sides are identified. A point on this cylinder is located by two coordinates x ∈ R/wZ and y ∈ [0, h]. The flat torus of width w and height h is the rectangle [0, w] × [0, h] where both pairs of opposite sides are identified. A point on this torus is located by two coordinates x ∈ R/wZ and y ∈ R/hZ. Assume from now on that w and h are positive integers. For a cylindric triangulation G, a periodic straight-line drawing of G of width w and height h is a crossing-free straight-line drawing (edges are drawn as segments, two edges can meet only at common end-points) of G on the flat cylinder of width w and height h, such that the vertexcoordinates are in Z/wZ × [0..h] (i.e., are integers). Similarly for a toroidal triangulation G, a periodic straight-line drawing of G of width w and height h is a crossing-free straight-line drawing (edges are drawn as segments, two edges can meet only at common end-points) of G on the flat torus of width w and height h, such that the vertex-coordinates are in Z/wZ × Z/hZ (i.e., are integers).

We insert the edges in the tambourine T in the unique way such that, looking from bottom to top, at least one edge in T goes strictly to the right, and all edges going strictly to the right have x-span at most w; in this way it is easy to check that the x-span of all edges in T is at most w.

Acknowledgments. The authors thank D. Gonçalves and B. Lévêque, and (independently) B. Mohar for interesting discussions about toroidal Schnyder woods, and N. Bonichon for explanations on the computation of a tambourine in a toroidal triangulation.

* This work has been supported by European project ExploreMaps -ERC StG 208471.


Ghislain Landry

Tsafack Chetsa

Laurent Lefevre
email: laurent.lefevre@ens-lyon.fr

Jean-Marc Pierson
email: pierson@irit.fr

Patricia Stolf
email: stolf@irit.fr

Georges Da Costa
email: dacosta@irit.fr

A Runtime Framework for Energy Efficient HPC Systems Without a Priori Knowledge of Applications

Keywords: energy efficiency, green leverage, phase recognition, execution vector

The rising computing demands of scientific endeavours often require the creation and management of High Performance Computing (HPC) systems for running experiments and processing vast amounts of data. These HPC systems generally operate at peak performance, consuming a large quantity of electricity, even though their workload varies over time. Understanding the behavioural patterns (i.e., phases) of HPC systems during their use is key to adjust performance to resource demand and hence improve the energy efficiency.

In this paper, we describe (i) a method to detect phases of an HPC system based on its workload, and (ii) a partial phase recognition technique that works cooperatively with on-the-fly dynamic management. We implement a prototype that guides the use of energy saving capabilities to demonstrate the benefits of our approach. Experimental results reveal the effectiveness of the phase detection method under real-life workload and benchmarks. A comparison with baseline unmanaged execution shows that the partial phase recognition technique saves up to 15% of energy with less than 1% performance degradation.

I. INTRODUCTION

The increasing reliance of scientific endeavours on computing has made High Performance Computing (HPC) mainstream in many scientific areas, including climate research, disease control, and drug discovery. Although HPC systems deliver tremendous peak performance for solving challenges faced by the humanity, their electrical consumption and high operational cost have become major concerns. As a consequence, different from a decade ago when supercomputers were ranked only by their peak performance 1 , nowadays they are also assessed based on their energy efficiency 2 .

Improving the energy efficiency of HPC is not trivial as it cannot use classical energy saving capabilities employed in other types of systems because they often reduce the peak performance; which is usually unacceptable in HPC. For example, reducing the CPU frequency or switching off nodes may reduce energy consumption, but depending on the workload this reduction comes at the expense of decreasing performance. It is indispensable to choose the right energy saving capability to the workload being executed.

During its life-cycle an HPC application commonly exhibits several behaviours, hereafter referred to as phases or regions, during which techniques can be applied to reduce the electrical consumption of the computing infrastructure. For example, by knowing that the system is running a compute intensive task, a resource manager can spin disks down or put the network interfaces in low power consumption mode. Discovering the task dependency graph can enable optimisations such as reducing processor frequency when running a task that is not in the critical path. This a priori knowledge about applications is often regarded important to improve the energy efficiency of HPC systems. In reality, however, HPC systems are generally shared by multiple applications with heterogeneous computational behaviours, where an optimisation for saving energy considering a single application is likely to impact the performance of others. In addition, the increasing complexity of HPC applications makes it challenging to optimise the energy efficiency per application (a detailed discussion is provided in Section II).

In this work, we overcome this complexity by improving the energy efficiency considering the HPC infrastructure itself (which can be shared by multiple applications) rather than focus on optimising individual applications. We propose and implement an online methodology for phase detection and identification in HPC systems. We then introduce a partial phase recognition technique that guides the usage of green capabilities/leverages. We define a green capability as any action that can save energy in an HPC system, such as: CPU frequency scaling, spinning down disks, scaling the speed of network interconnections, and task migration. An important feature of our methodology is that it does not require prior knowledge of the applications running on the system.

The paper is organised as follows: Section II presents related work. Section III describes our online framework, whereas Section IV presents and analyzes experimental results, and validates the system prototype. Finally, Section V concludes the paper and discusses future work.

II. RELATED WORK

Over the past few years, several energy saving capabilities have been proposed; they can work either at the hardware level or at the software level. Considering hardware, most vendors provide equipments with techniques to modify the performance of processors, network interfaces, memory, and I/O during their operation, thus allowing for reducing the energy consumed by HPC subsystems. For instance, modern processors are provided with Dynamic Resource Sleeping (DRS), which makes components hibernate to save energy and wakes them up on demand. Although progress has been made, advances in hardware solutions to reduce energy consumption have been slow due to the high costs of equipment design and the increasing demand for performance. Our work combines a set of hardware technologies to reduce the energy consumption of HPC systems.

Unlike hardware approaches, software solutions for reducing the energy usage of HPC systems have gained momentum with the research community. Rountree et al. [START_REF] Rountree | Adagio: making DVS practical for complex HPC applications[END_REF] use node imbalance to reduce the overall energy consumption of a parallel application. They track successive MPI communication calls to divide the application into tasks composed of a communication portion and a computation portion. A slack occurs when a processor waits for data to arrive during a task execution. The speed of the processor can be slowed down with almost no impact on the overall execution time of the application. Rountree et al. developed Adagio to track task execution slacks and compute the appropriate frequency at which a processor should run. Although the first instance of a task is always run at the highest frequency, further instances of the task are executed at the frequency computed after the first execution. In addition, a tool called Jitter [START_REF] Kappiah | Just in time dynamic voltage scaling: exploiting inter-node slack to save energy in mpi programs[END_REF] was developed to detect slack moments in performance to perform inter-node imbalance and use DVFS to adjust the CPU frequency. Different from Adagio, our fine-grained data collection can differentiate not only computation-intensive and communication-intensive execution portions (which we call phases/regions) but also memory-intensive portions. Memory-intensive phases can be run on a slower core without paying a significant performance penalty [START_REF] Kotla | Characterizing the Impact of Different Memory-Intensity Levels[END_REF].

Online techniques have been used to detect and characterise application execution phases, and set the appropriate CPU frequency [START_REF] Isci | Live, runtime phase monitoring and prediction on real systems with application to dynamic power management[END_REF], [START_REF] Choi | Fine-grained dynamic voltage and frequency scaling for precise energy and performance tradeoff based on the ratio of off-chip access to onchip computation times[END_REF]. The techniques rely on hardware monitoring counters to compute runtime statistics such as cache hit/miss ratio, memory access counts, and retired instructions counts; which are then used for phase detection and characterisation. The developed policies tend to be designed for single task environment [START_REF] Isci | Live, runtime phase monitoring and prediction on real systems with application to dynamic power management[END_REF], [START_REF] Choi | Fine-grained dynamic voltage and frequency scaling for precise energy and performance tradeoff based on the ratio of off-chip access to onchip computation times[END_REF]. We overcome this limitation by treating each node of a cluster as a black box, which means that we do not focus on any application, but instead on the platform. The flexibility provided by this assumption enables us to track not the applications/workloads execution phases, but nodes' execution phases. Our work also differs from previous work by using partial phase recognition instead of phase prediction. Online recognition of communication phases in an MPI application was investigated by Lim et al. [START_REF] Lim | Adaptive, transparent frequency and voltage scaling of communication phases in mpi programs[END_REF]. Once a communication phase is recognised, the authors apply CPU DVFS to save energy. They intercept and record the sequence of MPI calls during program execution and consider a segment of program code to be reducible if there are highly concentrated MPI calls or if an MPI call is long. The CPU is then set to run at the appropriate frequency when the reducible region is recognised again.

To the best of our knowledge, our work differs from those above in two major ways. First, our phase detection approach does not rely on a specific HPC subsystem or MPI communication calls. Second, our model does not focus on saving only the processor energy; it also takes advantage of other power saving capabilities available on HPC subsystems.

The closest research to this work is probably our previous paper [START_REF] Tsafack | DNA-inspired Scheme for Building the Energy Profile of HPC Systems[END_REF] in which we developed a cross platform methodology for detecting and characterising phases in HPC systems. Although both efforts use performance counters, they differ significantly; in our previous work we modelled the entire system's runtime as a state graph associated with a transition matrix. The transition matrix helps determine the next state in the graph. The problem with our previous approach is that it takes time to get a useful transition matrix since the system may enter different configurations (which we call letters) from one execution to another, which may prevent it from working well under a system with varying workloads. Moreover, the entire model is built offline.

III. METHODOLOGY A. Global Approach

The rationale behind our methodology is that it is possible to save energy while maintaining performance by selecting the most suitable green capabilities for the system at a phase. Approaches described in Section II set the CPU frequency according to the estimated use of the processor over a time interval. In addition to scaling the CPU frequency, an HPC system can potentially benefit from other energy reduction schemes conceived to adjust the system to the actual demand,. These schemes include adjusting the speed of the network interconnect; switching off memory banks; spinning down disks; migrating tasks across nodes.

A typical HPC system exhibits many different phases or behaviours during runtime, complicating the choice of green capabilities. Characterising phases at runtime, so that similar patterns can be easily identified in the future, is hence worthwhile. With this characterisation, a set of green capabilities that maintain the performance at a reasonable level (in terms of execution time) and improve the energy efficiency at given phase can be employed later during similar phases. Our approach associates each phase with a characterisation and a set of green capabilities. A detailed analysis of phase characterisation is provided in Section III-B.

As mentioned earlier, phase characterisation and identification enable reuse of green capabilities suitable to one phase by another phase. However, a phase can be considered similar to a previous phase once it completes, at which time applying green capabilities related to the former might not produce the expected results. The literature uses prediction to identify an upcoming phase, which requires a priori information about the workload being executed; not appropriate for systems shared by multiple applications. Therefore, as we do not have any a priori knowledge of the workload, we introduce a heuristic referred to as partial phase recognition or simply partial recognition.

Instead of trying to recognise a complete phase prior to adjusting the system (which would then be useless since the phase would have finished), we decide to adjust the system when a fraction of a phase has been recognised. This technique clearly gives false positives (i.e. sections incorrectly recognised as part of a phase), but we argue that adjusting the system is beneficial at least for a certain time. When the current phase diverges too much from the recognised phase, then another phase can be identified or a new phase characterised.

Further details on phase identification and partial recognition are given in Section III-C, whereas Section III-D details how green capabilities are selected, explains how they are linked to phases, and discusses how the proposed approach can be enforced in a distributed system with independent nodes.

B. Defining and Characterizing Phases

Our approach is based on partial recognition of phases during the system activity. A system phase also known as phase, represents the runtime behaviour of the system over a given interval during which specific metrics are stable; meaning that their corresponding values remain below a given threshold within the phase. These metrics can be application or system specific and are computed for each node using data about the computational state of the node.

To illustrate the idea, let us assume we only have one metric, the load of the system. Consider a system whose load, measured every second, is displayed in Figure 1; where the x-axis shows the time in seconds and the y-axis the load. With a threshold set to 5%, we can easily recognise 4 phases (following the vertical lines in the graphic, from time 0 to 2, from time 3 to 6, from time 7 to 9, and from time 10 to 13). Each of these phases is characterised by its duration and the different loads the system went through during that phase. The number of phases obviously varies with the threshold, for with a threshold of let's say 50% we would have had 2 phases (from time 0 to 9, and from time 10 to 13).

As a general rule, a single metric (the load of the system in our example) is not adequate for characterising useful phases. The set of data chosen will be small enough to save space, but adequately large so as to distinguish between different phases.

As stated beforehand, each phase is associated with its duration and data collected during the phase. The size of the data representing a phase is proportional to its duration, which implies that large amounts of data might be collected for very long phases. To prevent this, a digest is used to summarise the data. We define a signature as a representative for a phase, where the signature is computed from multiple values. In the following, the signature will be computed as the vector of averaged monitored values. In the previous example, the signature of a phase can be the average load during the phase.

C. Phase Identification and Partial Recognition

1) Execution Vectors: In order to apply any of the green capabilities, it is worth detecting and identifying appropriate phases. An Execution Vector (EV) is a vector of values that represent the behaviour of an application during a time interval. The length of these intervals is defined depending on the goal reactivity. It can be of the order of the mean time to perform an action (switching on or off a node is an action that takes minutes, whereas changing a processor's frequency is at the second time-scale). EVs are considered as points in the positive quadrant of the execution space spanned by the N vector dimensions; and their values must give a maximum insight on the current execution. These values can be taken from the system (load, network-or I/Orelated) or directly from the hardware using performance counters. Hardware performance counters are special hardware registers originally designed to performance tuning. However, they have recently been used for investigating the power usage of applications [START_REF] Contreras | Power prediction for intel xscale processors using performance moni-toring unit events[END_REF], [START_REF] Singh | Real time power estimation and thread scheduling via performance counters[END_REF]. According to the literature, the most important counters are bus cycles, number of instructions, cache misses, cache hits, cache references, and stall cycles (they may have different names depending on the architecture).

2) Phase Detection and Identification: Following the definition of a phase in Section III-B, we can assume without loss of generality that during a phase execution vectors are stable.

Figure 2 displays the similarity matrix using execution vectors collected when a cluster was running a synthetic benchmark which successively runs Lower-Upper symmetric Gauss-Seidel (LU), Embarrassingly Parallel (EP), Conjugate Gradient (CG) and Integer Sort (IS) from NPB-3.3 [START_REF] Bailey | The nas parallel benchmarks[END_REF] benchmark. A detailed explanation of similarity matrices using power vectors can be found in the literature [START_REF] Isci | Identifying program power phase behavior using power vectors[END_REF]. As for power vectors, they represent for the whole execution the distance between all pairs of EVs. At coordinate i, j in the upper matrix of Figure 2 the colour represents on a grey scale the distance between the EV at time i and j respectively.

Consequently, the colour at the point i, j tends to black when the distance between EV at i and EV at j tends to zero. Similarly, the colour at the position i, j tends to white as the distance between EV at i and EV at j increases. Above the diagonal line (the diagonal represents the execution timeline from the upper left corner to the lower right corner), we can report 7 triangular blocs. These blocs represent phases the system went through when executing the workload (our synthetic benchmark). From Figure 2, we can also observe that distances between execution vectors within the same phase approach zero, whereas in between phases the distances approach one.

The above observations led us to the premise that a phase change occurs when the distance between two consecutive EVs goes beyond a given threshold, which we refer to as the detection threshold.

While phase detection can serve any kind of system optimisation, keeping all execution vectors is not realistic. Hence, once a phase completes, it is associated a signature vector. The signature vector of a phase is a vector obtained by computing the arithmetic average of each metric (e.g., load, network, or a particular performance counter) of the vectors computed during the phase. It is used for phase identification, where a phase is considered similar to another if the distance between their signature vectors is below the detection threshold. If the distance between the signature vectors of the new phase and an already identified phase is below the detection threshold, then we do nothing; otherwise a new phase is added via its signature vector to the signature vectors list. Put simply, only the first occurrence of a phase is stored.

3) Partial Recognition: Partial recognition technique is used to identify phases before their completion. Once a phase is completed, in addition to its signature vector and its characteristics, we store its reference vector and its length. The reference vector of a phase is a randomly chosen execution vector among those pertaining to the set of vectors sampled during the phase. Thus, an ongoing phase (a phase started but not completed yet) P ongoing is recognised as an existing phase P exist if the manhattan distance between each execution vector of the already executed part of P ongoing and the reference vector of P exist is below a given threshold defining the percentage of dissimilarity between them. We state that the recognition threshold is X% if the already executed part of P ongoing equals X% of the length of P exist .

D. Using Phases to Save Energy

Having identified a phase based on information of a previously characterised phase, data on the latter is used for choosing the most relevant green capabilities. The following explains how this approach can be enforced in a distributed system with multiple nodes.

1) Link Between Phases and State: The stability of metrics within a phase implies that the behaviour of the system in a phase is more or less constant; meaning that there is a predominant behaviour within the phase. That predominant behaviour is considered as the state of the phase. Referring to the literature, a phase can be cpu-intensive, memoryintensive, and IO-intensive (including network accesses).

To find out the predominant behaviour of a phase or its state, we apply Principal Component Analysis (PCA) to the dataset made up by execution vectors belonging to the corresponding phase. PCA is typically issued when several underlying factors shape the data. Therefore, we use variables (each variable represents a metric) contributing less to the first principal axis of PCA to determine the predominant behaviour over the phase (components of the execution vector are access rate of our sensor set).

2) Enforcing Green Leverage: The goal is to optimise the energy consumption of HPC systems under the following constraints: (a) energy optimisation must be done without any specific information about the characteristics of the application prior to its execution; (b) management policies for reducing the overall system energy consumption must not have a significant impact on performance.

Given these constraints, the most suitable capabilities have to be selected for each phase. Knowing the state of a phase and its remaining time before completion, system adaptation 3) Global Architecture: The diagram of Figure 3 shows the steps of our detection making process. Partial recognition starts as soon as a new phase has begun. When the partial recognition threshold is reached (the ongoing phase is identified with an existing phase), the characteristics (including its duration) of the existing phase are used to adapt the system accordingly. In a distributed system with independent nodes we propose a centralised architecture with a global coordinator. The motivation for the centralised management is that task migration, which is one of our green capabilities, requires a global view of the system. Figure 4 depicts the global architecture. The decision making node implements the decision making process; it sends commands to appropriate nodes after decision making. On each node of the system a process samples one execution vector per second and sends it to the global coordinator.

IV. EXPERIMENTATION AND VALIDATION

A. Phase Detection and Partial Recognition Results 1) Phase Detection Results: Phase detection is critical to any management policy that intends to take advantage of workload variability. Analysis of phase detection using power vectors [START_REF] Isci | Identifying program power phase behavior using power vectors[END_REF] and basic bloc vectors [START_REF] Sherwood | Automatically Characterizing Large Scale Program Behavior[END_REF] have been provided in the literature. However, these techniques use offline phase analysis. We rely on EVs which are comparable to power vectors [START_REF] Isci | Identifying program power phase behavior using power vectors[END_REF] used in power phase behaviour identification.

The execution vector represents the access rate of a set of performance counters along with the disk access rate (read/write) and network access rate (send/receive) per second. Each new EV is normalised with respect to a maximum vector. The maximum vector is special in the sense that it comprises the maximum available value of the elements of all execution vectors.

For phase detection, we use the manhattan distance between two consecutive execution points; we then state that two executive EVs belong to the same phase or have the same behaviour if the manhattan distance between them is below a given threshold. The threshold is actually a fixed percentage of the maximum known distance between all consecutive EVs. This maximum is updated any time a greater value is found. Our experiments show that a threshold of 5% is effective in many cases. Figure 5 where doted vertical lines represent the identified phases shows the results of our online phase detection with a 5% threshold on one node of our cluster. During the experiment, our cluster was running the Advance Weather Research and Forecasting (WRF-ARW) [START_REF] Skamarock | A description of the Advanced Research WRF Version 2[END_REF] model. WRF-AWR is globally one phase broken up with very short transition phases. However, as our phase detection process looks at the runtime behaviour (via execution vectors) of each node of the platform and considers as a phase any behaviour between two consecutive behavioural changes, transition phases will cause WRF-AWR to break into many micro-phases.

To demonstrate the effectiveness of our online detection approach, we implemented using execution vectors, the offline phase detection approach presented by Isci et Online phase detection (with a 5% threshold) and partial recognition (with a 10% threshold) on one of the node running WRF-ARW.

al. [START_REF] Isci | Identifying program power phase behavior using power vectors[END_REF] on the dataset used for online phase detection. Figure 7 displays in a grey scale the upper diagonal matrix of distances between all pairs of execution vectors. Execution vectors are grouped considering an empirical threshold of 5%, which means that the colour at the position i,j in the graphic tends to white when the manhattan distance between execution vectors at i and j is greater than 5% of the maximum existing distance. Following the phase detection methodology described in Section III-C2; it can be seen from Figure 7 that our online phase detection is as effective as the offline approach.

2) Partial Phase Recognition: We assign to each newly identified phase a reference vector that is further used for partial phase identification employing the algorithm described in Section III-C. The reference vector is a random vector chosen from the set of vectors belonging to the phase. In Figure 6, similar phases are represented by rectangles filled with the same pattern. For clarity, phases that were partially identified (or identified using partial recognition) start with a small rectangle filled in in black. The length of the rectangle filled in black represents the part of the phase that was used for partial recognition. For this experiment, the recognition threshold was set to an empirical value of 10%.

B. Real-life System Implementation: Measurement and Evaluation Platform

We implemented our online phase detection and recognition and tested them on a cluster system. The implemented prototype monitors the overall cluster's behaviour via performance monitoring counters along with network bytes sent/received counts and disk read/write counts. For convenience, we use the term "sensors" to designate performance monitoring counters along with network bytes sent/received counts and disk read/write counts. Each node of the cluster sends a sample of data to the decision making node on a per-second basis. Samples collected corresponding to execution vectors are further used by the decision maker node to perform phase detection and partial recognition for any individual node. Partial recognition is more suitable in our case because phases are not detected at fixed intervals; besides, the system may be shared by many applications.

The use of green capabilities available on the managed cluster is guided by our partial phase recognition. When partial recognition is successful, the characteristics of the phase the newly started phase is identified with are translated into a predefined set of "green settings". The system adaptation implemented so far uses CPU frequency scaling as a proof of concept, but we plan to implement many other green leverages/capabilities including task migration, disk sleep state, and network interconnect speed scaling.

1) Dynamic Power Management via Green Leverage/Capabilities: Our phase characterisation relies on the assumption that sensors contributing less to the first principal axis of PCA give more information about what did not happen during the phase. Thus we arbitrarily select 5 sensors among those contributing less to the first principal axis of PCA for phase characterisation. In addition, we compute the average instruction rate of each phase, which is then used for idle phase characterisation. We need such a metric for characterising idle phases as over the course of idle phases almost all sensors are likely to be constant. Table I shows an example of phases translation to green settings. The left column specifies sensors appearing in the set of sensors selected for phase characterisation; whereas, the right column specifies the associated decision. For instance, when sensors selected from PCA for phase characterisation only include I/O related sensors, the decision associated is scaling the CPU frequency to its maximum.

2) Evaluation Platform: Our evaluation platform is a fifteen-node cluster set up on the French large-scale experimental platform called Grid5000 3 . Each node is an Intel Xeon X3440 with 4 cores and 16 GB of RAM. Available frequency steps for each core are: 2.53 GHz, 2.40 GHz, 2.27 GHz, 2.13 GHz, 2.00 GHz, 1.87 GHz, 1.73 GHz, 1.60 GHz, 1.47 GHz, 1.33 GHz and 1.20 GHz. In our experiments, low computational level always sets the CPU frequency to the lowest available frequency, which is 1.20 GHz, whereas high and medium computational levels set the CPU frequency to the highest available (2.53 GHz) and 2.00 GHz respectively. Each node uses its own hard drive which supports active, ready and standby states. Infiniband-20G is used for interconnecting nodes. The Linux kernel 2.6.35 is installed on each node where perf event is used to read the hardware monitoring counters. MPICH is used as MPI library. LU, SP, BT from NPB-3.3 [START_REF] Bailey | The nas parallel benchmarks[END_REF] and a real life application the Advance Research WFR (WRF-ARW) model are used for the experiments. Class C of the NAS benchmark is used (compiled with default options). WRF-ARW is a fully compressible conservative-form nonhydrostatic atmospheric model. It uses an explicit timesplitting integration technique to efficiently integrate the Euler equation. We monitored each node power usage on a per second basis using a power distribution unit.

C. Experimental Results: Analysis and Discussion

To evaluate our management policy, we consider 3 basic configurations of the monitored cluster: (a) on-demand configuration in which the Linux "on-demand" governor is enabled on all the nodes of the cluster; (b) performance configuration where each node's CPU frequency scaling governor is set to "performance"; (c) phase-detect configuration in which we detect phases identified using partial recognition and apply green capabilities accordingly.

Histograms of Figure 8 present the normalised average energy consumption of the overall cluster for each application under the three cluster configurations. Figure 9 shows their respective execution time. These figures show that our management policy (phase-detect) consumes in average 15% less energy than "performance" and "on-demand" while offering the same performance for our real-life application, WRF-ARW. For LU, BT and SP the average energy gain is between 3% and 6%, as we could not achieve as much energy savings as with WRF-AWR since many phases detected for these applications are compute intensive. From Figure 9, we notice a performance loss of less than 3% for LU and BT (performance is evaluated in terms of execution time). No need to mention that "on-demand" performs better than "performance" since all the applications we use are CPU intensive. Hence, the CPU load is always above 95% , a situation wherein the on-demand governor remains at the highest frequency.

We do better than a Linux governor because Linux ondemand governor will not scale the CPU frequency down unless the system load decreases below a threshold. The problem at this point is that the CPU load generally remains very high for memory intensive workloads/phases that do not require the full computational power. Thus, the energy gain made with these applications and particularly with WRF-AWR results in slowing the CPU frequency down when a node is suspected to be in a memory-bound phase.

Above results demonstrate the effectiveness of our system's energy management scheme based on phase detection and partial recognition. However, there are several underlying issues that need to be addressed in the future. Some of . Phase detection and partial recognition guided green leverage results: the chart shows average performance loss incurred by our management policy with respect to "performance" and "on-demand".

these issues are highlighted in this paragraph. Performance counters are not the same on all architectures, which means one may need to find out which sensors are useful on his architecture. Despite this limitation, the system management scheme we introduce here is a general purpose model and can be easily scaled to integrate previous research efforts.

V. CONCLUSION

In this paper, we present an approach to optimise the energy consumed by high performance computing systems at runtime. We propose to detect phases of an HPC system based on its workload. Then, by using a partial phase recognition technique, we predict the green capabilities than should be applied to save energy. This detection and improvements are achieved during runtime without knowing any details about the application. A software prototype has been implemented and validated. Experimental results reveal the effectiveness of our phase detection approach under real-life workload and benchmarks. Comparison with baseline unmanaged execution shows that our partial phase recognition algorithm can guide the use of green capabilities and save up to 15% of energy with less than 1% performance degradation. Future work includes integrating workload consolidation along with task migration as core functionalities of our prototype.

Figure 1 .

 1 Figure 1. Example of phase characterization based on load

Figure 2 .

 2 Figure 2. Matrix of distance (manhattan) between execution vectors, where the diagonal represent the execution time line. The darker, the more similar.

Figure 3 .

 3 Figure 3. Decision process.

Figure 4 .

 4 Figure 4. Global coordinator overview.

Figure 5 .

 5 Figure 5. Online phase detection with a 5% threshold for WRF-ARW.

Figure 7 .

 7 Figure 7. Offline phase detection with a 5% threshold for WRF-ARW (the matrix corresponds to half of the execution time). The diagonal represents the execution time-line from the upper left corner to the lower right corner.

Figure 8 .

 8 Figure 8. Phase detection and partial recognition guided leverage results: the chart shows the average energy consumption by applications under the three system's configurations.

Figure 9

 9 Figure 9. Phase detection and partial recognition guided green leverage results: the chart shows average performance loss incurred by our management policy with respect to "performance" and "on-demand".

Table I TRANSLATION

 I OF PHASE TO GREEN SETTINGS (IO RELATED SENSORS INCLUDES NETWORK AND DISK ACTIVITIES).

	Sensors selected from PCA Decisions
	for phase characterisation
	bus cycles &	CPU frequency sets to its maximum
	cache misses &	spin down the disk
	IO related sensors	network speed scaled down
	no IO related sensor	CPU frequency sets to its lowest
		network speed scaled up
	instructions &	CPU frequency set to its minimum
	last level cache misses (llc)	network speed scaled up
	instructions or llc &	CPU frequency sets to its average value
	IO related sensors	network speed scaled down
		spin down the disk
	IO related sensors	CPU frequency sets to its maximum

http://www.top500.org

http://www.green500.org

ACKNOWLEDGMENT

This work is supported by the INRIA large scale initiative Hemera focused on "developing large scale parallel and distributed experiments". Some experiments of this article were performed on the Grid5000 platform, an initiative from the French Ministry of Research through the ACI GRID incentive action, INRIA, CNRS and RENATER and other contributing partners (http://www.grid5000.fr). Authors would also like to thank Marcos Dias de Asuncao for his help in improving the readability of this paper.