Julien Bigot
email: julien.bigot@cea.fr

Hélène Coullon
email: helene.coullon@inria.fr

Christian Pérez
email: christian.perez@inria.fr

From DSL to HPC Component-Based Runtime: A Multi-Stencil DSL Case Study

High performance architectures evolve continuously to be more powerful. Such architectures also usually become more difficult to use efficiently. As a scientist is not a low level and high performance programming expert, Domain Specific Languages (DSLs) are a promising solution to automatically and efficiently write high performance codes. However, if DSLs ease programming for scientists, maintainability and portability issues are transferred from scientists to DSL designers. This paper deals with an approach to improve maintainability and programming productivity of DSLs through the generation of a component-based parallel runtime. To study it, the paper presents a DSL for multi-stencil programs, that is evaluated on a real-case of shallow water equations.

INTRODUCTION

Programming complex scientific applications on high performance architectures, such as clusters of multi-core nodes containing accelerators, requires for a scientist to be an expert in his own domain, but also in low level parallel and high performance programming. It can take months to years for a scientist to learn HPC programming, as it is not his main activity. As hardware continue to evolve to become more powerful but also usually more complex and difficult to use efficiently, the learning is never over. Domain Specific Languages (DSLs) are a promising solution to hide intricate details of HPC programming from non experts. By being specific to a well-defined set of problems, those languages enable the generation of well parallelized and optimized code for various architectures while requiring only a limited amount of information from the developer; part of the knowledge is directly embedded in the compiler. A critical question however concerns the specificity of DSLs. Actually, the less specific the DSL is, the more information its users must provide, and the less automatic optimizations can be embedded in the compiler. Conversely, the more specific the DSL is, the less applications can use it, and the more difficult it is to amortize the implementation cost of the DSL (language, compiler, and runtime). Programming complexity, as well as non-maintainability and non-portability often seem to be simply transferred from scientists to DSL designers.

For this reason, solutions to ease the design and implementation of DSLs have recently appeared [START_REF] Fernández | A data flow language to develop high performance computing dsls[END_REF]. The main idea of those solutions is to propose a single DSL or a framework to write new DSLs. In this paper, we study a first contribution to another approach to improve maintainability and portability in DSLs. This approach proposes a transformation from a DSL to a component-based runtime. Componentbased Software Engineering has proved many times good properties for code re-use, separation of concerns, maintainability and productivity of codes. For this reason, we think that combining DSLs and component-based runtimes could enable those properties, by inheritance, in DSLs. This paper proposes a first evaluation of this approach though a usecase study that is a DSL for multi-stencil programs, i.e., a subclass of numerically solved partial differential equations (PDEs), and its compiler.

Many solutions to automatically optimize and parallelize numerical simulations have been proposed in the literature, either as specific frameworks [START_REF] Balay | Efficient Management of Parallelism in Object Oriented Numerical Software Libraries[END_REF][START_REF] Coullon | The SIPSim implicit parallelism model and the SkelGIS library[END_REF] or as DSLs [START_REF] Christen | PATUS: A Code Generation and Autotuning Framework for Parallel Iterative Stencil Computations on Modern Microarchitectures[END_REF][START_REF] Devito | SC '11[END_REF][START_REF] Giles | Performance Analysis of the OP2 Framework on Many-core Architectures[END_REF][START_REF] Tang | The pochoir stencil compiler[END_REF]. Most of the time, DSLs focus on the parallelization and the optimization of a single numerical computation, also called a stencil kernel. However, a real case numerical simulation (such as fluid dynamics, magneto-hydrodynamics, molecular dynamics etc.) is most of the time not composed of a single stencil kernel, but of a set of stencil kernels and a set of additional local auxiliary computations [START_REF] Ragan-Kelley | Halide: A language and compiler for optimizing parallelism, locality, and recomputation in image processing pipelines[END_REF]. The DSL proposed in this paper, called Multi-Stencil Language or MSL, is a DSL to generate a parallel component-based structure of a numerical simulation, that we call a multi-stencil program. As discussed in the related work, this work is complementary to the optimization and parallelization of a single stencil kernel and is independent from implementation choices (as for example the form of the studied mesh).

The rest of this paper is organized as follows. The concepts of stencil kernels and multi-stencil programs are formalized in Section 2. Section 3 gives an overview of the proposed solution while Section 4 details the MSL language and its compiler which automatically generates the parallel computation part. Section 5 introduces component models and describes the proposed component-based runtime of MSL. Section 6 focuses on a real-case simulation, solving the Shallow-water equations, and analyzes preliminary performance and usability results. Sections 7 states on related work, while Section 8 concludes the paper.

STENCIL AND MULTI-STENCIL 2.1 Stencil Kernels

To numerically solve a set of PDEs, iterative methods (finite difference, finite volume, finite element methods etc.) are frequently used to approximate the solution through a discretized (step by step) process. Thus, the continuous time and space domains are discretized so that a set of numerical computations are iteratively (time discretization) applied onto a mesh (space discretization). In other words, in a mesh-based numerical simulation, the PDEs are transformed to a set of numerical computations applied at each time step on elements of the discretized space domain (the mesh). This paper focuses on one specific category of such numerical schemes based on stencil kernels [START_REF] Tang | The pochoir stencil compiler[END_REF], also called in applied mathematics explicit numerical schemes. This section introduces some definitions related to stencil kernels.

A mesh is the discretization of a physical domain. It is a connected undirected graph without bridges (an edge is a bridge if its removal results in two disconnected graphs), where nodes and edges are linked to form cells (closure). Most of the time, nodes of a mesh are linked to coordinates in the real space domain. An example of structured mesh is illustrated in Figure 1a. Each cell contains four nodes and four edges. Mesh entities are elements of the mesh, such as the center of cells, edges or nodes. In Figure 1a, two mesh entities are illustrated, the center of the cells (named Cells, in red) and edges on the horizontal axis (named Edgex, in blue on left and right of each cell). A data is a simulated quantity on which computations are performed. Each data is mapped onto mesh entities. For example, in Figure 1b, data A and B are mapped onto the center of cells, while, in Figure 1c, C is mapped onto the edges.

A stencil kernel computes the value of one data or a subpart of it (the kernel computation domain) using a numerical expression which takes as input one or more data. For example, the stencil kernel illustrated in Figure 1b computes A using B, while the one in Figure 1c computes A using C. Thus, a stencil kernel is defined by its numerical expression, a set of input data (only one in the examples), and its unique output data, the result. The computation domain is a subset of the mesh entities on which the output data is mapped. For example, in the first stencil, the computation of A is performed on elements represented with full lines, while dotted elements are not computed. On the other hand, on the second computation, the computation domain of A contains all the mesh entities on which it is mapped. This computation domain defines the elements over which the space loop iterates.

The numerical expression of a stencil kernel has the particularity to compute each element of the result independently using some elements of the inputs in a given neighborhood. Accessed elements form a neighborhood of the output known as the stencil shape. For example, the stencil shape of the first computation in Figure 1b contains direct neighbors on the right, left, top and bottom. Sometimes, the neighbor-

Mesh

Cells Edgex (a) Mesh and mesh domains. hood can also access different mesh entities, as for example in Figure 1c. Actually, in this computation, the neighborhood contains edges on the left and on the right of a cell.

As an example, the numerical expression of the first example could be:

A(x, y) = B(x+1, y)+B(x-1, y)+B(x, y +1)+B(x, y -1),

and the numerical expression of the second kernel could be:

A(x, y) = A(x, y) + C(x1, y1) + C(x1 + 1, y1).
As it can be seen from these definitions, stencil kernels implementations share many properties from an algorithmic point of view that can be used to apply well known optimization and parallelization strategies. Many solutions (languages or frameworks) thus propose to ease their programming by producing optimized and parallelized code from a simple description of the local computation to apply on each element. These are further discussed in Sections 7.

While stencil kernels have been studied a lot, the formalization of real overall mesh-based numerical simulations is poorly studied and not exactly for the same kind of appication [START_REF] Ragan-Kelley | Halide: A language and compiler for optimizing parallelism, locality, and recomputation in image processing pipelines[END_REF]. Actually, paying attention to complex numerical simulations, it appears that most of them are composed of more than one stencil kernel, with one or more stencil shapes and of additional local computations. For example, we would like to formalize and parallelize a numerical simulation which chains the stencil kernels of Figures 1b and1c. The next section formalizes concepts of stencil kernel, local kernel and multi-stencil program.

A Multi-Stencil Formalization

This section introduces some formal definitions used in the rest of the paper. Let ∆ be the set of data of the simulation A stencil kernel s is defined as the quadruplet:

s(R, w, exp, d), (1)
where R is a set of pairs (r, n), with r ∈ ∆ a data read by the computation, and n the stencil shape (neighborhood) used to access r. The data written by the kernel is denoted w ∈ ∆. exp is the numerical expression of the stencil kernel.

Finally the numerical expression is applied on the computation domain d.

For example, in Figure 1b, assuming the computation domain (full lines) is dc1 and the stencil shape is n1, the stencil kernel can be defined as:

R : {(B, n1)}, w : A, d : dc1, exp : A(x, y) = B(x+1, y)+B(x-1, y)+B(x, y+1)+B(x, y-1).

On the other hand, in the example of Figure 1c, assuming the computation domain is dc2 and the stencil shape is n2, the stencil kernel is defined as:

R : {(C, n2), (A, local)}, w : A, d : dc2, exp : A(x, y) = A(x, y) + C(x1, y1) + C(x1 + 1, y1).
One can notice that the input data A is associated to the stencil shape local. This means that the stencil shape applied on A is limited to the element with the exact same coordinate as the result element.

The specificities of kernels whose stencil shapes are all local makes it possible to apply specific optimizations. We therefore provide a specific definition in this case. A local (or auxiliary) kernel l is defined as the quadruplet:

l(R l , w, exp, d), (2)
where R l is the set of input data locally accessed, and other parameters are the same as for normal stencil kernels.

We define a multi-stencil program as a sextuplet:

MSP(T, M, E, D, ∆, Γ), ( 3
)
where T is the set of time iteration to run the simulation, M is the mesh of the simulation, E is the set of mesh entities, D is the set of computation domains used for computations, ∆ is the set of data of the simulation, each one mapped onto a mesh entity, and finally Γ is the set of computations. The set of computations Γ is an ordered list of stencil and local kernels. One can notice that this work, for now, is limited to a single mesh type for a given simulation.

MSL AND MSCAC OVERVIEW

The solution proposed in this paper is based on both a domain specific language to describe a multi-stencil program, i.e., a complete numerical simulation, and on a compiler which parses the language, applies transformations for automatic parallelization, and dumps the representation to a component-based runtime, ready-to-fill with the user code.

The proposed DSL, named MSL (Multi-Stencil Language), is an agnostic descriptive language for multi-stencil simulations. Agnostic means that the description of a numerical simulation does not depend on implementation choices as, for example, the type of mesh (structured or unstructured), or its associated interfaces (how to define a stencil shape, etc.). Descriptive, on the other hand, means that MSL is not made to handle the expression of numerical computations, but only the description of a multi-stencil simulation, in a coarse-grain fashion. A MSL description is made of six main sections that match the six tuples elements of the definition of a multi-stencil program of Equation 3.

As illustrated in Figure 2, the compiler, named MSCAC for Multi-Stencil Component Assembly Compiler, is composed of three different phases. The first phase parses the MSL description into an intermediate representation (IR). The second phase called MSC transforms the computation part (Γ) of the IR into a series-parallel tree decomposition [START_REF] Valdes | The recognition of series parallel digraphs[END_REF] to introduce parallelism in the overall numerical simulation. Finally, the third phase dumps the IR to an actual implementation of the parallel pattern of the simulation.

The parallelization techniques handled by MSC are data parallelism, where data is split among processors or cores while applying the same computation on every piece of data, and task parallelism where the program is transformed into a graph of tasks with dependencies. Different strategies can be used to schedule the tasks of this graph, i.e., to plan their execution in a valid order. The MSC compiler presented in this paper handles the transformation of the dependency graph to a static scheduling.

The final phase of compilation, the dump, uses the static scheduling of the tasks, created from Γ, and additional information on T , M, ∆, E and D. It generates a program in the back-end language which embeds the parallel pattern or skeleton of the overall simulation. Various models could be used as back-end, such as OpenMP [START_REF] Dagum | Openmp: an industry standard api for shared-memory programming[END_REF] for example. However this paper adopts a novel approach where the dump targets a component-based runtime. This approach makes it possible to leverage good properties of software component models such as increased code re-use, eased separation of concerns, and increased maintainability and productivity. These properties counter some known limitations of DSLs.

Thanks to the software components based approach, the generated code is made of four independent parts that loosely match the elements of the multi-stencil program definition. the impact on the other parts, as discussed in Section 5.

LANGUAGE AND PARALLELIZATION 4.1 Multi Stencil Language

As mentioned in the previous section, the Multi Stencil Language (MSL) is an agnostic descriptive language for multistencil simulations. Six main sections are required in a MSL description that match the six elements of the formal definition; they describe: 1) the mesh, 2) the mesh entities, 3) the computation domains with their inter-dependencies, 4) the data, 5) the time iterations, and finally 6) the computations. Figure 3 shows the grammar of MSL. Lines 1 to 9 define what is a MSL program with the different parts mentioned above. The remaining of this section describes these different parts of a MSL program and presents an example. One can notice that terminals are the string identifiers meshid, meshdom, compdomid, dataid, compid and neighborid. The terminal num is an integer terminal.

Mesh and Mesh Entities

As illustrated in the first line of the grammar in Figure 3, a mesh is simply defined by a string identifier. The different kind of mesh entities are described as a list of identifiers (Lines 2, 10 and 11 in Figure 3). Figure 4 gives an example of a mesh named cart used to represent a Cartesian mesh, and three kinds of mesh entities: cell, edgex, and edgey.

Computation Domains and Their Dependencies

A computation domain is a subset of a given kind of mesh entities, used for at least one computation. For each computation domain, the kind of mesh entities on which it is mapped is specified (lines 3-4 and 12-14 on Figure 3). By default, two computation domains of the same entities are handled as intersecting except when an explicit independence relation is mentioned in the independent section of the language (lines 5-6 and 15-17 in Figure 3). On lines 4 to 6 of Figure 4, computation domains are defined for each kind of mesh entities: allcell, alledgex and alledgey; they are used to represent the whole domain. On lines 7 and 8, two other domains of edgex are defined: part1edgex and part2edgex; they could be used to represent the boundaries of the domain, for example, or any other sub-part of the whole domain. It is the user's responsability to define the needed domains for his simulation. On line 10 the independence relation between part1edgex and part2edgex is specified.

Data and Time

A data element is a quantity to simulate. It is mapped on a given kind of mesh entities, as illustrated in lines 7 and 18-20 of Figure 3. The time section simply indicates a number of iterations to perform in the simulation. In the example of Figure 4 ten data elements are defined and 500 iterations (lines [START_REF] Fernández | A data flow language to develop high performance computing dsls[END_REF][START_REF] Ferrari | A new two-dimensional shallow water model including pressure effects and slow varying bottom topography[END_REF][START_REF] Gautier | Xkaapi: A runtime system for data-flow task programming on heterogeneous architectures[END_REF][START_REF] Giles | Performance Analysis of the OP2 Framework on Many-core Architectures[END_REF][START_REF] Mitchell | Creating minimal vertex series parallel graphs from directed acyclic graphs[END_REF][START_REF] Nieplocha | Advances, applications and performance of the global arrays shared memory programming toolkit[END_REF][START_REF]Corba component model 4.0 specification[END_REF][START_REF] Ragan-Kelley | Halide: A language and compiler for optimizing parallelism, locality, and recomputation in image processing pipelines[END_REF][START_REF] Richard | Towards application variability handling with component models: 3d-fft use case study[END_REF][START_REF] Schmitt | Exaslang: A domain-specific language for highly scalable multigrid solvers[END_REF][START_REF] Schoenmakers | A new algorithm for the recognition of series parallel graphs[END_REF][START_REF] Szyperski | Component Software: Beyond Object-Oriented Programming[END_REF].

Computation Description

The last part of the language contains the specifications of stencil and local kernels as defined in Section 2.2. This description does not contain a direct expression of the numerical expressions exp. Instead, the term exp of Equations 1 and 2 corresponds to an identifier that references an implementation in an external language (cf. Section 5). The description starts by the identifier of the computed data element (w in the formal description) with the computation domain between brackets (d). Then, after the equal sign, the kernel identifier (exp) is specified. Finally, between the parenthesis is the list of identifiers of the data elements read by the computation with their stencil shape between brackets (R). If the computation is local, no brackets appear (R l ).

In the example of Figure 4, nine computations are described. On line 24, the data element b is computed on the domain allcell by the kernel expression c_0 which reads data a without any neighborhood shape.

MSC Parallelization

The MSC parallelization is a subpart of the overall compiler MSCAC. It makes use of the ordered list of computations Γ, which can directly be extracted from the parser, to build a parallel representation of the computations of the overall multi-stencil program. This parallelization phase is divided in five different steps to transform Γ to a series-parallel tree decomposition [START_REF] Valdes | The recognition of series parallel digraphs[END_REF]. This sections describes these steps.

The Ordered List of Computations: Γ Γ is directly obtained from the parser. Actually, the list of computations in the MSL program is already ordered: Γ is a direct map of this list. In the example of Figure 4,

Γ = [c0, c1, c2, c3, c4, c5, c6, c7, c8].
The Synchronized Ordered List of Computations: Γsync Data parallelization splits data among processors and the same program is applied on each sub-part of the data. However, this parallelization technique requires additional synchronizations between processors. The required synchronizations can be automatically computed from the ordered list of computations Γ. A synchronization is needed each time a data read by a stencil computation has been written by a previous computation. In such a case, a computation is added before the stencil computation. This synchronization computation reads the data to synchronize, and write the same data. The computation domain of such a synchronization computation encompasses the whole mesh entities on which the data is declared. As a result, Γ is transformed to a synchronized ordered list of computations Γsync. For example, in Figure 4, the stencil computation c1 reads the data element b which has been written by c0.

For this reason the sublist [c0, c1] of Γ is transformed into [c0, sync1, c1] in Γsync. The new computation sync1 reads and writes b and is applied on edgex. As a result, a dependency is kept between c0 and sync1 and between sync1, and c1. The same transformation is performed for the two other stencils: c4 and c8. Thus, for the example of Figure 4, Γsync = [c0, sync1, c1, c2, c3, sync4, c4, c5, c6, c7, sync8, c8].

The Dependency Graph: Γ dep

From the synchronized ordered list of computations Γsync, a dependency graph Γ dep is built. A dependency exists between computations a and b (including synchronizations) if and only if a data element read by b has been written by a in Γsync with intersecting domains. Nodes of the dependency graph represent computations of Γsync, while edges are dependencies between them. The dependency graph Γ dep is a directed acyclic graph (DAG). For example, the dependency DAG Γ dep of the example of Figure 4 is presented in Figure 5. The Minimal Series-Parallel Graph: Γmsp Once a dependency graph is built, many solutions can be used to build a parallel application, as for example dynamic schedulers [START_REF] Augonnet | StarPU: a unified platform for task scheduling on heterogeneous multicore architectures[END_REF][START_REF] Gautier | Xkaapi: A runtime system for data-flow task programming on heterogeneous architectures[END_REF]. In this work a static scheduling of the dependency graph is built. To do so the dependency graph is transformed to a minimal series-parallel graph [START_REF] Valdes | The recognition of series parallel digraphs[END_REF]. As it has been shown in [START_REF] Valdes | The recognition of series parallel digraphs[END_REF], the transitive reduction of a DAG is a minimal series-parallel graph if and only if the forbidden shape called N-shape, illustrated in Figure 6a, is not found in the DAG. Thus, to transform Γ dep to the minimal seriesparallel graph Γmsp, an algorithm is applied to remove all Nshapes by an over-constraint [START_REF] Mitchell | Creating minimal vertex series parallel graphs from directed acyclic graphs[END_REF], as illustrated in Figure 6b. The Tree Decomposition: Γtsp Many works [START_REF] Schoenmakers | A new algorithm for the recognition of series parallel graphs[END_REF][START_REF] Valdes | The recognition of series parallel digraphs[END_REF] explains how to build a series-parallel tree decomposition from a minimal series-parallel graph. This transformation decomposes the minimal series-parallel graph as a tree where internal nodes are sequences, indicated by a S node, and parallel sections, indicated as a P node, while leaves of the tree are the nodes of Γmsp. For example, the series-parallel tree decomposition Γtsp of Figure 5 is illustrated in Figure 7.

Data parallelism optimizations: Γ data

The series-parallel tree decomposition Γtsp is used to build the hybrid parallelization of the simulation, but it can also be used to optimize the data parallelization. In fact, it first seems natural to create the data parallel structure directly from Γsync, as the program is parallelized from data decomposition and applying the same sequential tasks on each processor. However, as a computation kernel writes a single data w (to finely manage dependencies between tasks), the same computation domain might be iterated more than once, while some numerical expressions could be computed inside a joint space loop. For this reason, an optimized and merged ordered list of computation Γ data is built from Γtsp. Actually, in Γtsp, computations of a parallel subtree with a same computation domain can be computed in the same space loop, in any order. Moreover, in Γtsp, two consecutive computations of a sequence subtree with a same computation domain can be computed in the same space loop, if the computation order of the sequence is kept. For example, from Γtsp of Figure 7, as kernels c3 and c5 are executed inside a parallel section, and as their computation domains are the same (alledgey), a single component K can be written peforming both computations in the same domain loop.

S c0 sync1 c1 P c7 sync8 c8 S S c2 sync4 c4 c6 c3 c5

A COMPONENT BASED BACK-END

This section describes the projection of a MSL program to a component-based runtime, including the components used to implement the static scheduling of Γtsp. The section first gives an overview of component models, and especially the Low Level Components (L 2 C) used in this work.

Component-Based Software Engineering

Component-based software engineering (CBSE) is a domain of software engineering [START_REF] Szyperski | Component Software: Beyond Object-Oriented Programming[END_REF] which promotes code re-use, separations of concerns, and thus maintainability. An application is made of a set of component instances. A component is a black box that implements an independent functionality of the application, and which interacts with its environment only through well defined interfaces: its ports. A port can, for example, specify services provided or required by the component. With respect to high performance computing, some works have also shown that component models can achieve the needed level of performance, and scalability while also helping in application portability [START_REF] Allan | A component architecture for high-performance scientific computing[END_REF][START_REF] Bigot | Increasing reuse in component models through genericity[END_REF][START_REF] Richard | Towards application variability handling with component models: 3d-fft use case study[END_REF] Many component models exist, each of them with its own specificities. Well known component models include, for example, the CORBA Component Model (CCM) [START_REF]Corba component model 4.0 specification[END_REF], and the Grid Component Model (GCM) [START_REF] Baude | Gcm: a grid extension to fractal for autonomous distributed components[END_REF] for distributed computing, while the Common Component Architecture (CCA) [START_REF] Allan | A component architecture for high-performance scientific computing[END_REF], and Low Level Components (L 2 C) [START_REF] Bigot | A low level component model easing performance portability of hpc applications[END_REF] are HPC-oriented. This work makes use of L 2 C for the experiments. instances and connections between their ports. Such an assembly is described in LAD, an XML dialect.

L

MSP Component Runtime Overview

As described in Section 3, the compiler generates four independent parts in a component assembly in addition to a static part as represented in Figure 9: each rounded box represents one component instance, drawn with simple lines, or more, drawn with double lines. For example, the component Data is instantiated more than once. Driver is static and drives the whole application. The four other parts are generated based on the content of the different sections of the MSL program: Computations, from the parallelization of Γ (done by MSC as explained in Section 4.2), DDS from M and E, Data from ∆ and D, and Time from T .

DDS manages the structure and the partitioning of the mesh, a single instance is generated as MSL only supports a single mesh for now. The back-end presented in this paper uses the distributed data structure proposed in SkelGIS library [START_REF] Coullon | The SIPSim implicit parallelism model and the SkelGIS library[END_REF]. and Γ data , the computation components instances (denoted K), provided by the user, are connected to Data. As shown in Figure 10, a K component provides a port to launch the computation and exposes a use port for each data element manipulated in the computation; a star is added on the use port to Data to denote it.

DDS is used by

A complete assembly is obtained by replicating the component assembly of Figure 9 on available processors (or cores).

Connections is done with MPI ports of Data components.

As a result of this approach, each part of the generated code is rather independent. Changing the main time loop to use a convergence criteria rather than a fixed number of iterations would only require changing the Time component. Similarly, changing the approach used for scheduling would only impact the code generated in Computations. The technology used for data parallelism can also be changed by replacing the DDS and Data components, and by using the adequate interfaces in K. As a matter of fact, we have started evaluating an alternative based on Global Arrays [START_REF] Nieplocha | Advances, applications and performance of the global arrays shared memory programming toolkit[END_REF] instead of SkelGIS.

Control Components and Dump

The series-parallel tree decomposition Γtsp represents a static scheduling of the simulation. Γ data , on the other hand, is an optimized merged ordered list of computations. Both are dumped to components by generating an assembly that exactly matches their structure. We introduce what we call control components to represent all nodes of Γtsp or Γ data . These components can be used for any simulation, which increases code reuse between simulations. A control component exposes a single provide port containing a control-only method, without any parameter. Three control components have been needed: SEQ, PAR, and SYNC. They are represented in Figure 11.

Sequence component (SEQ) It is the direct representation of a sequence node of Γtsp or Γ data . This component sequentially calls an ordered list of other components. It exposes an ordered use-multiple port to be connected to the components to call in sequence.

Parallel component (PAR)

It is the direct representation of a parallel node of Γtsp (not Γ data ). This component simultaneously calls a set of other components. It exposes a use-multiple port to be connected to the components to call in parallel.

Synchronization component (SYNC)

It is the direct representation to an update leaf of Γtsp or Γ data . This component calls the synchronization of a given data. It exposes a use port to be connected to the data to be updated.

Then, from Γtsp or Γ data , and using the control components described above, a direct dump can be done to a component assembly for a data or hybrid (data and task) parallelization of a simulation. Figure 12 displays the assembly part corresponding to Γtsp of Figure 7. In this figure, the ports linked to data (use and use-multiple ports of SYNC and K) are represented but are not connected for readability. Moreover, each computation component is an instance of a specific K component using the identification name of the computation of an MSL instance. Thus, this figure represents what the component Computations of Figure 9 contains if an hybrid parallelization is performed. On the other hand, if a data parallelization is performed, a single SEQ component is used by Time, which invokes in the appropriate order the computation kernels and synchronizations. Among those kernels, a single merged kernel would represent c3 and c5.

REAL-CASE EVALUATION

Navier-Stokes equations is a well known set of PDEs in fluid dynamics to simulate a flow evolution in time. At the University of Orléans, France, the MAPMO laboratory works on a software, called FullSWOF1 , which solves the Shallow water equations obtained from the three dimensional Navier-Stokes equations, by averaging on the vertical direction [START_REF] Ferrari | A new two-dimensional shallow water model including pressure effects and slow varying bottom topography[END_REF]. Those equations are solved using a two-dimensional Cartesian discretization of the space domain, and a finite volume numerical methods more described in [START_REF] Coullon | The SIPSim implicit parallelism model and the SkelGIS library[END_REF]. We have developed a MSL version of FullSWOF that contains 3 mesh entities, 7 computation domains, 48 data and 98 computations. These computations are made of 32 stencil kernels and 66 local kernels.

Compiler Evaluation

The series-parallel tree decomposition Γtsp of this simulation, extracted by the MSC transformation, is composed of 17 sequence nodes and 18 parallel nodes. Figure 13 represents for a given level of parallelism, i.e., the number of tasks to perform concurrently, the number of time this level is observed in the final component assembly. One can notice that the level of task parallelism extracted from the Shallow water equations is limited by two sequential parts in the application (level 1). As a level of 16 parallel tasks is reached two times, and also five times for the fourth level, sequential restrictions could be amortized. Moreover, each identified task can itself be parallelized, introducing data parallelism on loops. Thus, it is rather clear that the task parallelization technique should be used to improve the data parallelism when reaching its limits, but not to use alone. Moreover, as the level of parallelism in the application is heterogenous, the number of threads to launch for task parallelism, and the number of cores to keep for data parallelism is difficult to choose. One can notice that the transformation of Γ dep to a minimal series-parallel graph is the longest step of MSC, because of the removal of the forbidden N-shapes in the graph. Thus, the fact that many N-shape are removed in Γ dep shows that the creation of a static schedule of tasks may not be the best solution for complex simulations. It has to be noticed, however, that the all computation of Γtsp is still usefull to detect how space loops can be merged in Γ data .

Preliminary Performance Evaluation

To evaluate performance of the generated component-based parallel structure, we have proceeded as follows. Our current implementation of the MSCAC compiler generates a backend code using the SkelGIS distributed data structure as an underlying library. For this reason, our evaluation compare the shallow water equations first parallelized with a pure SkelGIS code (data structures, applicators and interfaces of SkelGIS [START_REF] Coullon | The SIPSim implicit parallelism model and the SkelGIS library[END_REF]), and second parallelized with MSL (which uses the SkelGIS data structure only). Moreover, as SkelGIS library handles data parallelization for distributed memory architectures, we have limited for now our evaluations to the merged data parallelization of MSCAC (dumped from Γ data ).

Evaluations have been performed on the cluster Edel of Grid'5000 2 for which each node is equiped with two 8-cores Intel Xeon E5520, 24GB of memory, and is interconnected with an Infiniband 40G network. On Edel, three different experiments have been performed that differ in the size of 2 https://www.grid5000.fr/ the Cartesian mesh and the number of time iterations: (1) 5, 000 × 5, 000, 500 iterations; (2) 10, 000 × 10, 000, 500 iterations; (3) 10, 000 × 10, 000, 2, 000 iterations. Figure 15 reports the results of the three experiments. The first row of graphs is the number of iterations computed per second for various number of cores. The second row of graphs represents the execution time for various number of cores, using a logarithmic scale.

One can notice that, for the three different experiments, execution times of MSL, which itself uses the SkelGIS distributed data structure, are improved compared to a pure SkelGIS parallelization. This result can be explained by the automatic generation of the parallel component-based runtime. Actually, in a pure SkelGIS code, communications between processes are hidden from the user through specific interfaces called applicators [START_REF] Coullon | The SIPSim implicit parallelism model and the SkelGIS library[END_REF]. Using MSL, on the other hand, communications are directly managed by SY N C components, which are automatically placed in the componentbased runtime by MSCAC. Thus, no overhead is introduced to hide communications. In addition to this, the scalability slope, shown by the number of iterations per second, is close to SkelGIS. This is an important result as the execution time of MSL is at the same time improved compared to SkelGIS. Actually, the less the execution time is, the less the scalability capacity is. Thus, these results confirm the result of [START_REF] Bigot | A low level component model easing performance portability of hpc applications[END_REF] that using a component-based runtime does not damage performance of the back-end code. Finally, we notice that SkelGIS has proved its scalability on the Shallow water equations compared to an MPI parallelization [START_REF] Coullon | The SIPSim implicit parallelism model and the SkelGIS library[END_REF], that makes our performance results relevant.

RELATED WORK

Many specific solutions exist to ease parallel programming of numerical simulations, as for example PATUS [START_REF] Christen | PATUS: A Code Generation and Autotuning Framework for Parallel Iterative Stencil Computations on Modern Microarchitectures[END_REF] and Pochoir [START_REF] Tang | The pochoir stencil compiler[END_REF] which propose a simple description to write stencil codes for general structured meshes, and OP2 [START_REF] Giles | Performance Analysis of the OP2 Framework on Many-core Architectures[END_REF] and Liszt [START_REF] Devito | SC '11[END_REF] for general unstructured meshes. Those solutions can be considered as stencil compilers. Actually they produce optimized (cache tiling, etc.) or parallel (CUDA, OpenCL, etc.) codes for a single stencil computation. Recently, domain specific languages for very specific numerical methods such as ExaSlang [START_REF] Schmitt | Exaslang: A domain-specific language for highly scalable multigrid solvers[END_REF] for multigrid solvers, have also been proposed.All those solutions let the user imple- ment their own stencil codes in a sequential programming style, and produce high performance codes for shared memory architectures and sometimes distributed memory architectures. Most stencil compilers, though, handle the parallelization or the optimization of a single stencil kernel, considering that it represents the main computation time of numerical simulations. However, most real case numerical simulation are composed of more than one stencil kernel, involving one or more stencil shapes, and additional local kernels, as explained in Section 2.2.

Other solutions, as SkelGIS library [START_REF] Coullon | The SIPSim implicit parallelism model and the SkelGIS library[END_REF] or Global Arrays [START_REF] Nieplocha | Advances, applications and performance of the global arrays shared memory programming toolkit[END_REF] propose more or less specialized distributed data structures (specific meshes for SkelGIS, general arrays for Global Arrays), which simplifies or totally hide communications and synchronizations between processes.

Halide [START_REF] Ragan-Kelley | Halide: A language and compiler for optimizing parallelism, locality, and recomputation in image processing pipelines[END_REF] is a DSL which focusses on the automatic parallelization and optimization of stencil pipelines in image processing. Thus, this work is the closest to MSL. However, at least one important difference can be noticed: the specific domain it is applied on. Actually, image processing is applied onto 2D or 3D grids, while explicit numerical schemes to solve PDEs can be applied on many different meshes (grids, but also unstructured or hybrid meshes). This main difference is the reason why MSL is a descriptive and agnostic language, and why it does not target specific optimizations, as tiling for example.

The work presented in this paper takes place at a different abstraction level and can be seen as complementary to stencil compilers or distributed data structures. The MSCAC compiler only needs the agnostic description of the overall simulation to generate a parallel pattern of the program. This parallel pattern then needs to be filled with implementation parameters and computation codes. Those computation codes could be generated by stencil compilers or could be written using existing distributed data structures (as SkelGIS in this paper). As a result, MSL is not a new contribution to stencil compilers, to tiling optimizations, or to distributed data structures, but proposes a coarse-grain automatic parallelization of an overall multi-stencil program. Finally, MSL is a case study to show that component-based runtimes could improve maintainability adn flexibility of DSLs with a low impact on performance.

CONCLUSION

This paper has presented the domain specific language MSL and its compiler MSCAC designed to produce a parallel (data and hybrid) component-based runtime of an overall multi-stencil program, i.e., a mesh-based numerical simulation reduced to explicit schemes. After details on the language and its compiler, MSL has been evaluated by the description of a real case numerical simulation: the shallow water equations. The component-based data parallelization of the simulation has been compared to a pure SkelGIS parallelization, and has shown improved execution times as well as a promising scalability. Those results demonstrate that component-based runtimes may be relevant back-end codes for DSLs as they do not introduce performance damage. Moreover, components bring software engineering benefits such as separation of concerns, code re-use, improving maintainability.

Although MSL is a promising case study from DSLs to component-based runtimes, many works in progress aims at improving this first contribution. First, to more clearly show the improvement of DSLs maintainability using componentbased back-end, an alternative DDS component is under study, using Global Arrays [START_REF] Nieplocha | Advances, applications and performance of the global arrays shared memory programming toolkit[END_REF]. In addition to this, alternatives for the Computations component, computed by MSC, are under study such as a dump to a pure OpenMP [START_REF] Dagum | Openmp: an industry standard api for shared-memory programming[END_REF] code or the use of dynamic schedulers [START_REF] Augonnet | StarPU: a unified platform for task scheduling on heterogeneous multicore architectures[END_REF][START_REF] Gautier | Xkaapi: A runtime system for data-flow task programming on heterogeneous architectures[END_REF].

  2-neighborhood 2-entity stencil.

Figure 1 :

 1 Figure 1: (a) a Cartesian mesh and two kind of mesh entities, (b) an example of stencil kernel on cells, (c) an example of stencil kernel on two different entities of the mesh.

Figure 2 :

 2 Figure 2: The three phases of the MSCAC compiler: i) the parser of the MSL description; ii) the MSC parallelization; iii) the dump to the Component-Based Runtime (CBR).

  Those are the components implementing: a) the computations scheduling (Γ), b) the data structures (M and E), c) data of the simulation (∆ and D), and d) the time loop (T ). The actual implementation of the numerical expressions are provided by the user and used in the computations scheduling. This component-based parallel structure enables to change the technology used for each part while limiting

Figure 3 :

 3 Figure 3: Grammar of the Multi-Stencil Language.

Figure 4 :

 4 Figure 4: Example of a MSL program.

Figure 5 :

 5 Figure 5: Γ dep of the example of Figure 4.

  Over-constraint on the forbidden N shape.

Figure 6 :

 6 Figure 6: Forbidden N subgraph shape for a DAG to be minimal series-parallel.

Figure 7 :

 7 Figure 7: Γtsp of Figure 5. Nodes S represent sequences, nodes P represent parallel sections.

Figure 8 :Figure 9 :

 89 Figure 9: MSCAC Component-based Runtime Overview.

*Figure 10 :

 10 Figure 10: A kernel component K.

Figure 11 :

 11 Figure 11: The three control components.

Figure 14 illustratesFigure 12 : 7 Figure 13 :Figure 14 :

 141271314 Figure 14 illustrates the execution time for each step of the MSC transformation for an overall execution time of ten seconds. Execution times have been computed on a laptop with a bi-core Intel Core i5 1.4 GHz, and 8 GB of DDR3.

  10k × 10k mesh, 2k iterations

  5k × 5k mesh, 500 iterations

  10k × 10k mesh, 500 iterations

Figure 15 :

 15 Figure15:The first line represents, for each experiment, the number of iteration computed per second for a given number of cores. The second line represents, for each experiment, the execution time for a given number of cores, using a logarithmic scale.

  2 C L 2 C[START_REF] Bigot | A low level component model easing performance portability of hpc applications[END_REF] is a minimalist C++ based HPC-oriented component model where a component extends the concept of class. The services offered by the components are specified trough provide ports, those used either by use ports for a single service instance, or use -multiple ports for multiple service instances. Services are specified as C++ interfaces. L 2 C also offers M P I ports that enable components to share MPI communicators. Components can also have attribute ports to be configured. As illustrated in Figures8, aprovide port is represented by a white circle, a use port with a black circle, a use -multiple port by a black circle with a white m in it. MPI port are connected with a black rectangle. A L

2

C-based application is a static assembly of components

http://www.univ-orleans.fr/mapmo/soft/FullSWOF/

ACKNOWLEDGMENTS

This work has partially been supported by the PIA ELCI project of the French FSN. Experiments presented in this paper were carried out using the Grid'5000 testbed, supported by a scientific interest group hosted by Inria and including CNRS, RENATER and several Universities as well as other organizations (see https://www.grid5000.fr).


PhD Marta Rava
email: mrava@cnio.es

PhD Nicole Le Moual

Xavier Dumont

Stefano Guerra

Valerie Siroux

Bénédicte Jacquemin

Francine Kauffmann

Alfred Bernard

Rachel Nadif

Benedicte Jacquemin

Serum club cell protein 16 is associated with asymptomatic airway responsiveness in adults: Findings from the French epidemiological study on the genetics and environment of asthma

Keywords: Uteroglobin, asthma, lung function test, airway hyper-reactivity, biological marker

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

INTRODUCTION

Club cell protein (CC-16, also called CC10, SCGB1A1 or CCSP) is mainly secreted by the club cells along the tracheobronchial tree, diffuses passively into the bloodstream, and is finally rapidly eliminated by glomerular filtration 1 .

While the exact function of CC-16 is not known, it is suggested to be a protective mediator in the airway inflammatory process 2 and to have protective effects from oxidative stress on the respiratory tract 3,4 . Serum CC-16 level is considered as a sensitive biomarker of airways epithelium integrity, transiently increasing following acute exposure to pulmonary irritants and decreasing after chronic exposures associated with an impairment or reduction of club cells. Overall the assay of CC-16 is a sensitive test to allow the detection of early damage to the respiratory epithelium 3 and because of its presumed relationship to inflammation, CC-16 has gained increasing interest as a biological marker in the study of chronic pulmonary diseases 5,6 .

Asthma is a chronic inflammatory disorder of the airways characterised by increased airway responsiveness (AR) and associated with accelerate decline in lung function 7 . Although AR is a key characteristic of asthma 8 it may occur without asthma-like symptoms or previous diagnosis of asthma, which is called asymptomatic AR. 9 .

Findings on the association between CC-16 and asthma are conflicting. No association with asthma was observed in adults from a large cross-sectional study 10 , while a negative association was observed in participants with asthma in clinical studies 11,12 . Only few large epidemiological studies evaluated the association between CC-16 level and lung function, providing evidence of a negative association 10,[13][14][15] . A negative association between serum CC-16 level and AR was described only in one cross-sectional study 14 , and no study evaluated the specific association with asymptomatic AR. Overall there is a lack of large epidemiological studies considering together the relations of serum CC-16 with asthma, lung function and AR. Our hypothesis is that serum level of CC-16 decreases in participants with asthma, impaired lung function and airway responsiveness. The aim of the present paper is to study the associations between serum CC-16 levels and asthma, impaired lung function and airway responsiveness in the French epidemiological study on the genetics and environment of asthma (EGEA) whom participants with or without asthma have been extensively characterized regarding respiratory phenotypes and lung function tests.

.

METHODS

Study population

EGEA is a cohort study based on an initial group of asthma cases, their first-degree relatives, and controls (first survey, n=2047) 16 . At the second survey 17 , detailed phenotyping was performed with lung function tests including methacholine challenge (see Figure 1 and the online supplemental population and protocol section). This analysis is based on participants who were adults at the second survey (≥ 16 years old, n=1570 adults) with available data on asthma, serum CC-16 level and who did not smoke 1 hour before the examination (N=1298).

Biological phenotypes

Serum CC-16 level was determined with latex immunoassay using a rabbit anti-CC-16 antibody (Dakopatts, Glostrup, Denmark) 18 and purified CC-16 protein as standards (see the online supplemental biological phenotypes section).

Serum creatinine was determined with standard traceable IDMS (isotope dilution mass spectrometry traceable) assay and its clearance was calculated by the CKD-EPI (Chronic Kidney Disease Epidemiology Collaboration) formula 19 which depends on sex and age of participants.

Clinical phenotypes

Asthma and asthma related phenotypes

Current asthma referred to the report of asthma attacks or asthma treatment or asthma-like symptoms in the past 12 months 20 . Participants without current asthma and that never reported ever asthma were considered as controls. Based on five asthma symptoms over the past 12 months (wheezy breathlessness, woken up by chest tightness and by an attack of shortness of breath, attack of shortness of breath at rest and after exercise), the asthma symptom score was computed, as previously proposed 21,22 .

Asthma control was assessed from survey questions and lung function, closely adapted from the 2006-2009 Global Initiative for Asthma guidelines 2006, as previously described in 20 (see online supplemental population and protocol section for more details).

Pulmonary function

Pulmonary function test was performed using a standardized protocol with similar equipment across centres according to the ATS/ERS guidelines. FEV 1 and FVC percent predicted values were based on Stanojevic et al. reference equations 23 .

Airway responsiveness

Methacholine bronchial challenge test was performed in participants with no medical contraindication, unless baseline FEV 1 <80% predicted or FEV 1 post dilution was lower than 90% of the baseline FEV 1, following the same protocol the European Community Respiratory Health Survey (ECRHS) 24 . Airway hyper-responsiveness (AHR) was defined as a 20% fall in FEV 1 from the highest FEV 1 postdiluent (PD 20 ) during the test with an accumulated dose of methacholine. "High AHR" refers to a PD 20 related to an accumulated methacholine dose of 1 mg (PD 20 ≤1mg) and "any AHR" refers to a PD 20 ≤4 mg.

The degree of airway responsiveness was quantified by a transformed dose-response slope obtained regressing the percentage fall in FEV 1 on log10 methacholine dose 25 . Low values of slope correspond to high degree of bronchial responsiveness. We referred to participants without asthma but with airway responsiveness as "asymptomatic AR".

Ethics Statement

Ethical approval was obtained from the relevant institutional review board committees at Cochin Port-Royal Hospital and Necker-Enfants Malades Hospital, Paris. Written informed consent was signed by kin or guardians of the minors/children and all other adult participants.

Statistical analyses

Since CC-16 level strongly depends on glomerular filtration rate (eGFR) level (Table 1), we obtained eGFR-adjusted residuals by regressing log-transformed CC-16 on eGFR and used the standardized residuals (CC-16 z-score, Mean= 0, SD=1) in the subsequent analysis as independent variable. Association between CC-16 and asthma was estimated with a logistic regression model. Since the expression of CC-16 is known to be inducible by corticosteroids 26 , association with asthma was repeated after excluding participants who used inhaled corticosteroids during the previous 12 months (n=265). To account for the specific design of the EGEA study analyses on lung function and AR were conducted separately in participant with and without asthma and estimates were obtained with linear or logistic regression model, where appropriate. All parameters were estimated with generalized estimating equations (GEE) with an exchangeable working correlation to account for the potential clustering within the families. Estimates were adjusted for sex, smoking status and pack-years, age and height when needed. Since CC-16 has shown a time-dependent daily variation, 27 blood sampling time (range 7 -20h) was included as covariate. Further, to test whether the association between CC-16 and AR was independent from the bronchial calibre of participants, the level of FEV 1 % predicted was added as a covariate in the multivariate model. Finally, since tobacco smoke is associated with decreased level of CC-16, as well as with impaired lung function and increased AR, the analyses were stratified by smoking status. All statistical analyses were performed using R version 2.14. 28

RESULTS

Characteristics of the sample (Table 1)

The study population included 1298 participants with a mean age of 43 years (range 16-80 years). Lung function test was performed by 1278 participants: among them 416 did not perform the methacholine test because of poor lung function (n=113) or other reasons (n=303) and data on AR were available for 862 participants (see online result section for more details).

CC-16 Serum level

Overall, serum CC-16 geometric mean (GM) level was 12.4 µg/L, ranging from 2.2 to 70.6µg/L (Table 1). CC-16 level was higher in men, increased with age, decreased in smokers and with pack-years smoked decreased with blood sampling time and was negatively correlated with eGFR (Table 2). After adjustment for eGFR and blood sampling time, the negative associations with age and smoking were confirmed only in participants without asthma (beta=-0.04, p=0.06 and beta for trend=-0.08, p=0.006, respectively).

Associations between CC-16 and asthma

Univariate association between CC-16 level and current asthma (as reported in Table 1) was no longer confirmed after adjustment for eGFR (beta = 0.88 (95%CI 0.74, 1.05)) and when eGFR adjusted CC-16 z-score was considered as predictive variable (OR adj =0.92 (95%CI: 0.82, 1.02)). No association was observed between CC-16 z-score and asthma in participants not treated by inhaled corticosteroids (n=265, OR adj =0.94 (95%CI 0.82, 1.09)), nor by oral corticosteroids or anti-leukotriene (alone or combined, data not shown).

In participants with asthma, no association was observed between CC-16 z-score and asthma control (p>0.10) or asthma score (p=0.08 for score≥2) and CC-16 levels did not vary with ICS use (GM=11.7 (95%CI 8.03, 15.4) µg/L and GM=11.9 (95%CI 8.35, 15.5) µg/L, p =0.16, in participants no-ICS and ICS users respectively).

Associations between CC-16 and lung function

In participants with asthma, no association was observed between CC-16 and lung function (table E1, online), and this result was confirmed also in non-ICS users (n=253, p=0.10).

In participants without asthma (Table 3), reduced FEV 1 % predicted (borderline significant) and FEV 1 /FVC were associated with reduced CC-16 z-score, considered both on a continuous scale and categorized into quartiles (Figure E1, online). Association with lung function was higher in current smokers than in never smokers (Figure 2E, online), with statistically significant interactions for FEV 1 % predicted and FEV 1 /FVC.

Associations between CC-16 and AR

CC-16 GM level decreased from 12.2 (95%CI: 11.7, 12.7) µg/L in participants who performed the methacholine challenge test (N=862) to 11.1 (95%CI 9.75, 12.5) µg/L in 113 participants who did not performed the test because of poor lung function.

In participants with asthma no association was observed between CC-16 z-score and AR (table E1 online), and this result was confirmed also in non-ICS users.

In participants without asthma lower CC-16 level was associated with lower level in the doseresponse slope (Figure 2a) borderline higher risk of any AHR (PD20≤4mg) and high AHR (PD 20 ≤1mg) (Table 3). These results were confirmed also when FEV 1 % predicted was included in the multivariate model (beta for the slope: 0.20 (95%CI: 0.05, 0.35)) and when only participants without any asthma symptom (asthma score<2) were considered (n=465, beta for the slope =0.23 (95%CI 0.07, 0.39)). Similar results were observed when CC-16 zscore was grouped in quartiles (Figure E1, online). Participants in the highest CC-16 z-score quartile were associated with lower slope (beta adj =0.50 (95%CI: 0.08, 0.92)) and a lower risk of high AHR (OR adj = 0.30 (95%CI: 0.12, 0.75)) than participants in the lowest CC-16 z-score quartile.

Interaction between smoking and CC-16 z-score was statistically significant for the slope (p=0.02) (Figure 2E online). After adjustment for covariates, lower CC-16 z-score was associated with lower slope both in current as well as in ex-smokers (beta adj =0.50 (95%CI: 0.21, 0.79), beta adj =0.24 (95%CI: 0.01, 0.46), respectively).

DISCUSSION

Using data from a large epidemiological study on asthma, we observed that reduced lung function and higher AR were associated with decreased CC-16 level only in participants without asthma. No significant association was observed between CC-16 level and current asthma.

The major assets of this study are the sample size, the standardization of the sample collection and the availability of clinical parameters. Definition of asthma case is precise 29 leading to a very limited risk of false positives. Serum CC-16 level were similar to the reference values in healthy non-smokers (between 10 and 15 µg/L on average) 3 . Prevalence of airway hyperresponsiveness, measured by a methacholine challenge test was quite high in participants without asthma. A possible explanation is that part of the participants without asthma is first degree relatives of asthma cases. Nevertheless, this result is consistent with the relatively considerable number of asymptomatic participants with AHR reported in cross-sectional epidemiologic studies, ranging from 19.3 to 62.4%. 30 CC-16 level was standardized on eGFR so that all the associations were not affected by eGFR's sex and age-related variation. CC-16 decreased with smoking and pack-years smoked, varied with blood sampling time and was negatively associated with eGFR, consistent with previously published results 31 . However, with a cross-sectional approach, causal association between CC-16 and respiratory diseases could not be established.

We found no association between CC-16 level and asthma. Similarly, no association was observed with asthma in 859 participants to population-based multi-centre ECRHS 10 neither with asthma and bronchitis in 6531 men aged 67-77 year from a cross-sectional study 32 . To date, most evidence of an association between CC-16 and asthma comes from small (less than 100 participants) clinical studies, both in children 33,[START_REF] Laing | Association between plasma CC16 levels, the A38G polymorphism, and asthma[END_REF] and in adults 11,12 . A possible explanation for this apparent discrepancy is that participants from clinical studies are more likely to have severe asthma than participants from population-based studies. Further, clinical samples are known to be vulnerable to selection bias that may compromise the generalizability of findings.

We found that FEV 1 and FEV 1 /FVC were associated with serum CC-16 level only in participants without asthma. Previously, a positive association between serum CC-16 and FEV 1 and FVC was observed by in 402 fire-fighters from a cross-sectional study 14 and with FEV 1 % predicted and FEV 1 /FVC ratio in adults from the general population in the ECHRS study 10 , after adjustment for the same large range of covariates, including blood sampling time. The association between lung function parameters and CC-16 found in our study is consistent with previous published results on an association between lung function and markers of inflammation, such as C-Reactive Protein 35 .

We observed a negative association between CC-16 level and AR in adults without asthma that was confirmed in participants without any asthmatic symptom and was also independent from bronchial calibre. To the best of our knowledge this is the first study showing an association between CC-16 and asymptomatic AR. Previously, a negative association between AR and CC-16, apparently unnoticed by the authors, was observed in 4724 smokers with mild to moderate airflow limitation from the Lung Health Study 13 and in 402 fire-fighters, after adjustment for smoking and atopic status 14 . These findings may be related to mild inflammatory changes, that are present in airway mucosa of patients with AHR 36 that may lead to damages of airway epithelium and consequent decreased number of club cells in small airways 37 . Alternatively, genetically determined lower level of CC-16 and subsequent fragile epithelium may contribute to the development of asymptomatic AHR. In such case, methacholine or any other bronchoconstrictive substance may easily penetrate the lung epithelium and induce bronchial hyperreactivity. This hypothesis is supported by the findings of Taniguchi and colleagues 38 , that showed that the A38G polymorphism in the CC-16 gene was associated with both plasma CC-16 level and AR in 154 asymptomatic, young, healthy adults. In the EGEA study, we previously found that methacholine AHR was significantly related to asthma incidence in asymptomatic participants 9 . It could be postulated that asymptomatic AR might precede development of respiratory symptoms and asthma 36,37 through mechanisms related with airway inflammation , damage to epithelium integrity and structural changes in the airway mucosa 39 that are associated with decreased CC-16 levels in the airway. Future longitudinal studies, especially on paediatric population, will be useful to evaluate the role of CC-16 as a prognostic marker for asthma.

The lack of association between CC-16 level and lung function or AR in asthmatic group could be partly explained by the complex mechanisms, including inflammation and oxidative stress that characterizes this disease, leading to less variability of CC-16 level in participants with asthma and less sensitivity to variation in AR. Further, the differential association between CC-16 and lung function and AR between participants without and with asthma support the evidence of heterogeneity of AR in asymptomatic and symptomatic participants 30 .

In conclusion, we observed that decreased CC-16 serum level were associated with asymptomatic airway responsiveness and impaired lung function, only in participants without asthma. Those results suggested that serum level of CC-16 may help detecting early damages of the lung epithelium integrity. Further longitudinal studies will help in understanding whether this information will be used for early detecting asthmatic trait. * FEV 1 % predicted and FVC% predicted ranged from 0 to 1 and were obtained using Stanojevic reference equations 40 . † Beta estimates (95%CI) for a unitary increase in the GRF adjusted CC-16 z-score, estimated with GEE linear regression models ‡ Odds ratios (95%CI) for a unitary increase in the GRF adjusted CC-16 z-score, estimated with GEE logistic regressions models All estimates are reported as unadjusted and adjusted for age, smoking status and pack-years (never smokers, ex-smokers that smoked ≤15 pack-years, ex-smokers that smoked ≥15 pack-years, current smokers with <15 pack-years and current smokers with ≥15 pack-years), blood sampling time, sex (for FEV 1 /FVC, any and high AHR and AR) and height (for FEV 1 /FVC)

16-Aug-16 1

Serum CC-16 level is associated with asymptomatic Airway responsiveness in adults from the EGEA study SUPPLEMENTAL DATA

METHODS

Population and protocol

The EGEA combines a case-control and family study of adult and childhood asthma (http://egeaintranet.vjf.inserm.fr/). The first EGEA survey (EGEA1) was conducted from 1991 to 1995, and the protocol and descriptive characteristics have been described elsewhere

RESULTS

Characteristics of the sample

Participants with asthma were younger and smoked less frequently than participants without asthma. Lung function parameters were lower and AHR (PD 20 ≤4mg) was more frequent in participants with asthma than in those without asthma.

Among participants who performed the methacholine test (N=862), 800 (514 with asthma and 286 without) had information about their reduction in FEV1 after a provocative dose of methacholine of 4mg (PD20<4mg). 339 out of 800 (42.4%) had PD20<4mg and among those 138 (138/514=26.8)

were without asthma and 201 (201/286=70.3) were with asthma. For what concern PD20<1mg, data were available for 815 participants (523 without asthma and 292 with asthma). Among those, 67 (12.8) with asthma and 162 (55.4) controls reported PD20<1mg.

Participants who performed the methacholine test were more often females (54% vs. 45%, p=0.007), were younger (mean age=41 vs. 50 years, p<0.001), smoked more frequently (24.3 % vs.

19.6%, p=0.06), had higher level of eGFR (mean value: 103 vs. 95 mL/min*1.73m 2 , p<0.001, age and sex adjusted p<0.001) and higher FEV 1 (mean value: 3.5 vs. 2.9 L, p<0.001, age and sex adjusted p<0.001) than participants who did not (N=416).

  Inserm), Evry (Centre National de Genotypage) and Annemasse (Etablissement francais du sang) where biological samples are stored. They are indebted to all the individuals who participated, without whom the study would not have been possible.EGEA cooperative groupCoordination: V Siroux (epidemiology, PI since 2013); F Demenais (genetics); I Pin (clinical aspects); R Nadif (biology); F Kauffmann (PI 1992-2012). Respiratory epidemiology: Inserm U 700, Paris: M Korobaeff (Egea1), F Neukirch (Egea1); Inserm U 707, Paris: I Annesi-Maesano (Egea1-2); Inserm CESP/U 1018, Villejuif: F Kauffmann, N Le Moual, R Nadif, MP Oryszczyn (Egea1-2), R Varraso; Inserm U 823, Grenoble: V Siroux. Genetics: Inserm U 393, Paris: J Feingold; Inserm U 946, Paris: E Bouzigon, F Demenais, MH Dizier; CNG, Evry: I Gut (now CNAG, Barcelona, Spain), M Lathrop (now Univ McGill, Montreal, Canada). Clinical centers: Grenoble: I Pin, C Pison; Lyon: D Ecochard (Egea1), F Gormand, Y Pacheco; Marseille: D Charpin (Egea1), D Vervloet (Egea1-2); Montpellier: J Bousquet; Paris Cochin: A Lockhart (Egea1), R Matran (now in Lille); Paris Necker: E Paty (Egea1-2), P Scheinmann (Egea1-2); Paris-Trousseau: A Grimfeld (Egea1-2), J Just.

Figure 1 .Figure 2 .

 12 Figure 1. Flow chart of the EGEA population and the participants included in the study

Figure E1 . 8 Figure E2 .

 E18E2 Figure E1. Association between lung function, AR and GFR-standardized CC-16 z-score grouped in quartiles in participants without asthma

Table 1 :

 1 Function Tests, European Community for Steel and Coal. Official Statement of the European Respiratory Society. Eur Respir J Suppl 1993;16:5-40. Socio-demographic and clinical characteristics of the participants

	Characteristics	All	Without	With asthma p
			asthma
	N	1298	803	495
	Females, %	51.4	53.4	48.1	0.07
	Age, mean (SD), years	43.8 (16.4)	46.4 (15.8)	39.6 (16.5)	<0.001
	Smoking status
	and intensity, %
	Lifetime non-smoker	50.1	50.9	48.7	<0.001
	Ex-smoker, <15 pack-years	19.2	20.2	17.5
	Ex-smoker, ≥15 pack-years	8.31	8.53	7.94
	Smoker, <15 pack-years	16.5	12.9	22.4
	Smoker, ≥15 pack-years	5.9	7.4	3.46
	Asthma symptom score
	0	13.1	7.33	22.4	<0.001
	1	86.9	92.7	77.6
	≥2	32.9	63.3	9.78
	Lung function test, N	1278	791	487
	FEV 1 % predicted*, Mean (SD)	97.0 (17.1)	100.8 (15.2)	90.8 (18.2)	<0.001
	FVC% predicted*, Mean (SD)	99.7 (14.7)	100.8 (14.5)	98.0 (14.8)	0.001
	FEV 1 /FVC, Mean (SD)	77.9 (9.32)	79.8 (7.37)	74.9 (11.2)	<0.001
	FEV 1 %predicted*>80%, %	86.9	92.7	77.6	<0.001
	Methacholine	862	540	322
	challenge test, N
	Dose-response slope,
	N, Mean (SD)	827, 6.34(2.03)	527, 6.95 (1.81) 300, 5.27 (1.96) <0.001
	Any AHR (PD20≤4mg), N	800	514	286

AR: Airway responsiveness; AHR: Airway Hyper-responsiveness; CC-16: Club cell secretory protein; FEV 1 : forced expiratory volume in 1 second; FVC: forced vital capacity; GFR: Glomerular Filtration Rate; GM: Geometric Mean; PD 20 : 20% reduction in FEV 1 after methacholine challenge test (1 mg or 4 mg cumulative dose); SD: Standard Deviation; * Stanojevic reference equations.

Table 3 .

 3 Association between airway responsiveness and lung function and GFR-standardized CC-16 z-score in participants without asthma

		N	Unadjusted beta/	p	Adjusted beta	p
			OR (95%CI)		/OR (95%CI)
	Lung function measurements
	FEV 1 % pred. Values *, †	791	1.15 (0.05,2.24)	0.04	1.01 (-0.07,2.08)	0.07
	FVC % pred. Values *, †	791	0.57 (-0.45,1.60)	0.27	0.42 (-0.61,1.45)	0.43
	FEV 1 /FVC †	762 0.62 (0.12,1.12)	0.02	0.50 (0.06,0.95)	0.02
	Methacholine challenge test, N,
	Dose-response slope, Mean (SD)	527 0.25 (0.09,0.41)	0.001	0.21 (0.06,0.37)	0.005
	No AHR (reference)	376 1		1	-
	Any AR (PD20≤4mg) ‡ , %	138	0.83 (0.68,1.00)	0.048 0.82 (0.66,1.02)	0.07
	High AR (PD20≤1mg) ‡ , %	67	0.70 (0.54,0.90)	0.006	0.65 (0.48,0.87)	0.004
	AR: Airway responsiveness; AHR: Airway Hyper-responsiveness; CC-16: Club cell secretory protein; CI: Confidence Interval; FEV 1 : forced
	expiratory volume in 1 second; FVC: Forced Vital Capacity; GFR: Glomerular Filtration Rate; OR: Odds Ratios; PD 20 : 20% reduction in
	FEV 1 after methacholine challenge test (1 mg or 4 mg cumulative dose); SD: Standard Deviation;

  Glomerular filtration rate (eGFR) was estimated using the CKD-EPI (Chronic Kidney Disease Epidemiology Collaboration) formula4 (see below) and expressed as mL/min*1.73m 2 : eGFR = 141 X min(Scr/κ,1) α X max(Scr/κ,1) -1.209 X 0.993 Age X 1.018 [if female] X 1.159 [if black],

	16-Aug-16 2 16-Aug-16 3
	quality management approach was followed for the implementation of the EGEA2 data collection,
	and an International Organization for Standardization (ISO) 9001:2000 certification was obtained
	(http://www.afaq.org/certification5262711141114).
	Phenotype definitions where Scr is the serum creatinine level (mg/dL), κ= 0.7, if female or 0.9 if males, α=-0.329, if
	Inhaled corticosteroid use was defined by a positive answer to the question, ''Have you used females, or -0.411, if males
	inhaled corticosteroids to help your breathing at any time in the past year?'' with an exhaustive list
	of medications available in France at the time of the study.
	Participants were defined with (1) controlled asthma if all the following features were present: no
	1 . Briefly, 388 more than 1 instance per week of trouble breathing (defined by the answer to ''How often have you
	asthmatic cases, recruited in 5 chest clinics, 1244 first-degree relatives of cases, and 415
	population-based controls were recruited (total, n= 2047). A 12-year follow-up of this population
	was conducted from 2003 to 2007 (EGEA2) in 5 centres in France (Grenoble, Paris, Lyon,
	Marseille, Montpellier). Participants were contacted by postal questionnaire, and non-responders
	were further contacted by telephone (step 1). Among the alive cohort (n = 2002), 92% (n = 1845)
	completed the short self-questionnaire. Responders were invited to participate in the second phase
	of the study (step 2). The complete examination was mostly performed in clinical centers (n = 1316;
	85.3%). However, to improve the follow-up rate, some participants were examined at home (n = 72;
	4.7%) or answered questionnaires by phone (n=78; 5.0%) or mail (n=77; 5.0%). Finally, 77% of the
	alive cohort (n= 1543) completed at least a detailed questionnaire. In addition, 58 new family Biological phenotypes
	members where included in the study at the second survey. Written consent was obtained from all Serum CC-16 level was determined with latex immunoassay using a rabbit anti-CC-16 antibody
	participants at both surveys. Ethical approval to carry out the study was obtained for both surveys (Dakopatts, Glostrup, Denmark) and purified CC-16 protein as standards (home made at the
	from the relevant committees (Cochin Royal Hospital, Paris, for EGEA1, and Necker-Enfants Laboratory of Toxicology Applied Pharmacology, Faculty of Medicine, Catholic University of
	Malades Hospital, Paris, for EGEA2). Examination procedures included a detailed questionnaire, Louvain, Belgium). When pooled normal sera is fractionated by fast protein liquid chromatography
	with questions on asthma and respiratory symptoms, treatment, allergic rhinitis, active smoking, on Sephacryl S-200 (Pharmacia Biotechnology, Uppsala, Sweden), CC-16 elutes as a single
	and exposure to environmental tobacco smoke. Participants had blood samples that allowed to component with an apparent M r of ~16 kD, and is indistinguishable from the native protein. To
	measure total IgE in a centralized laboratory (n= 1421; 88.8%). Spirometry was performed by using avoid possible interference by complement, rheumatoid factor, or chylomicrons, the plasma
	a standardized protocol with similar equipment across centers according to the American Thoracic samples were pretreated by heating at 56º C for 30 min and by the addition of polyethylene glycol
	Society / European Respiratory Society guidelines 2 to measure FEV 1 (n= 1414; 88.3%). Skin prick (16%, 1:1 [vol/vol] and trichloroacetic acid (10%, 1:40 [vol/vol]). After overnight precipitation at
	tests to 12 aeroallergens (cat, Dermatophagoides pteronyssinus, Blattela germanica, olive, birch, 4ºC, the samples were centrifuged (2,000x g for 10 min) and CC-16 was determined in the
	Parieteria judaica, timothy grass, ragweed pollen, Aspergillus, Cladospo-rium herbarum, supernatants. All samples were analyzed in duplicate at two different dilutions. The assay has a
	Alternaria tenuis, Cypress) were performed in 1326 participants (82.8%). Sensitization was defined
	by the presence of at least 1 positive skin test (mean wheal diameter ≥3 mm). Strong efforts were
	made to standardize all the examination procedures across centers and to minimize missing data. A

had trouble with your breathing because of your asthma in the past 3 months?'') and no asthma attack in the last 3 months (defined by the answer to ''How many asthma attacks have you had in the past 3 months?''), no nocturnal symptoms (woken up because of asthma or by an attack of shortness of breath) in the last 3 months, use of short-acting β2-agonist inhalers ≤2/wk in the last 3 months, no use of oral corticosteroids in the past year, FEV 1 ≥80%predicted; (2) partly controlled asthma if 1 or 2 of these features were absent; and (3) uncontrolled asthma if ≥3 of these features were absent or if respiratory problems had caused hospital or emergency admissions in the past year or use of oral corticosteroids in the past year or ≥12 asthma attacks in the past 3 months. detection limit of 0.5 ng/ml and an average analytical recovery of 95%, with the intra-and interassay coefficients of variation ranging from 5 to 10% 3 . Serum creatinine was determined with standard traceable IDMS (isotope dilution mass spectrometry traceable) assay.

Table E1 .

 E1 Association between asthma, airway responsiveness and lung function and GFR-

	standardized CC-16 z-score in participants with asthma
		N	Unadjusted beta/	p	Adjusted beta	p
			OR (95%CI)		/OR (95%CI)
	Lung function measurements
	FEV 1 % predicted* ,a	791 1.53 (-0.21,3.27) 0.09 1.53 (-0.21,3.27) 0.19
	FVC % predicted* ,a	791 0.44 (-0.82,1.71) 0.49 0.44 (-0.82,1.71) 0.67
	FEV1/FVC

a 762 0.86 (-0.28,2.00) 0.14 0.86 (-0.28,2.00) 0.15

Methacholine challenge test, N

Acknowledgments

The authors thank all those who participated in the study and in the various aspects of the examinations and all those who supervised the study in all centers. The authors are grateful to the three CIC-Inserm units of Necker, Grenoble and Marseille which supported the study and where subjects were examined. They are also grateful to the three biobanks in Lille (CIC-

Funding: Research funded by the European Respiratory Society (Fellowship LTRF fellowship n°123-2011) and the NETA project (French Agency of Health Safety, Environment and Work (AFSSET, EST-09-15).

Data and quality management: Inserm ex-U155 (Egea1): J Hochez; Inserm CESP/U 1018, Villejuif: N Le Moual; Inserm ex-U780: C Ravault (Egea1-2); Inserm ex-U794: N Chateigner (Egea1-2); Grenoble: J Quentin-Ferran (Egea1-2).


Olivier Devillers

Louis Noizet

Walking in a Planar Poisson-Delaunay Triangulation: Shortcuts in the Voronoi Path

Keywords:

Let X n be a planar Poisson point process of intensity n. We give a new proof that the expected length of the Voronoi path between (0, 0) and (1, 0) in the Delaunay triangulation associated with X n is 4 π 1.27 when n goes to infinity; and we also prove that the variance of this length is O(1/ √ n). We investigate the length of possible shortcuts in this path, and defined a shortened Voronoi path whose expected length can be expressed as an integral that is numerically evaluated to 1.16. The shortened Voronoi path has the property to be locally defined; and is shorter than the previously known locally defined path in Delaunay triangulation such as the upper path whose expected length is 35/3π 2 1.18.

Introduction

The Delaunay triangulation is one of the most classical object of computational geometry and searching for paths in a point set using Delaunay edges is useful, e.g. for point location, nearest neighbor search [START_REF] Devillers | Walking in a Triangulation[END_REF], or routing in networks [START_REF] Bose | Online routing in triangulations[END_REF].

If the points are random, several walking strategies have been studied [START_REF] Bose | On the stabbing number of a random Delaunay triangulation[END_REF][START_REF] Broutin | Efficiently navigating a random Delaunay triangulation[END_REF][START_REF] Chenavier | Stretch Factor of Long Paths in a planar Poisson-Delaunay Triangulation[END_REF][START_REF] Devillers | The worst visibility walk in a random Delaunay triangulation is O( √ n)[END_REF][START_REF] Devroye | Expected time analysis for Delaunay point location[END_REF], in this paper we consider variations of a particular strategy called Voronoi path that consists in linking in order all the nearest neighbors of a point moving linearly from s to t where s and t are two points in the plane. We analyze these paths when s and t are two fixed points and when the point set is a Poisson point process of density n, possibly augmented by the two points s and t.

The Voronoi path is known to have an expected stretch factor 4 π 1.27 when n → ∞ [START_REF] Baccelli | Markov paths on the Poisson-Delaunay graph with applications to routing in mobile networks[END_REF], we provide an alternative proof of this result and prove that this length is quite stable by showing that the variance is small. Then we explore improvements on the Voronoi path by using some shortcuts. The length of one of this improved path can be expressed as an integral that we compute numerically getting an expected length of 1.16.

Any path in the Delaunay triangulation obviously yields an upper bound for the length of the shortest path. The best known upper bound being 35 3π 2 1.182 which is obtained as the length of a path called upper path [START_REF] Chenavier | Stretch Factor of Long Paths in a planar Poisson-Delaunay Triangulation[END_REF]. We say that a path is locally defined, if it can be decided if an edge belongs to the path between s and t by just knowing the neighborhood of the edge, s and t. Analyzing non locally defined paths, such as the shortest path is much more difficult than locally defined ones such as the upper path. Our improved Voronoi path is locally defined and gives a shorter alternative to the upper path.

Notations and Definitions

For a point set χ we define its Delaunay triangulation Del(χ) as the set of edges [p, q] with p, q ∈ χ such that there exist a disk D with D ∩ χ = {p, q}. One can remark that if there is such a disk, there is also such a disk so that p and q are on the boundary of the disk (shrink the first disk staying inside up to a position where the points are on the boundary).

The Voronoi diagram associated with χ is the tuple (R i ) i∈χ where ∀p ∈ χ; R p = {q ∈ R 2 /d(q, χ) = d(q, p)} (with d the Euclidean distance). R p is the Voronoi cell of seed p.

The Voronoi Path V P χ (s, t) between two points s and t is defined as the path formed by the seeds of the Voronoi cells intersecting the segment st (see Figure 1 for an example of Voronoi path). If s, t ∈ χ this path links s to t, otherwise it links the nearest neighbor of s to the nearest neighbor of t.

We denote M (p 1 , p 2 ) the intersection point between the bisector of p 1 and p 2 and the line (st). The ball centered at M (p 1 , p 2 ) passing through p 1 and p 2 is denoted B(p 1 , p 2 ) and its radius is denoted R(p 1 , p 2 ).

In the sequel our point set will be a Poisson point process X n of intensity n or the same set augmented by two points X = X n ∪ {s, t} where s = (0, 0) and t = (1, 0).

We denote p i:j the tuple of points (p i , p i+1 , . . . , p j ), and p i = ..j the same tuple of points verifying ∀k, l ∈ [i, j], p k = p l .

Expectation of Stretch Factor of the Voronoi Path

The first lemma states that the fact that s and t belong to the point set has a small influence on the length of the Voronoi path when n is big:

s t V P Figure 1: The Voronoi path Lemma 1. Let X := X n ∪ {s,
t} where X n is a Planar poisson point process of intensity n and s, t ∈ R 2 . Let (V P χ (s, t)) be the length of the Voronoi Path from s to t in Del(χ). Then

E [ (V P X (s, t))] = E [ (V P Xn (s, t))] + O Å s -t √ n ã Proof.
First, we remark that with very high probability 1 -e -n π 4 , the disk of diameter [st] contains some points of X n and thus any disk centered on [st] of the form B(•, •) cannot contain both s and t; we first assume this is the case, and no such ball does contain both s and t.

V P X (s, t) and V P Xn (s, t) only differ by few edges around s and t. Let sp i be the first edge of V P X (s, t) and p 1 , p 2 , . . . , p i be the i first vertices of V P Xn (s, t). First we remark that all p j are neighbors of s in Del(X). actually, by definition of the paths, there is a disk D i centered on a point in [st] with p i on its boundary, s inside and no points of X n nor t inside, this disk witnesses that sp i is a Delaunay edge in Del(X). Thus, to go from V P Xn (s, t) to V P X (s, t) we have to add one Delaunay edge incident to s: sp i and to remove few edges between neighbors of s in Del(X). The length variation can be bounded using triangular inequality

p 1 p 2 + p 2 p 3 + . . . + p i-1 p i + sp i ≤ sp 1 + 2 sp 2 + 2 sp 3 + . . . + 2 sp i-1 + 2 sp i which is O Ä s-t √ n ä [6, Prop. 2.2].
The same applies to the end of the path around t. In the rare case with an empty disk of diameter [st] almost the same reasoning applies except that the two parts of V P Xn (s, t) to be removed may overlap. O Ä

s-t √ n
ä is still an upper bound on the length of the removed part. Now the added part is just edge st of length one, but since it arises only with probability e

-n π 4 = o Ä s-t √ n ä the result still holds.
Theorem 2. X n is a Planar poisson point process of intensity n and s, t ∈ R 2 . Let (V P Xn (s, t)) be the length of the Voronoi Path from s to t in Del(X n ). Then

E î (V P Xn (s,t)) s-t ó = 4 π .
Proof. Without loss of generality, we may assume that s, t = (0, 0), (1, 0)

(V P Xn (s, t)) = 1 2 p 1 = ..2 ∈X 2 n 1 [p1p2∈V P Xn (s,t)] ||p 2 -p 1 ||, Inria (V P Xn (s, t)) = 1 2 p 1 = ..2 ∈X 2 n 1 [M (p1,p2)∈st] 1 [B(p1,p2)∩Xn=∅] ||p 2 -p 1 ||,
Using Slivnyak-Mecke formula, we transform this sum in an integral [START_REF] Schneider | Stochastic and Integral Geometry[END_REF]Theorem 3.3.5]:

E [ (V P Xn (s, t))] = n 2 2 (R 2 ) 2 1 [M (p1,p2)∈st] P [B(p1, p2) ∩ X n = ∅] ||p 2 -p 1 ||dp 1:2 .
Let Φ be the function

Φ : R × R + × [0, 2π) 2 -→ R 2 × R 2 (x, r, α 1 , α 2 ) -→ (p 1 , p 2 ),
where for i=1,2 we let p i = (x, 0) + r(cos α i , sin α i ).

As long as p 1 and p 2 do not have the same absciss, which occurs with probability 1, x is the absciss of M (p 1 , p 2 ). r is the distance between this point and p 1 . So Φ is a C 1 -diffeomorphism up to a null set. Its Jacobian is

det(J Φ ) = 1 cos α 1 -r sin α 1 0 0 sin α 1 r cos α 1 0 1 cos α 2 0 -r sin α 2 0 sin α 2 0 r cos α 2 = r 2 (cos α 2 -cos α 1 ). Since p 2 -p 1 = 2r sin α1-α2 2 and P [B(0, r) ∩ X n = ∅] = e -nπr 2 , we get E [ (V P Xn (s, t))] = n 2 2 ∞ -∞ ∞ 0 [0,2π) 2 1 [0<x<1] e -nπr 2 2r sin α 1 -α 2 2 |det(J Φ )| dα 1 dα 2 drdx = n 2 ∞ -∞ 1 [0<x<1] dx ∞ 0 e -nπr 2 r 3 dr [0,2π) 2 sin α 1 -α 2 2 |cos α 2 -cos α 1 | dα 1 dα 2 = n 2 1 2π 2 n 2 [0,2π) 2 sin α 1 -α 2 2 |cos α 2 -cos α 1 | dα 1 dα 2 = n 2 1 2π 2 n 2 × 8π = 4 π .
The above trigonometric integral is invariant by substituting

(α 2 , α 1 ) or (2π -α 1 , 2π -α 2 ) to (α 1 , α 2 ). Thus: [0,2π) 2 sin α 1 -α 2 2 |cos α 2 -cos α 1 | dα 1 dα 2 = 4 π 0 2π-α2 α2 sin α 1 -α 2 2 (cos α 2 -cos α 1 ) dα 1 dα 2 = 8π. RR n°8946
4 Variance of Stretch Factor of the Voronoi Path

Theorem 3. V ï (V P Xn ) s -t ò = O(n -1 2 ).
Proof. Once again, we assume that s = (0, 0) and t = (1, 0).

V [ (V P Xn )] = E (V P Xn ) 2 -E [ (V P Xn )] 2 , (V P Xn ) 2 = Ö 1 2 p 1 = ..2 ∈Xn 2 V P (p 1:2 ) è 2
,

where V P (p 1:2 ) = 1 [B(p1,p2)∩Xn=∅] 1 [M (p1,p2)∈st] p1 -p2 (V P Xn ) 2 = 1 4 Ö 2 p 1 = ..2 ∈Xn 2 V P (p 1:2 ) 2 è + 1 4 Ö 4 p 1 = ..3 ∈Xn 3 V P (p 1:2 ) V P (p 2:3 )) è + 1 4 Ö p 1 = ..4 ∈Xn 4 V P (p 1:2 ) V P (p 3:4 )) è . E    Ö p 1 = ..4 ∈Xn 4 V P Xn (p 1:2 ) V P Xn (p 3:4 )) è   = n 4 (R 2 ) 4 E V P X (p 1:2 ) V P X (p 3:4 ) dp 1:4 , where X = X n ∪ {p 1:4 } E    Ö p 1 = ..4 ∈Xn 4 V P Xn (p 1:2 ) V P Xn (p 3:4 )) è   = n 4 (R 2 ) 4 E 1 [B(p1,p2)∩X =∅] 1 [B(p3,p4)∩X =∅] 1 [M (p1,p2)∈st] 1 [M (p3,p4)∈st] p 1 -p 2 p 3 -p 4 dp 1:4 ≤ n 4 (R 2 ) 4 E 1 [B(p1,p2)∩Xn=∅] 1 [B(p3,p4)∩Xn=∅] 1 [M (p1,p2)∈st] 1 [M (p3,p4)∈st] p 1 -p 2 p 3 -p 4 dp 1:4 .

Inria

With the same substitution as previously, done twice, we get:

E    Ö p 1 = ..4 ∈Xn 4 V P Xn (p 1:2 ) V P Xn (p 3:4 )) è   ≤ n 4 [0,1] 2 R 2 + [0,2π) 4
e -n A(B((x,0),r)∪B((x ,0),r )) 2r sin

α 1 -α 2 2 2r sin α 3 -α 4 2 r 2 | cos α 2 -cos α 1 |r 2 | cos α 4 -cos α 3 |dα 1:4 drdr dxdx .
By rewriting the exponential as:

e -n A(B((x,0),r)∪B((x ,0),r )) = e -nπ(r 2 +r 2 ) + e -n A(B((x,0),r)∪B((x ,0),r )) -e nπ(r 2 +r 2 ) and applying Fubini's theorem, we get:

E    1 4 Ö p 1 = ..4 ∈Xn 4 V P Xn (p 1:2 ) V P Xn (p 3:4 )) è   ≤ E [ (V P Xn )] 2 + r n ,
where

r n = n 4 [0,1] 2 R 2 + [0,2π) 4
e -n A(B((x,0),r)∪B((x ,0),r )) -e -nπ(r 2 +r 2 ) 2r sin

α 1 -α 2 2 2r sin α 3 -α 4 2 r 2 | cos α 2 -cos α 1 |r 2 | cos α 4 -cos α 3 |dα 1:4 drdr dxdx .
Breaking the symmetry between r and r , we get

r n = 8n 4 [0,1] 2 R+
∞ r e -n A(B((x,0),r)∪B((x ,0),r )) -e -nπ(r 2 +r 2 ) r 3 r 3 drdr dxdx

× [0,2π) 4 sin α 1 -α 2 2 sin α 3 -α 4 2 (cos α 2 -cos α 1 )(cos α 4 -cos α 3 ) dα 1:4 .
Since we now have r ≤ r, we get e -n A(B((x,0),r)∪B((x ,0),r )) -e -nπ(r 2 +r 2 ) ≤ e -nπr 2 ,

r n ≤ 8(8π) 2 n 4 [0,1] 2 ∞ 0 r 0 r 6 e -nπr 2 1 [B(z,r)∩B(z ,r ) =∅] dr drdxdx ≤ 512π 2 n 4 1 0 x +2r x -2r ∞ 0 r 7 e -nπr 2 drdxdx ≤ 512π 2 n 4 ∞ 0 4r 8 e -nπr 2 dr ≤ 512π 2 n 4 • 4 105 32π 4 n 4 √ n = O(n -1 2 ). E    p 1 = ..2 ∈Xn 2 V P Xn (p 1:2 ) 2    = n 2 (R 2 ) 2 E î 1 [B(p1,p2)∩Xn=∅] 1 [M (p1,p2∈st)] p 1 -p 2 2 ó dp 1 dp 2 = n 2 1 0 ∞ 0 [0,2π) 2 e -nπr 2 4r 2 sin 2 α 1 -α 2 2 r 2 | cos α 1 -cos α 2 |dα 1 dα 2 drdx = 4n 2 1 0 dx ∞ 0 e -nπr 2 r 4 dr [0,2π) 2 sin 2 α 1 -α 2 2 | cos α 1 -cos α 2 |dα 1 dα 2 = 4n 2 3 4πn 5 2 16 3 = O(n -1 2 ). E    p 1 = ..3 ∈Xn 3 V P Xn (p 1:2 ) V P Xn (p 2:3 )    ≤ n 3 (R 2 ) 3 E î 1 [B(p1,p2)∩Xn=∅] 1 [B(p2,p3)∩Xn=∅] 1 [M (p1,p2∈st)] 1 [M (p2,p3∈st)] p 1 -p 2 p 2 -p 3 ó dp 1 dp 2 dp 3 ≤ 2n 3 (R 2 ) 3 e -nπR(p1,p2) 2 1 [ p2-p3 ≤ p1-p2 ] 1 [M (p1,p2∈st)] p 1 -p 2 p 2 -p 3 dp 1 dp 2 dp 3 ≤ 2n 3 (R 2 ) 3 e -nπR(p1,p2) 2 1 [p3∈B(M (p1,p2),3R(p1,p2))] 1 [M (p1,p2∈st)] 4R(p 1 , p 2 ) 2 dp 1 dp 2 dp 3 = 8n 3 (R 2 ) 2 e -nπR(p1,p2) 2 9πR(p 1 , p 2 ) 2 1 [M (p1,p2∈st)] R(p 1 , p 2 ) 2 dp 1 dp 2 = 72πn 3 1 0 dx ∞ 0 e -nπr 2 r 2 r 2 r 2 dr [0,2π) 2 | cos α 2 -cos α 1 |dα 2 dα 1 ≤ 72πn 3 • 15 16π 3 n 7 2 • 8 = O(n -1 2 ).
Combining these terms in the definition of V [ (V P Xn )] terminates the proof.

Improvement upon the Voronoi Path

We call shortcut of V P χ a triangle (p 1 , p 2 , p 3 ) such that (p 1 , p 2 ) and (p 2 , p 3 ) are in the Voronoi Path, and (p 1 , p 2 , p 3 ) is in Del(χ) (see Figure 2). Notice that it may exist other shortcuts replacing more than two edges in the Voronoi path, but the probability of existence decrease with the length of the replaced chain. In this paper we limit our interest to the above defined simple shortcuts.

Let SC (p 1 , p 2 , p 3 ) be defined as the length saved by taking the shortcut (p 1 , p 2 , p 3 ), i.e

SC (p 1 , p 2 , p 3 ) = p 1 -p 2 + p 2 -p 3 -p 1 -p 3 .
As shown on Figure 2 some shortcuts are incompatible, but the set of shortcuts can be divided in two sets, the shortcuts above the Voronoi path and the shortcuts below the Voronoi path that are compatible. By symmetry, the expected length of the above shortcuts is equal to the one of below shortcuts and is equal to half the total length of all shortcuts. Let gain Xn denote this expected saving in the Voronoi path for a Poisson point process X n .

E [gain

Xn ] = E    1 2 p 1 = ..3 ∈Xn 3 1 [p1:3∈Del(Xn)] 1 [p1:2∈V P Xn ] 1 [p2:3∈V P Xn ] 1 [xp 1 <xp 2 <xp 3 ] SC (p 1:3 )    .
We use the Slivnyak-Mecke formula:

E [gain Xn ] = n 3 2 (R 2 ) 3 E î 1 [p1:3∈Del(Xn)] 1 [p1:2∈V P Xn ∪{p 3 } ] 1 [p2:3∈V P Xn ∪{p 1 } ] 1 [xp 1 <xp 2 <xp 3 ] SC (p 1:3 ) ó dp 1:3 = n 3 2 (R 2 ) 3 P [(B(p 1:3 ) ∪ B(p 1 , p 2 ) ∪ B(p 2 , p 3 )) ∩ X n = ∅] 1 [M (p1,p2)∈st] 1 [M (p2,p3)∈st] 1 [p1 ∈B(p2,p3)] 1 [p3 ∈B(p1,p2)] 1 [xp 1 <xp 2 <xp 3 ] SC (p 1:3 )dp 1:3 = n 3 (R 2 ) 3 e -n A(B(p1:3)∪B(p1,p2)∪B(p2,p3)) 1 [M (p1,p2)∈st] 1 [M (p2,p3)∈st] 1 [yp 2 <0] 1 [p1 ∈B(p2,p3)] 1 [p3 ∈B(p1,p2)] 1 [xp 1 <xp 2 <xp 3 ] SC (p 1:3 )dp 1:3 .
We have limited our attention to half of the shortcuts using the assumption that y p2 < 0.

The assumption x p1 < x p2 < x p3 ensure that each triangle is counted only once. Under these two hypotheses, we consider only a subset of the possible shortcuts that verify y Ω > 0 (with Ω the center of B(p 1:3 )) and p 1 , p 2 , p 3 counterclockwise (ccw) then p 1 ∈ B(p 2 , p 3 ) and p 3 ∈ B(p 1 , p 2 ). Since we do not

B(p 1 , p 2 ) B(p 1:3 ) p 2 Ω p 3 p 1 B(p 2 , p 3 )
r consider all shortcuts any longer, in the sequel, we only have an upper bound on gain Xn . Actually our experiments show that this is a reasonably good upper bound.

E [gain Xn ] ≥ n 3 (R 2 ) 3
e -n A(B(p1:3)∪B(p1,p2)∪B(p2,p3))

1 [M (p1,p2)∈st] 1 [M (p2,p3)∈st] 1 [yp 2 <0] 1 [yΩ>0] 1 [p1,p2,p3 ccw] 1 [xp 1 <xp 2 <xp 3 ] SC (p 1:3 )dp 1:3 = E 1 .
Let r be the radius of B(p 1:3 ), S( yΩ r , α 1 , α 2 , α 3 ) = A(B(p1:3)∪B(p1,p2)∪B(p2,p3))

r 2
be the area of the union of the three balls normalized by r 2 , and h = yΩ r the normalized distance from Ω to line (st). Since Ω is assumed above line (st) and p 2 below it, h ∈ [0, 1]. We define a region for the angles α i :

1 h α 2 α 1 α 3 I h α 1 α 3 α 2 Ω I h =      (α 1 , α 2 , α 3 ) ∈ R 3      2π -α 2 < α 1 < α 2 π + arcsin h < α 2 < 2π -arcsin h α 2 -2π < α 3 < 2π -α 2     
. E 1 becomes:

E 1 = n 3 (R 2 ) 3 e -nr 2 S( yΩ r ,α1,α2,α3) 1 [M (p1,p2)∈st] 1 [M (p2,p3)∈st] 1 [0<yΩ<r] 1 [α1:3∈I h ] SC (p 1:3 )dp 1:3
Basic trigonometry gives:

x M (p1,p2) = x Ω - yΩ tan( α1+α2 2
)

and x M (p2,p3) = x Ω - yΩ tan( α2+α3 2 )
.

Defining

J yΩ,α1:3 =    ï yΩ tan( α1+α2 2 )
, 1 +

yΩ tan( α2+α3 2 ) ò , if yΩ tan( α1+α2 2 ) ≤ 1 + yΩ tan( α2+α3 2 )
∅ , otherwise

,

we have 1 [M (p1,p2)∈st] 1 [M (p2,p3)∈st] = 1 [xΩ∈Jy Ω ,α 1:3 ]
. We are now ready to substitute the variables using Blaschke-Petkantschin formula [START_REF] Schneider | Stochastic and Integral Geometry[END_REF]Theorem 7.3.1 ]. We get:

E 1 = n 3 [0,2π) 3 ∞ 0 R R e -nr 2 S( yΩ r ,α1,α2,α3) 1 [xΩ∈Jy Ω ,α 1:3 ] 1 [0<yΩ<r] 1 [α1:3∈I h ] 2r sin α1-α2 2 + sin α2-α3 2 -sin α1-α3 2 2r 3 A (α 1:3 ) dy Ω dx Ω drdα 1:3 , = 2n 3 [0,2π) 3 ∞ 0 r 0 Ç Jy,α 1:3 dx å e -nr 2 S( y r ,α1,α2,α3) 1 [α1:3∈I h ] g(α 1:3 )r 4 dydrdα 1:3 , where g(α 1:3 ) = 2 A (α 1:3 ) sin α1-α2 2 + sin α2-α3 2 -sin α1-α3 2 , E 1 ≥ 2n 3 [0,2π) 3 1 0 ∞ 0 (1 -rh g (α 1:3 )) e -nr 2 S(h,α1,α2,α3) r 5 1 [α1:3∈I h ] g(α 1:3 )drdhdα 1:3 , Inria θ 12 φ 12 r 12 r = 1 p 1 p 3 p 2 B(p 1:3 ) B(p 1 , p 2 )
By a very similar reasoning we get

θ 23 = 2 arccos( Ã 1 - Ç sin α2-α3 2 r 23 å 2 sign(2h + sin α 2 + sin α 3 )).

Value of the gain

Using the above expression for S(h, α 1:3 ), the integral of Equation (1) has been numerically approximated using Maple giving E [gain Xn ] 0.108. So the expectation of the length of the new path is 1.165.

Maple file is available with this research report.

Alternative Paths

As a concluding remark, we mention several possibilities of paths from s to t that can be defined in a Delaunay triangulation:

-the shortest path, -compass routing (vertex following v minimizing the angle with vt), -upper path (edges vw of triangles uvw with u below line (st) and v and w above), -greedy angle (vertex following v minimizing the angle with the horizontal through v admist the vertices of edges intersecting line (st)), -closest neighbor (vertex following v minimizing the distance to t), -the Voronoi path, -the Voronoi path with all possible shortcuts taken greedily, and -the Voronoi path with half the shortcuts (i.e. ccw shortcuts). Some of these paths are locally defined. i.e., the fact that vw belongs to the path can be decided knowing only s, t and some neighborhood of vw. Some are incremental, i.e., the vertex following v can be decided knowing that v is on the path, s, t, and some neighborhood of v.

Using CGAL [START_REF] Cgal | Computational Geometry Algorithms Library[END_REF] we experiment on the length of these paths (with st = 1) for a random set of points. The length of number of edges that we obtained after 1000 experiments with point

Inria

3 Figure 2 :

 32 Figure 2: Shortcuts Voronoi path of Figure 1 (in red).

  Inria density 10 6 are in the following table:

	Path	Experimental Number	Path	Theoretical bound
		Length	of edges	properties
	Shortest path Compass routing	1.041 1.068	927 956	incremental	∈ [1 + 10 -11 , 1.182] [6] Θ( √ n) edges [7]
	Greedy angle	1.097	997	incremental
	V P greedy shortcuts	1.130	995	incremental
	V P ccw shortcuts Closest neighbor walk	1.164 1.167	1081 873	locally defined incremental	this paper numerical integration] n) edges [7] 1.165 [ Θ( √
	Upper path	1.177	1072	locally defined	35 3π 2	1.182 [6]
	V P	1.274	1273	incremental locally defined	4 π	1.273 [1]

RR n°8946

rue du Jardin Botanique CS20101 54603 Villers-lès-Nancy Cedex Publisher Inria Domaine de Voluceau -Rocquencourt BP 105 -78153 Le Chesnay Cedex inria.fr ISSN 0249-6399

where g (α

-

. The length of J y,α1:3 is 1 -rh g (α 1:3 ) when the interval is nonempty; otherwise, 1 -rh g (α 1:3 ) is negative and still bounds the interval's length from below.

Since

the contribution of the term rhg to E 1 is negligible and we get

Determination of S(h, α 1 , α 2 , α 3 )

We will look at the different terms of this sum. First we remark that when α 1:3 ∈ I h , the two last terms disappear since B(u 1 , u 2 ) ∩ B(u 2 , u 3 ) ⊂ B(u 1:3 ) (the apexes of B(u 1 , u 2 ) ∩ B(u 2 , u 3 ) are p 2 and its symetric with respect to the line y = -h).

The first term is just π the area of the unit circle.

The second and third terms are πr 2 12 and πr 2 23 with r 12 = R(u 1 , u 2 ) and r 23 = R(u 2 , u 3 ). The fourth term is

The fifth term is similar to the fourth one.

Above undefined quantites can be expressed in term of h and α 1:3 . Since the angle of

) 2 . The angles φ are easy to compute: φ 12 = α 2 -α 1 and φ 23 = α 3 -α 2 + 2π. RR n°8946

The angle θ 12 verifies cos