T Boukhobza
email: taha.boukhobza@cran.uhp-nancy.fr

F Hamelin

Observability analysis and sensor location study for structured linear systems in descriptor form with unknown inputs

Keywords: Descriptor systems, right-hand side observability, sensor location, graph theory

This paper deals with the state and input right-hand side observability and the sensor location problem to ensure this property for structured linear systems in descriptor form. The proposed method is based on a graph-theoretic approach. It is applicable to a large class of descriptor systems including regular and nonregular systems even if they are underdetermined, overdetermined or square. After defining a bipartite graph which represents such systems, we deduce some graphic necessary and sufficient conditions to achieve the right-hand side observability of the state and the input. Next, we propose a sensor placement procedure to recover this property. More precisely, we provide the minimal number of the required sensors and either their pertinent location or a necessary and sufficient condition to check if a given sensor configuration is adequate or not to ensure state and input right-hand side observability. Our approach uses algorithms with polynomial complexity orders and assumes only the knowledge of the system's structure and so is suited for large scale systems.

Introduction

Descriptor systems result from a convenient and natural modelling process [START_REF] Luenberger | ADynamic equation in descriptor form[END_REF][START_REF] Lewis | A tutorial on the geometric analysis of linear timeinvariant implicit systems[END_REF] and their applications can be found in various fields [START_REF] Müller | Descriptor systems: Pros and cons of system modelling by differential-algebraic equations[END_REF] such as robotics, electrical circuit networks, biologic and economic systems. Many of the structural properties of descriptor systems are revisited in [START_REF] Aplevich | Implicit Linear Systems[END_REF][START_REF] Dai | Singular Control Systems[END_REF][START_REF] Yip | Solvability, controllability and observability of continuous descriptor systems[END_REF]. Among these properties, the observability analysis of the state and the unknown input of descriptor systems is a great interest mainly in the design of control laws, fault detection and isolation schemes, fault tolerant control strategies, supervision proposals and so on. In most cases, the studies concerning the observability of descriptor systems deal with algebraic and geometric tools, which assume the exact knowledge of the state space matrices characterizing the model of the system. In many modelling problems, these matrices have a number of fixed zero entries, determined by the physical structure, while the remaining entries are not exactly known. To study the properties of these systems in spite of the poor knowledge we have of them, the idea is that we only keep the zero/nonzero entries in the matrices of their state-space representation. Many interesting works on these models, called structured models, are related to the graph-theoretic approach and aim to analyse properties such as controllability, observability . . . [START_REF] Dion | Generic properties and control of linear structured systems: A survey[END_REF]. It results from these works that the graph-theoretic approach provides quite simple and elegant analysis tool. Concerning descriptor systems, the studied structured form is

Σ Λ :            E λ ẋ = A λ x + B λ u y = C λ x + D λ u z de f = E λ x (1)
where x ∈ R n , u ∈ R m and y ∈ R p are respectively the state vector, the unknown input vector and the output vector. E λ ∈ R q×n , A λ ∈ R q×n , B λ ∈ R q×m , C λ ∈ R p×n and D λ ∈ R p×m represent matrices which elements are either fixed to zero or assumed to be nonzero parameters noted λ i . The vector of these parameters is Λ = (λ 1 , λ 2 , . . . , λ h ) T and it is assumed that Λ can take any value in R h . If all parameters λ i are numerically fixed, we obtain a so-called admissible realization of Σ Λ . We say that a property is true generically for Σ Λ if it is true for almost all its realizations or equivalently for almost all parameters λ i . Here, " for almost all its realizations " is to be understood as " for all parameter values (Λ ∈ R h ) except for those in some proper algebraic variety in the parameter space ". Among the studies analysing structured systems in descriptor form by means of a graph-theoretic approach, we can cite (Reinschke and Wiedemann, 1997;[START_REF] Wiedemann | Strukturelle Zugänge Zur Analyse und Synthese Linearer Regelungssysteme in Deskriptorform[END_REF][START_REF] Murota | Observer design for general linear time-invariant systems[END_REF][START_REF] Van Der Woude | The generic canonical form of a regular structured matrix pencil[END_REF][START_REF] Boukhobza | Observability of structured linear systems in descriptor form: a graph-theoretic approach[END_REF]. In [START_REF] Murota | Observer design for general linear time-invariant systems[END_REF], the authors use bipartite graphs to characterize structural controllability, observability of descriptor systems. However, they make some assumptions about the studied systems which must be generically regular i.e. matrices E λ and A λ have to be square (q = n) and generically det[sE λ -A λ ] 0. In (Reinschke and Wiedemann, 1997), the authors use digraphs to characterize structural controllability of descriptor systems which must be generically regular and satisfy another constraint: no main diagonal element of [sE λ -A λ ] equals zero. In [START_REF] Wiedemann | Strukturelle Zugänge Zur Analyse und Synthese Linearer Regelungssysteme in Deskriptorform[END_REF], structural properties of regular descriptor systems, including the observability, are studied by assuming also that there is no zero element in the main diagonal of [sE λ -A λ ]. More recently, (van der Woude, 2002) studies the generic number of invariant zeros of a structured regular descriptor system. Many descriptor systems do not satisfy all these assumptions like the underdetermined ones which are very helpful to represent systems with unknown inputs. This fact motivates the works presented in [START_REF] Boukhobza | Observability of structured linear systems in descriptor form: a graph-theoretic approach[END_REF], where necessary and sufficient conditions for the total or partial state observability of possibly non-regular descriptor systems are given using a digraph.

Unfortunately, to check these conditions, a recursive algorithm with an exponential complexity order must be used.

In this context, the aim of this paper is, first, to extend the works of [START_REF] Murota | Observer design for general linear time-invariant systems[END_REF] to establish the right-hand side observability for as well regular and non-regular, which can be possibly non-square, descriptor systems with unknown inputs, but contrary to [START_REF] Boukhobza | Observability of structured linear systems in descriptor form: a graph-theoretic approach[END_REF] by the use of algorithms having polynomial complexity orders. In a second part of the paper, using also a graph-theoretic approach, we study the number and the location of the additional sensors which are necessary and useful to recover the right-hand side input and state observability conditions. We answer to the following question: when a descriptor system is not input and state right-hand side observable, where must we place a minimal number of additional sensors to recover this observability property ? At this aim, we give the necessary and sufficient number of additional sensors we must add. All the proposed results are based on algorithms with polynomial complexity order. This may be an important criterion when we deal with large scale systems. Moreover, since we consider structured systems, our approach can be used during a conception stage. The paper is organized as follows: after section 2, which is devoted to the problem formulation, a bipartite graph representation of structured systems is given in section 3. Our main results are presented in Section 4. A brief conclusion ends the paper.

Problem statement

First of all, let us precise the assumptions made on the system and the type of observability we study. Obviously, u(t) and Ex(0 -) are assumed to be admissible and are such that system Σ Λ is solvable [START_REF] Yip | Solvability, controllability and observability of continuous descriptor systems[END_REF][START_REF] Geerts | Solvability conditions, consistency and weak consistency for linear differential-algebraic conditions and timeinvariant singular systems: The general case[END_REF]. Moreover, without loss of generality, we only assume that E λ does not have any row identically zero. Note that no assumption is done on the rank of [sE λ -A λ ] neither on its main diagonal. Therefore, contrary to the works which assume the regularity of the system, our approach is adapted to analyse the state and also the input observability [START_REF] Hou | Input observability and input reconstruction[END_REF] of systems with unknown inputs since we do not made any hypothesis on the regularity of Σ Λ . The kind of observability we study is related to the right-hand side observability (Hou and Müller, 1999a).

Hereafter, we recall the definition of the right-hand side observability. This definition extends the one given in (Hou and Müller, 1999a) for non-structured descriptor systems: Definition 1 Structured system Σ Λ is generically state and input right-hand side observable, if, for almost all the realizations of structured system Σ Λ , y(t) = 0 for t ≥ 0 implies x(0 + ) = 0 and u(0 + ) = 0.

It is proved in [START_REF] Hou | Observer design for descriptor systems[END_REF], that it ensures the existence of a generalized observer which allows to reconstruct the state and the input. In fact, Σ Λ is state and input right-hand side observable iff for t > 0, state x(t) and input u(t) can be uniquely expressed using the measurements and their derivatives [START_REF] Müller | On the observer design for descriptor systems[END_REF][START_REF] Hou | Observer design for descriptor systems[END_REF].

To summarize, considering a structured system in descriptor form Σ Λ , we propose first to answer the question whether or not state x(t) and input u(t) are generically right-hand side observable. The second and more original part of this paper is dedicated to the location of additional sensors to recover the right-hand side observability conditions. More precisely, starting from a system which is not generically state and input right-hand side observable, we propose a systematic procedure for additional sensor placement in order to recover the generic state and input right-hand side observability.

3 Graphic representation of a structured descriptor system

We associate to structured system Σ Λ a bipartite graph noted

B(Σ Λ ) = (V + , V -, W)
, where V + and V -are two disjoint vertex subsets and W is the edge set. The vertices are associated to the internal state, the input and the output components of Σ Λ and the edges represent links between these variables. More precisely,

V + = X ∪ U and V -= Y ∪ Z, with X = {x 1 , x 2 , . . . , x n }, U = {u 1 , u 2 , . . . , u m }, Z = {z 1 , z 2 , . . . , z q }, Y = {y 1 , y 2 , . . . , y p }. Edge set W = A-edges ∪ B-edges ∪ C-edges ∪ D-edges ∪ E-edges, where A-edges = (x j , z i ) | A λ (i, j) 0 , B-edges = (u j , z i ) | B λ (i, j) 0 , C-edges = (x j , y i ) | C λ (i, j) 0 , D-edges = {(u j , y i ) | D λ (i, j)
0 and E-edges = (x j , z i ) | E λ (i, j) 0 . Here, M λ (i, j) is the (i, j)th element of matrix M λ and (v 1 , v 2 ) denotes a directed edge from vertex v 1 ∈ V to vertex v 2 ∈ V. Note that to differentiate the state, the output and the dynamic variables from their corresponding vertices, the latter are written in bold font. Hereafter, we illustrate the proposed digraph representation with an example.

Example 2 Consider the structured system defined by:

E λ =                         λ 1 λ 2 0 0 0 0 λ 3 λ 4 0 0 0 0 λ 5 0 0 0 0 0 λ 6 λ 7                         , A λ =                         0 0 0 λ 8 λ 9 0 0 λ 10 0 λ 11 λ 12 0 0 0 0 0 0 0 0 λ 13                         , B λ =                         0 0 λ 14 0 λ 15 0 λ 16 0                         , C λ =                 λ 17 λ 18 0 0 0 0 0 λ 19 0 0 0 0 λ 20 0 0                 and D λ =                 0 0 0 0 0 λ 21                
. To such a model, we associate the digraph in Figure 1.

Let us now give some useful definitions and notations.

• Two edges are disjoint if they have no common vertex. In Example 2, edges (u 2 , y 3 ) and (x 3 , y 2 ) are disjoint but not (u 2 , y 3 ) and (x 3 , y 3 ) neither (x 3 , y 3 ) and (x 3 , y 2 ).

Let V 1 and V 2 denote two subsets of V. The cardinality of V 1 is noted card V 1 . • θ V 1 , V 2 is the maximal number of disjoint edges from V 1 to V 2 . • A matching in a bipartite graph B(Σ Λ ) = (V + , V -, W) is an edge set M ⊆ W such
that all the edges of M are disjoint. A matching is maximal if it has a maximal cardinality which is equal to θ V + , V -. We can associate to each maximal matching M, a non bipar-

tite digraph noted B M (Σ Λ ) = (V + , V -, W) where (v 1 , v 2 ) ∈ W ⇔ (v 1 , v 2 ) ∈ W or (v 2 , v 1 ) ∈ M. We denote by ∂ + M (resp. ∂ -M) the set of vertices in V + (resp. in V -) covered by the edges of M. We note S + 0 = V + \ ∂ + M and S - 0 = V -\ ∂ -M . In B M (Σ Λ )
, two vertices v i and v j are said to be strongly connected if it exists a path from v i to v j and a path from v j to v i . The relation "is strongly connected to" is an equivalence relation and we can define its equivalence classes. We call each equivalent class a strongly connected component. These strongly components are well-known in the graph theory [START_REF] Murota | Observer design for general linear time-invariant systems[END_REF]. They can be ordered using a partial order relation " " defined between two strongly connected components V i and V j as V i V j if there exists a path from an element of V j to an element of V i . We can decompose uniquely the obtained bipartite graph into partially ordered irreducible components denoted B i = (V + i , V - i , W i ), using Dulmage-Mendelsohn decomposition [START_REF] Dulmage | Coverings of bipartite graphs[END_REF][START_REF] Murota | Observer design for general linear time-invariant systems[END_REF] described below to make the paper self-contained:

Find a maximal matching r, let V i be the strongly connected component of the graph obtained from B M (Σ Λ ) after deleting the vertices of V 0 and V ∞ .

M in B(Σ Λ ), Let V + 0 = S + 0 ∪ {v ∈ V + , ∃ a path in B M (Σ Λ ) from S + 0 to v} V - 0 = {v ∈ V -, ∃ a path in B M (Σ Λ ) from S + 0 to v}. V 0 de f = V - 0 ∪ V + 0 . Let V + ∞ = {v ∈ V + , ∃ a path in B M (Σ Λ ) from v to S - 0 } V - ∞ = S - 0 ∪ {v ∈ V -, ∃ a path in B M (Σ Λ ) from v to S - 0 }. V ∞ de f = V - ∞ ∪ V + ∞ . For i = 1, . . . ,
For i = 1, . . . , r,

V + i = V i ∩ V i and V - i = V -∩ V i .
It is important to note that the obtained subsets V + i and V - i for i = 0, 1, 2, . . . , r, ∞ are the same whatever the choice of the maximal matching M [START_REF] Dulmage | Coverings of bipartite graphs[END_REF][START_REF] Murota | Observer design for general linear time-invariant systems[END_REF]. B 0 is called the minimal inconsistent part and B ∞ is the maximal inconsistent part.

In Example 2, since M covers all the vertices of V - and V + , we have that V 0 = V ∞ = ∅. Furthermore, we construct B M (Σ Λ ) by adding the following edges to B(Σ Λ ): {(y 1 , x 1 ), (z 1 , x 2 ), (z 4 , x 4 ), (z 2 , x 5 ), (z 3 , u 1 ), (y 2 , x 3 ), (y 3 , u 2 )}. we have the following decomposition:

V + 1 = {u 2 }, V - 1 = {y 3 }, V + 2 = {x 3 }, V - 2 = {y 2 }, V + 3 = {u 1 , x 1 , x 2 , x 4 , x 5 } and V - 3 = {y 1 , z 1 , z 2 , z 3 , z 4 }. • We say that component B i contains an E-edge if there ex- ists an E-edge connecting an element of V + i to another one of V - i . • We define B(Σ Λ ) = (V + , V -, W) a bipartite graph as the subgraph of B(Σ Λ ) = (V + , V -, W) from which we re- move the E-edges set i.e. V + = X ∪ U, V -= Y ∪ Z and W = A-edges ∪ B-edges ∪ C-edges ∪ D-edges. • We denote by Bi = ( V+ i , V- i , Wi
), i = 0, 1, . . . , r, ∞ the irreducible components of B(Σ Λ ) obtained after using Dulmage-Mendelsohn (D-M) decomposition. For system of Example 2, the results of D-M decomposition are:

V0 = V∞ = ∅, V+ 1 = {u 2 }, V- 1 = {y 3 }, V+ 2 = {x 3 }, V- 2 = {y 2 }, V+ 3 = {x 2 }, V- 3 = {y 1 }, V+ 4 = {x 1 }, V- 4 = {z 3 }, V+ 5 = {x 4 }, V- 5 = {z 1 }, V+ 6 = {u 1 , x 5 } and V- 6 = {z 2 , z 4 }. 4 Main results

Recalls on the state right-hand side observability of regular descriptor systems without unknown input

Let us denote

P(s) =       A λ -sE λ B λ C λ D λ     
 the pencil matrix of system Σ Λ . Generic rank of matrix P(s), denoted g_rank(P(s)) depends on s. Hence, g_rank(P(s)) = r, ∀ s ∈ C means that for almost all parameter values λ ∈ R h , rank(P(s)) = r, ∀ s ∈ C. The generic right-hand side observability of structured sys- (Hou and Müller, 1999a). To establish the state and input right-hand side observability for the most general descriptor systems, we first recall the graphical results established in [START_REF] Murota | Observer design for general linear time-invariant systems[END_REF] for regular descriptor systems without unknown input:

tem Σ Λ is equivalent to g_rank       A λ -sE λ B λ C λ D λ       = n + m, ∀s ∈ C,
Σ R Λ :                E λ ẋ = A λ x y = C λ x z de f = E λ x (2)
where

E λ ∈ R n×n and A λ ∈ R n×n , with det(A λ -sE λ ) 0.
Then, we extend the obtained results to handle the state and input right-hand side observability of any descriptor system. The following result is given in [START_REF] Murota | Observer design for general linear time-invariant systems[END_REF]:

Theorem 3 Structured system Σ R Λ is generically right-hand side observable or equivalently g_rank       A λ -sE λ C λ      
= n iff the two following conditions are satisfied : Cond.1) after the D-M decomposition of B(Σ Λ ), V0 = ∅, Cond.2) after the D-M decomposition of B(Σ Λ ), ∀i = 1, . . . , r, component B i does not contain any E-edge i.e. there is no E-edge linking V + i to V - i .

Input and state right-hand side observability of descriptor systems

Consider now the more general case of structured descriptor systems with unknown inputs Σ Λ defined by Equations ( 1). The generic right-hand side observability of structured system Σ Λ is equivalent to

g_rank       A λ -sE λ B λ C λ D λ       = n + m, ∀s ∈ C (3)
Therefore, a necessary graphical condition to the generic right-hand side observability is that in bipartite graph B(Σ Λ ) associated to Σ Λ , we have 3) is satisfied for s = 0. When V 0 = ∅, then, since by assumption g_rank(E λ ) = q, there exists necessarily a complete matching in B(Σ Λ ) including all the E-edges. In this case, this complete matching includes necessarily n + mq edges associated to the submatrix C λ D λ . After some row permutations, we can decompose matrix (C λ D λ ) in order to obtain

V 0 = ∅ or equivalently, θ(X ∪ U, Z∪Y) = n+m. Indeed, this condition is necessary to have in B(Σ Λ ), V0 = ∅ which is equivalent to g_rank       A λ B λ C λ D λ       = n + m i.e. that condition (
      C λ s D λ s C λ r D λ r     
 such that, in the associated bipartite graph, a complete matching covers exactly the edges associated to the square matrix

      A λ -sE λ B λ C λ s D λ s       . Let us denote Ēλ =       E λ 0 0 0       , Āλ =       A λ B λ C λ s D λ s       and Cλ = C λ r D λ r .
The obtained system is then square and regular, from (van der Woude, 2002) (Theorem 6), because there exists a complete matching in the bipartite graph associated to Āλs Ēλ . Note that the Ē-edges set is exactly equal to the original E-edges set. All the hypotheses or usability conditions of Theorem 3 are satisfied and we can then establish: Proposition 4 Descriptor structured system with unknown inputs Σ Λ is generically input and state right-hand side observable or equivalently g_rank

      A λ -sE λ B λ C λ D λ       = n + m iff the two following conditions are satisfied : Cond.1) after the D-M decomposition of B(Σ Λ ), V0 = ∅, or equivalently θ(X ∪ U, Z ∪ Y) = n + m in B(Σ Λ ) (Complete matching condition) Cond.2) after the D-M decomposition of B(Σ Λ ), ∀i = 1, . . . , r, component B i does not contain any E-edge i.e. there is no E-edge linking V + i to V - i (s-condition). In Example 2, condition Cond.1 of Proposition 4 is satisfied but not Cond.2. Indeed, even if between V + 1 and V - 1 as between V +
2 and V - 2 there is no E-edge, there are edges between V + 3 and V - 3 ((x 1 , z 1 ) or (x 2 , z 1 ) for instance).

Input and state right-hand side observability recovering by additional sensor implementation

When the conditions of Proposition 4 are not satisfied, our proposed procedure to recover the input and state right-hand side observability consists of two steps which correspond one by one to the conditions of Proposition 4. In the sequel, we define a new output vector Υ representing the additional sensors collecting the new measurements Υ(t) = H λ x x(t) + H λ u u(t). Hence, we denote the completed system by Σ c Λ :

Σ c Λ :              E λ ẋ(t) = A λ x(t) + B λ u(t) y(t) = C λ x(t) + D λ u(t) Υ(t) = H λ x x(t) + H λ u u(t)
The additional sensor components can be represented by vertex set Υ and edge subsets H x -edges and H u -edges from respectively X to Υ and U to Υ. We subdivide vertex subset Υ into two subsets denoted Υ 1 and Υ 2 . Each subset Υ i , i = 1, 2 corresponds to the sensors possibly added at the i th stage of our procedure.

Complete matching condition recovering

First of all, note that the complete matching condition is also equivalent to have in B(Σ Λ )

∀Ω ⊆ X ∪ U, θ Ω, Z ∪ Y = card Ω (4) We say that there is a dilation in the digraph of the system when the latter condition is not satisfied for some vertex subset Ω ⊆ X∪U. The aim of additional sensors at this first stage is to eliminate all the dilations. Then, the fact that condition Cond.1 of Proposition 4 is not satisfied i.e. θ X∪U, Z∪Y < card(X) + card(U) = n + m = card(V + ) in B(Σ Λ ) or in other words θ V + , V -< card(V + ) in B(Σ Λ ), implies that some additional sensor vertices are needed to complete V -. This number is at least equal to card( V+ 0 )card( V-0 ). In addition to the number of required sensors, another problem is to precise as finely as possible their location because obviously adding anywhere card( V+ 0 )-card( V-0 ) sensors does not allow to recover the complete matching condition. Using the D-M decomposition of the bipartite graph B(Σ Λ ), we can state:

Proposition 5 Consider structured linear system Σ Λ represented by bipartite graph B(Σ Λ ). To recover the complete matching condition, the minimal number of additional sensors, noted

Υ 1 in bipartite graph B(Σ c Λ ), is equal to γ = n + m -θ V + , V -= card V+ 0 -card V- 0 .
These additional sensors must measure γ states and unknown inputs in V+ 0 such that we obtain a maximal matching of size n

+ m in B(Σ c Λ ) i.e. θ X ∪ U, Z ∪ Y ∪ Υ 1 = n + m in B(Σ c Λ )
. Proof: Sufficiency: The Dulmage-Mendelsohn decomposition characterizes all the maximal matchings [START_REF] Murota | Observer design for general linear time-invariant systems[END_REF] and so if card V+ 0card V-0 edges are added between S+ 0 and the additional sensors, in order to form a matching of size card V+ 0card V-0 , then in B(Σ c Λ ), S+ 0 = ∅. Thus, there is no minimal inconsistent part B0 in B(Σ c Λ ). In this case, we obtain a maximal matching of size n + m between V + and Z ∪ Y ∪ Υ 1 in B(Σ c Λ ). Necessity: Let us consider a solution which provides a maximal matching M c of size n + m in B(Σ c Λ ). M c covers all the state and input vertices of V + and in particular the elements of V+ 0 . Thus, there are at least card( V+ 0 )card( V-0 ) edges from V+ 0 to the additional sensors. It can be useful to specify more precisely the additional sensor location to recover the complete matching condition. At this aim, let us define by BN (Σ Λ ) the non-directed graph corresponding to B(Σ Λ ). We call V0 -path every simple path of BN (Σ Λ ) which covers only vertices of V0 . For each V0 -path P, we define a vertex subset ϑ(P) such that ϑ(P) = {v ∈ V+ 0 , such that P covers v}. Finally, we say that a V0 -path P is maximal, if there does not exist a V0 -path P such that ϑ(P) ⊂ ϑ(P ). It is quite easy to see that to recover the complete matching condition, it is necessary to have in B(Σ C Λ ): for each maximal V0 -path P, θ ϑ(P), Υ 1 0 (5)

It is also possible to prove that, to recover the complete matching condition, it is necessary to have, in B(Σ c Λ ), for each maximal V0 -path P:

θ V+ 0 , Υ 1 ∪ V- 0 -θ V+ 0 \ ϑ(P), Υ 1 ∪ V- 0 > 0 (6)
Conditions ( 5) and ( 6) ensure that there is at least one sensor dedicated to take a measurement in each subset ϑ(P). In fact, these above conditions can be seen as complementary tools which allow to have effortlessly a better precision on the sensor location, since we obtain a kind of repartition of the required γ sensors.

Example 6 Consider the structured linear descriptor system represented by the bipartite graph depicted in Figure 2, where E-edges set is {(x 1 , z 1 ), (x 1 , z 4 ), (x 2 , z 2 ), (x 3 , z 3 ), (x 4 , z 4 ), (x 5 , z 4 )} and A-edges set is {(x 2 , z 1 ), (x 4 , z 1 ), (x 3 , z 2 ), (x 3 , z 3 ), (x 1 , z 4 ), (x 4 , z 4 ), (x 5 , z 4 )}. We have the following decomposition in B(Σ Λ ):

V+ 0 = {u 1 , x 3 , x 6 , u 2 , x 5 }, V-
In this paper, we propose a graph-theoretic tool to analyze, at first, the state and input generic right-hand side observability of structured linear descriptor systems. The studied system can be underdetermined, overdetermined or square and possibly non-regular. Necessary and sufficient conditions for state and unknown input generic observability are expressed in graphic terms. These conditions extend the ones of [START_REF] Murota | Observer design for general linear time-invariant systems[END_REF] which deal with regular descriptor systems without input. Next, when these observability conditions are not satisfied, we study the minimal number of the required additional sensors and either their location or necessary and sufficient conditions to be satisfied by any acceptable location in order to recover the generic state and input observability. It is shown also that our approach uses classical programming techniques and is free from numerical dif-ficulties since it uses well-known combinatorial techniques. Indeed, all the algorithms used in the paper have polynomial complexity orders. This makes our proposed method well-suited to large scale systems.

Figure 1 .

 1 Figure 1. Bipartite graph B(Σ Λ ) associated to the system of Example 2

Figure 2 .

 2 Figure 2. Bipartite graph B(Σ Λ ) representing system of Example 6

x 4 } and V-2 = {y 2 , z 1 }. Furthermore, in B(Σ Λ ), we have that V+ 0 = {x 3 , x 6 , u 2 , x 5 }, V-0 = {z 2 , z 4 }. As V+ 0 ∅, condition Cond.1 of Proposition 4 is not satisfied. The number of necessary sensors to recover this condition is card( V+ 0 )card( V-0 ) = 2 and they must take their measurements in {x 3 , x 6 , u 2 , x 5 }. However, if we add two sensors such that the first one measures x 3 and the second one x 6 , condition Cond.1 of Proposition 4 will not be satisfied. If we compute the maximal V0 -paths, we have : P 1 = x 6 → z 2 → x 3 and P 2 = u 2 → z 4 → x 5 . Therefore, the set of additional sensors Υ 1 must satisfy:

In the case of Example 6, condition (6) is satisfied if the above constraints are also satisfied.

s-condition recovering

Assume that condition Cond.1 of Proposition 4 is verified by the system. If it is not the case, then we must add first some sensors as it is described in the previous subsection before dealing with the s-condition recovering stage. We denote by Ỹ = Y ∪ Υ 1 the output vertex subset at the beginning of this second and last stage of our procedure. On the one hand, if condition Cond.1 is satisfied then V0 = ∅ and so V 0 = ∅. Thus, adding a sensor on V i , i = 1, . . . , r, makes that this strongly connected component disappears and all its vertices will belong to V ∞ . On the other hand, if condition Cond.1 is satisfied, when two vertices belong to the same strongly connected component, then, after adding a new sensor anywhere in the bipartite graph, they remain in the same strongly connected component. Therefore, the only way to satisfy Cond.2 of Proposition 4 is to add some sensors in order to make that all the vertices of strongly connected components B i , i = 1, . . . , r, which contain Eedges, will belong to B ∞ in completed system B(Σ c Λ ). This reasoning is summarized in the following proposition: Proposition 7 Consider structured linear system Σ Λ represented by bipartite graph B(Σ Λ ). Assume that the complete matching condition is satisfied i.e. V0 = ∅. To recover the s-condition, the additional sensors Υ 2 must measure at least one state in each subset V + i , i = 1, . . . , r, if B i contains an E-edge i.e. f or i=1, . . . , r, if B i contains an E-edge then we must have θ(V + i , Υ 2 ) 0 Note that it is possible that only one sensor, which takes its measurement in each strongly connected component containing an E-edge allows to satisfy the requirements of the s-condition. From a global point of view, the minimal number of additional sensors, required to recover both the complete matching condition and the s-condition, is then less or equal to γ+1, where γ = card( V+ 0 )-card( V-0 ) = n+m-θ X∪U, Y∪Z , and θ X ∪ U, Y ∪ Z is computed here in B(Σ Λ ) (i.e. without E-edges). In fact, we have two cases, if V + ∞ includes an element of V+ 0 then the minimal number of necessary sensors is γ. In this first case, we can use a sensor to recover both complete matching condition by putting it on an element of V + ∞ ∩ V+ 0 and satisfying the additional constraint of Proposition 7. Otherwise, if V + ∞ ∩ V+ 0 = ∅, the sensors used in the first stage allow just to transform set V0 into several sets V i , i 0 and i ∞. Then, these sensors cannot be useful to recover the s-condition which goal is to transform some V i sets (i 0 and i ∞) into a part of V ∞ . Therefore, we need a new sensor to recover s-condition. For Example 2, a sensor measuring any component in {u 1 , x 1 , x 2 , x 4 , x 5 } is necessary and sufficient to recover 5 the s-condition. in this case,

where Υ is the added sensor such that θ({u 1 , x 1 , x 2 , x 4 , x 5 }, Υ) 0. In the case of system described in Example 6, let us assume that to satisfy the complete matching condition, we have added in the first stage of our procedure two sensors measuring respectively x 6 and u 2 denoted respectively Υ 1 1 and Υ 2 1 . In this case, V 0 will be an empty set and V ∞ remains empty. The Dulmage-Mendelsohn decomposition gives the following components:

There are E-edges between V + 3 and V - 3 and between V + 4 and V - 4 and so B 3 and B 4 contains E-edges. Therefore, the additional sensor Υ 2 , must satisfy θ({u 1 , x 3 , x 6 } , Υ 2 ) = 1 and θ({u 2 , x 5 } , Υ 2 ) = 1. For this example, we can take Υ 2 = α 1 x 3 + α 2 x 5 .

Complexity analysis

From a computational point of view, our proposed approach needs few information about the system and is quite easy to check by means of well-known combinatorial techniques or simply by hand for small systems. More precisely, let us denote by M the number of edges in the digraph and N = n + m + q + p the number of vertices of our graphs. The first step of our procedure requires a Dulmage-Mendelsohn decomposition which can be implemented using an algorithm with a complexity order O(M 2 ) = O(n 4 ) [START_REF] Lovasz | Matching Theory[END_REF][START_REF] Chen | Constrained minimum vertex cover in bipartite graphs: Complexity and parameterized algorithms[END_REF], assuming without loss of generality that n ≥ p and n ≥ m. Next, our procedure requires the calculation of the strongly connected components which can be done using an algorithm which complexity order equals O(Nlog(N)) = O(nlog(n)) [START_REF] Fleischer | On Identifying Strongly Connected Components in Parallel[END_REF]. Checking the observability conditions necessitates only the knowledge of V 0 issued from the Dulmage-Mendelsohn decomposition and at most n × q tests to know if each E-edge belongs to a strongly connected component. This can be done with a O(n 3 ) complexity order algorithm. The additional sensor location characterization, does not necessitate more complex computation and the proposed method can be implemented using a global algorithm with a polynomial complexity equal to O(n 4 ).


Imelda Vickers
email: .imelda.vickers@cuh.ie

Darina O'flanagan

Mary Cafferkey

Hilary Humphreys

Multiplex PCR to determine Streptococcus pneumoniae serotypes causing otitis media in the Republic of Ireland with further characterisation of antimicrobial susceptibilities and genotypes

Keywords: Streptococcus pneumoniae, otitis media, serotype, PCR

Purpose: To determine the serotypes, genotypes and antimicrobial susceptibilities of Streptococcus pneumoniae causing otitis media (OM) in children in Dublin, Ireland. Methods: S. pneumoniae isolates (n = 28) from spontaneously discharging OM were studied. Serotyping was performed using a previously undescribed multiplex PCR scheme in combination with serological methods. Multilocus sequence typing (MLST) was performed using standard procedures. Antimicrobial susceptibility testing was performed using the Etest method.

Results: Fourteen different S. pneumoniae serotypes were identified. The five most common serotypes were 3, 19F, 19A, 14 and 6A which accounted for 68% of all infections. The 7-valent pneumococcal conjugate vaccine (PCV7), 10-valent pneumococcal conjugate vaccine (PHiD-CV) and 13-valent pneumococcal conjugate vaccine (PCV13) provided potential coverage of 43%, 46% and 86%, respectively.

Reduced susceptibility to penicillin was evident for 25% of isolates and associated with serotypes 14, 19A, 19F and 9V. A total of 21 different sequence types (STs) were identified. , 54% of which were representatives of Pneumococcal Molecular Epidemiology Network (PMEN) clones or their or their variants represented 54% (15/28) of all isolates. and 40% of these demonstrated reduced susceptibility to at least one antimicrobial agent.

Conclusions: Continued monitoring and characterisation of S. pneumoniae causing OM in

Ireland is warranted in order to guide future vaccine and treatment policies.

INTRODUCTION

Streptococcus pneumoniae is one of the leading causes of otitis media (OM) in children and is responsible for approximately 38% of all cases worldwide [START_REF] Abdelnour | Streptococcus pneumoniae serotype 3 among Costa Rican children with otitis media: clinical, epidemiological characteristics and antimicrobial resistance patterns[END_REF]. The majority (83%) of children less than 3 years of age experience at least one episode of OM [START_REF] Zhou | Trends in acute otitis media-related health care utilization by privately insured young children in the United States, 1997-2004[END_REF]. The cost of outpatient visits and treatments due to OM is considerably high, with an estimated expense of $2-5 billion per annum in the United States (US). The seven-valent pneumococcal conjugate vaccine (PCV7; Prevenar 7™) was or continues to be routinely administered to young children in a number of countries. An overall reduction of 6% and 8.9% in clinical OM incidence was demonstrated in Finnish and US infants following immunisation with PCV7 [START_REF] Fletcher | Brief review of the clinical effectiveness of PREVENAR against otitis media[END_REF]. Vaccine efficacy of 10-50% was evident against recurrent OM or for the prevention of tympanostomy tube placement. A decline of 42.7% and 41.9% in OM associated ambulatory visits and antibiotic prescriptions since vaccine introduction in the US, respectively, was demonstrated by Zhou et al. [START_REF] Zhou | Trends in acute otitis media-related health care utilization by privately insured young children in the United States, 1997-2004[END_REF].

PCV7 was introduced to the infant immunizsation programme in Ireland in late 2008. In order to determine the effects and impact of PCV7 in Ireland and indeed elsewhere, it is important to have an understanding and knowledge of S. pneumoniae serotypes causing infection. To our knowledge, no data associated with the epidemiology of S. pneumoniae causing OM in Ireland has been documented. Accurate S. pneumoniae identification is essential for diagnosis and for optimal treatment. Reliable typing systems are important to detect changes in the dominant strains over time. Our group previously developed a multiplex PCR assay to identify the seven serotypes contained in PCV7 [START_REF] O'halloran D M | Multiplex PCR for identification of seven Streptococcus pneumoniae serotypes targeted by a 7-valent conjugate vaccine[END_REF]. Here we describe a new and improved multiplex PCR scheme for use on lysed single colony extracts that confirms S. pneumoniae identification and detects an expanded number of serotypes. In addition, we describe the serotypes, genotypes and antimicrobial susceptibilities of S. pneumoniae causing OM in children in Ireland.

METHODS

DNA preparation.

A single colony from a fresh overnight culture grown on Columbia blood (5% sheep blood) agar plates was immersed in 20µl Microlysis™ Plus (Microzone Ltd., Sussex, UK) and processed according to manufacturer's guidelines. PCR was then performed on the lysed single colony extracts.

Multiplex PCR assembly. All available relevant sequences were aligned using ClustalX, and the appropriate primers (9V-F, 3-R, lytA-F and lytA-R) determined using Primer Express ® (Applied Biosystems, Warrington, UK) software. The lytA primers were designed to target conserved nucleotides in typical lytA alleles, characteristic of pneumococci not present in lytA alleles of mitis group streptococci [START_REF] Llull | Characteristic signatures of the lytA gene provide a basis for rapid and reliable diagnosis of Streptococcus pneumoniae infections[END_REF][START_REF] Obregon | Molecular peculiarities of the lytA gene isolated from clinical pneumococcal strains that are bile insoluble[END_REF][START_REF] Whatmore | Genetic relationships between clinical isolates of Streptococcus pneumoniae, Streptococcus oralis, and Streptococcus mitis: characterization of "Atypical" pneumococci and organisms allied to S. mitis harboring S. pneumoniae virulence factorencoding genes[END_REF]. The primer sequences and amplicon size of products were as previously described, with the exception of final primer concentrations in each reaction mixture; serotype 4 (0.75µM) and 14 (0.5µM) [START_REF] Brito | Serotyping Streptococcus pneumoniae by multiplex PCR[END_REF], 6 (0.6µM) [START_REF] Lawrence | Evaluation of semiautomated multiplex PCR assay for determination of Streptococcus pneumoniae serotypes and serogroups[END_REF], 18 (0.75µM), 19F (0.6µM), 23F (0.75µM), 1 (0.4µM), 5 (0.48µM), 7F (0.64µM) and cpsA (0.08µM) [START_REF] Pai | Sequential multiplex PCR approach for determining capsular serotypes of Streptococcus pneumoniae isolates[END_REF]. A serotype 9V forward primer was designed (5'-agttcaccttgtaacagttatcccata-3') (0.6µM) and used in combination with a previously described 9V reverse primer [START_REF] Pai | Sequential multiplex PCR approach for determining capsular serotypes of Streptococcus pneumoniae isolates[END_REF] (0.6µM), which yielded a product of 107bp.

A new serotype 3 reverse primer was also designed (5'-ttctgcccaccttagttgct-3') (0.4µM) and used in combination with a previously described forward primer [START_REF] Lawrence | Evaluation of semiautomated multiplex PCR assay for determination of Streptococcus pneumoniae serotypes and serogroups[END_REF] (0.4µM), which yielded a product of 470bp. The newly designed lytA forward primer (5'caactaaagaagagttcatga-3') (0.1µM) and lytA reverse primer (5'-ctcgtgcgttttaattccagc-3') (0.1µM) yielded a product of 125bp.

The first multiplex PCR reaction was designed to screen for serotypes contained in PCV7 with the use of cpsA primers as an internal positive control (IPC) [START_REF] Brito | Serotyping Streptococcus pneumoniae by multiplex PCR[END_REF][START_REF] Pai | Sequential multiplex PCR approach for determining capsular serotypes of Streptococcus pneumoniae isolates[END_REF] (Figure 1). The second multiplex PCR reaction was designed to confirm S. pneumoniae identification (lytA) and screen for serotypes 1, 3, 5 and 7F, with cpsA as an IPC.

Serotyping.

Isolates were serotyped using a combination of the multiplex PCR scheme described above, and/or for those isolates for which the serotype was not determined using the above method, the serotype was determined by serological methods [START_REF] Smart | Pneumococcal serotyping[END_REF] using antisera obtained from the Statens Serum Institute (Copenhagen, Denmark).

DNA isolation and multilocus sequence typing (MLST).

DNA was purified from fresh overnight cultures grown on Columbia blood agar using the Gentra ® Puregene ® Yeast/Bact. kit (Qiagen) according to manufacturer's guidelines.

MLST was performed using the standard procedure to amplify the seven house-keeping genes as previously described [START_REF] Enright | A multilocus sequence typing scheme for Streptococcus pneumoniae: identification of clones associated with serious invasive disease[END_REF]. Allele numbers and sequence types (STs) were assigned and eBURST analysis performed using www.mlst.net. Genetically related clones were either single locus variants (SLVs) or double locus variants (DLVs) of a specific clone. STs were compared with those of the PMEN clones available at http://www.sph.emory.edu/PMEN/ .

Antmicrobial susceptibility.

Isolates were tested for susceptibility to penicillin, cefotaxime, erythromycin, clindamycin and tigecycline using the Etest (AB bioMérieux, Solna, Sweden) method and results were interpreted following CLSI guidelines [START_REF]Performance standards for antimicrobial susceptibility testing; eighteenth informational supplement[END_REF]. Susceptibility to penicillin was defined as follows: susceptible, MIC ≤0.06µg/ml; intermediate, 0.12-1µg/ml; and highlevel resistant, ≥2µg/ml; corresponding to the new oral penicillin breakpoints. Isolates with reduced susceptibility to three antibiotics or more were considered multidrug resistant (MDR).

RESULTS

Multiplex PCR validation.

The multiplex PCR reactions correctly identified 73.5% (97/132) of all relevant serotypes targeted in the control collection of isolates (Table 1).

Three discrepancies did occur. However, repeat assessment confirmed the PCR result as accurate. The cpsA primers (internal positive control) were consistently positive for all isolates and the 90 pneumococcal serotypes tested (Figure 1). Exceptions included lack of reactivity with isolates harbouring the type 25A, 25F and 38 capsule, as previously observed [START_REF] Yu | A rapid pneumococcal serotyping system based on monoclonal antibodies and PCR[END_REF]. The lytA primers consistently and specifically amplified all pneumococci tested, with no cross reactivity with other organisms identified.

Patient characteristics.

Twenty eight cases of OM caused by S. pneumoniae were identified in children, age range 23 days to nine years. The majority (68%) of infections occurred in children <3 years. The male: female ratio was 1.5:1.

Serotype distribution.

A total of 14 different S. pneumoniae serotypes were identified. The most common serotypes were 3 (n=5) and 19F (n=5), followed by 19A (n=3), 14 (n=3) and 6A (n=3).

These five serotypes together accounted for 68% of infections overall. Three isolates caused infection after PCV7 introduction, one expressed vaccine serotype 14, another expressed vaccine serotype 9V, and the third isolate expressed serotype 6A. The potential coverage provided by PCV7 within all age groups included was 43%, PHiD-CV (Synflorix™, GlaxoSmithKline plc., UK) was 46% and PCV13 (Prevenar 13, Pfizer Inc., USA) was 86%. For those less than two years of age the potential coverage provided by both PCV7 and PHiD-CV was 53%, increasing to 87% for PCV13.

Sequence Types and Antimicrobial susceptibility.

Reduced susceptibility to penicillin was observed for 25% of the isolates, 86% of which occurred in children <3 years. One isolate (serotype 14) demonstrated high level penicillin resistance and caused an infection in a one year old boy. Four serotypes were associated with reduced susceptibility to penicillin, the most common being serotype 14 (n=3), followed by 19A (n=2), 19F (n=1) and 9V (n=1). Intermediate resistance to cefotaxime was observed for 3.6% of isolates. Resistance to erythromycin and clindamycin was observed for 18% and 11% of isolates, respectively. All isolates were sensitive to tigecycline. MDR was evident for 11% of isolates and all were members of serogroup 19.

Twenty-one STs were identified in total (Table 2). The predominant clone associated with serotype 3 was ST180 (n=3) and the genetically unrelated serotype 3 ST458 (n=1) was also identified. A high level of variation was observed among the clones associated with serotype 19F (ST420, ST425, ST73, ST87 and ST3508), all of which were considered unrelated to each other following eBURST analysis (data not shown). ST87 demonstrated intermediate resistance to penicillin and a very high level of resistance to the macrolide class of antibiotics, namely, erythromycin and also clindamycin (Table 3). All other clones expressing the serotype 19F capsule were sensitive to all antimicrobials tested. The two ST276 (serotype 19A) isolates were MDR.

Each of the two ST15 isolates expressing serotype 14 capsules displayed similar antibiograms to each other i.e. intermediate resistance to penicillin and high level resistance to erythromycin (Table 3). High level penicillin resistance was evident for a ST156 isolate expressing a serotype 14 capsule, while another representative of the same clone expressing a serotype 9V capsule, demonstrated intermediate resistance to penicillin. Two STs were associated with serotype 6A, ST65 (n=2) and ST2611 (n=1), each from a different clonal lineage and susceptible to all antimicrobials tested.

PMEN clones.

Over half (54%) of the isolates in this study were representatives of PMEN clones or their variants while 40% of these demonstrated reduced susceptibility to at least one antimicrobial agent (Table 2; Table 3). The most common PMEN clone was Netherlands 3 -31 (n=4) followed by Spain 9V -3, a SLV of Denmark 14 -32 and a SLV of England 14 -9 with two representative isolates each. Other clones included Sweden 4 -38, a SLV of Tennessee 23F -4, a SLV of Tennessee 14 -18, and a SLV and DLV of Netherlands 15B -37.

DISCUSSION.

This is the first report of the serotype distribution of S. pneumoniae causing otitis media in children in Ireland. In order to determine the serotypes associated with these infections, a combination of multiplex PCR and traditional serological methods was employed. During the multiplex PCR validation process using control collections of isolates, no cross reactivity with organisms other than S. pneumoniae was identified and the serogroup/serotype of all S. pneumoniae isolates targeted by the scheme was accurately determined, although three exceptions were initially identified. Nonetheless, further analysis confirmed the PCR result as accurate.

For the control collection of isolates, associated with invasive infections in Ireland and whose serotype was previously determined by the Quellung reaction [START_REF] Bennett | Penicillin susceptibility and epidemiological typing of invasive pneumococcal isolates in the Republic of Ireland[END_REF], 74% (97/132) of the serogroups/serotypes were identified using the multiplex PCR scheme.

Due to For this reason and also the less laborious nature of the PCR method, in addition to the absence of a DNA purification step, it was considered feasible and more efficient to employ this method in combination with serological methods for serotype determination of isolates in our laboratory setting. A limitation of this PCR method is its inability to differentiate between the commonly occurring serotype 6A, 6B and the more recently described 6C [START_REF] Park | Discovery of a new capsular serotype (6C) within serogroup 6 of Streptococcus pneumoniae[END_REF]. In addition, this method does not identify serotype 19A, which is a frequently occurring serotype among OM isolates in this study. For these and the other serotypes not targeted with the multiplex PCR assays, confirmation to serotype level using traditional serological methods and/or PCR in the case of 6C, is required.

The WHO estimates that complications of OM result in approximately 51000 deaths in children younger than five years of age every year. This infection is also one of the most common reasons for sick children to visit their general practitioner/ physician and a major reason for antibiotic prescriptions, especially in the case of recurrent OM, whereby 40% will have six or more recurrences by seven years of age [START_REF] Vergison | Otitis media and its consequences: beyond the earache[END_REF]. Efforts are being made in a number of countries to reduce antibiotic consumption, as levels of antimicrobial resistance are developing more rapidly than the development of new antimicrobials. However, use of a pneumococcal conjugate vaccine reduced otitis associated antimicrobial consumption by 20% [178]. Of the four serotypes associated with reduced susceptibility to penicillin in this study, three of these are contained in PCV7 and PHiD-CV, while all are contained in PCV13. Reduced susceptibility to penicillin (25% of isolates) was more common than resistance to erythromycin (18%) and clindamycin (11%). This level of penicillin resistance is half that observed (50%) for OM isolates in Spain in the 1990's, where OM clinical failures were associated with 90% penicillin resistance [189]. Conversely, penicillin resistance in Ireland is considerably higher than the 1.5% level observed for S. pneumoniae isolates in the Netherlands [1920], a country with one of the lowest rates of antimicrobial prescribing in the world [201].

In the present study, the top five ranked serotypes were 3, 19F, 19A, 14 and 6A, which are also the most common causes of OM worldwide [212]. The potential coverage provided by PCV7 and PHiD-CV was similar at 43% and 46%, respectively. This increased to 86% for PCV13. In the US since PCV7 introduction a serotype 3 clone, namely, Netherlands 3 -31(ST180) has become more prevalent and is emerging as a common vaccine replacement serotype causing disease [START_REF] Abdelnour | Streptococcus pneumoniae serotype 3 among Costa Rican children with otitis media: clinical, epidemiological characteristics and antimicrobial resistance patterns[END_REF]212,223]. The present study indicates that this particular serotype 3 clone, susceptible to all antibiotics tested, was most commonly associated with OM in Dublin prior to PCV7 vaccine introduction. This finding is similar to reports from Costa Rica, Israel and Lithuania demonstrating serotype Overall, the phenotype and genotype of S. pneumoniae isolates causing OM in Dublin is similar to that observed in a number of other countries. However, a limitation of the study is the small number of isolates available for characterisation over a relatively short period confined to one region in Ireland. For this reason results should be interpreted with some degree of caution. However, this study has provided important baseline information, for the most part, prior to PCV7 introduction, when serotype replacement due to vaccine introduction is unlikely to have occurred. It is interesting that three of the most common serotypes causing OM are not contained in PCV7, although cross-protection against serotype 6A may be provided, as evidenced by the decline of this serotype since PCV7 introduction in the US [323]. However, serotype 3, 6A and 19A are contained in PCV13, so infections caused by these serotypes are potentially vaccine preventable.

In conclusion, this study reinforces the importance of continued monitoring of the epidemiology of S. pneumoniae causing OM in Ireland, especially in the post-PCV7 vaccine era, in order to determine if serotype replacement occurs and to monitor antibiotic resistance. Such monitoring will guide future vaccine and treatment policies.

  isolated from swabs of spontaneous otological drainage in children with OM, received by the routine microbiology laboratory of the Children's University Hospital, Temple St., Dublin, an acute and tertiary referral paediatric hospital, were included in the study. Between January 2007 and December 2008, 28 S. pneumoniae isolates were identified by traditional microbiological methods following culture on Colombia blood agar plates. PCV7 was introduced to the immunization program in Ireland in September 2008, thus, three isolates were recovered after PCV7 introduction. Bacterial isolates control collection. A control collection (n=144) of S. pneumoniae isolates, whose serotype had been previously determined [5], were retrieved from storage in 1:3 glycerol (nutrient broth containing 0.5% glucose and 10% glycerol)-horse serum storage medium at -70 o C. Of these, 132/144 (92%) isolates had remained viable and were available for reserotyping.identified using standard microbiological methods. An additional control collection of 51 Statens Serum Institute (SSI) (Copenhagen, Denmark) S. pneumoniae strains were also utilised to validate the PCR assays. To ensure accurate S. pneumoniae identification, other organisms used to test the specificity of the assay included well characterised isolates of 21 bacterial species including S. agalactiae, S. pyogenes, Neisseria meningitidis and Haemophilus influenzae. Additional clinical isolates of viridans group streptococci (n=20) were also assessed to ensure no cross reactivity occurred.

PCR reaction .

 reaction Each 25µl PCR reaction mixture contained the following: 1x reaction buffer (Tris-Cl, KCl, (NH4) 2 SO 4 , 1.5mM MgCl 2 ; pH 8.7; Qiagen), 1mM MgCl 2 (Qiagen), 1x Q-solution (Qiagen), 200µM of each deoxynucleoside triphosphate (Promega Inc.), 2U Hot Star Taq DNA polymerase (Qiagen) and 1µl of lysed Microlysis ® Plus colony extract as template. Thermocycling was performed on a Peltier Thermal Cycler 200 DNA engine (MJ Research, Massachusetts, USA) using the following conditions: 95 O C for 15min, followed by 35 cycles at 95 o C for 40s, 55 o C for 45s, 65 o C for 150s and a final extension step at 65 O C for 1min. PCR products (10µl) were analysed using 3% (w/v) agarose gel electrophoresis with 5µg/ml ethidium bromide and visualised under UV transillumination. PCR product sizes were compared with the molecular size standard (100bp ladder; New England Biolabs, Beverley, MA).

3

  as a common cause of non-invasive infections [223]. Unlike the majority of S. pneumoniae serotypes, serotype 3 isolates have a mucoid appearance on colonial morphology. The generation of vaccine efficacy against serotype 3 OM was not demonstrated in a study involving a candidate 11-valent pneumococcal conjugate vaccine [23]. This may have been due to the atypical immune response generated or the high abundance of capsular polysaccharide expression against which antibodies may be less effective [234]. Nonetheless, the recently licensed PCV13 does contain serotype 3. The serotype 19A strains came from two different clonal lineages. ST276 (serotype 19A), which is a SLV of Denmark 14 -32 with a similar antibiogram to the clone that existed prior to PCV7 introduction in France, was also identified. ST276 has undergone clonal expansion in France to represent the predominant 19A replacement clone since PCV7 introduction. A number of OM treatment failures have been associated with serotype 19A [245]. Hence, close monitoring is necessary to determine whether or not the same level of expansion of this serotype occurs in Ireland [256]. The extensively characterised Spain 9V -3 (ST156) clone was present expressing a serotype 14 and 9V capsule. Spain 9V -3 is a widely disseminated clone commonly associated with penicillin resistance. It seems to have emerged in Spain in the 1980's and has established itself as a common cause of infections in a number of countries such as France, Denmark, Germany and the Netherlands [267-301]. The two serotype 6A clones, ST65 and ST2611, were also from diverse lineages and unlike the multidrug-resistant 6A clones in OM that are emerging in Israel [312], neither clone displayed resistance to the antimicrobials tested.

Figure legend Fig. 1

 legend1 Figure legend

  Fig.1

Table 1 .

 1 Validation results of Multiplex PCRs. Repeat serology confirmed serotype 14, 4 and 1 PCR results for isolates initially identified as 9V, 19F and 9V.

		Quellung	Multiplex PCR	Multiplex PCR
	Serotype*	Number	Number	cpsA	lytA
	14	28	29 a	+	+
	9V/(9A)	21	19 a	+	+
	4	12	13 a	+	+
	23F	4	4	+	+
	19F	4	3 a	+	+
	6A/6B	4	4	+	+
	18C/(18A,18B,18F) 3	3	+	+
	7F/(7A)	3	3	+	+
	3	9	9	+	+
	1	4	5 a	+	+
	5	5	5	+	+

*Rarely occurring serotypes not resolved by PCR method are indicated in parentheses. a

Table 2 .

 2 Properties and STs of S. pneumoniae isolates causing OM in Ireland

	ST			Allele number			Serotype(s) Total	PMEN clone
		aroE gdh gki recP spi xpt ddl	no.
									of
									isolates
	180 7	15 2	10	6	1	22	3	4	Netherlands 3 -31
	458 2	32 9	47	6	21 17	3	1
	73	2	13 2	15	6	1	1	19F	1
	87	5	5	7	7	8	5	4	19F	1
	420 15	38 19 5	36 20 6	19F	1
	425 1	5	1	5	1	1	8	19F	1
	3508 47	5	15 5	15 20 8	19F	1
	66	2	8	2	4	6	1	1	9N	1	SLV
											Tennessee 14 -18
	276 2	19 2	17	6	22 14	19A	2	SLV Denmark 14 -
											32
	2081 8	5	14 4	17 1	14	19A	1	DLV
											Netherlands 15B -
											37
	15	1	5	4	5	5	3	8	14	2	SLV England 14 -
											9
	156 7	11 10 1	6	8	1	14, 9V	2	Spain 9V -3
	65	2	7	4	10	10 1	27	6A	2

ACKNOWLEDGEMENTS

We are grateful to the staff of the Children's University Hospital, Temple St., Dublin for their assistance with the collection of isolates used in this study. We also wish to thank the Health Protection Surveillance Centre which financially supports the laboratory surveillance of pneumococcal disease in Ireland. We acknowledge the use of the pneumococcal MLST database, which is located at Imperial College London and funded by the Wellcome Trust.


Constantinos Makassikis

Stéphane Vialle

Xavier Warin

FT-GReLoSSS: a Skeletal-based approach towards application parallelization and low-overhead fault tolerance

Keywords: distributed fault tolerance, application-level checkpointing, SPMD paradigm, framework, skeletons

is a C++/MPI framework to ease the development of fault-tolerant parallel applications belonging to a SPMD family termed GReLoSSS. The originality of FTG is to rely on the MoLOToF programming model principles to facilitate the addition of an efficient checkpoint-based fault tolerance at the application level. Main features of MoLOToF encompass a structured application development based on fault-tolerant "skeletons" and lay emphasis on collaborations. The latter exist between the programmer, the framework and the underlying runtime middleware/environment. Together with the structured approach they contribute into achieving reduced checkpoint sizes, as well as reduced checkpoint and recovery overhead at runtime. This paper introduces the main principles of MoLOToF and the design of the FTG framework. To properly assess the framework's ease of use for a programmer as well as fault tolerance efficiency, a series of benchmarks were conducted up to 128 nodes on a multicore PC cluster. These benchmarks involved an existing parallel financial application for gas storage valuation, originally developed in collaboration with EDF company, and a rewritten version which made use of the FTG framework and its features. Experiments results display low-overhead compared to existing system-level counterparts.

FT-GReLoSSS: approche fondée sur des squelettes pour la parallélisation d'application et la tolérance aux pannes à faible coût Résumé : FT-GReLoSSS (FTG) est un framework C++/MPI pour faciliter le développement d'applications parallèles tolérantes aux pannes et appartenant à une famille d'algorithmes SPMD nommée GReLoSSS. L'originalité de FTG est de reposer sur les principes du modèle de programmation MoLOToF pour faciliter l'ajout d'une tolérance aux pannes de niveau applicatif fondée sur la réalisation de points de reprise (checkpointing). MoLOToF se caractérise d'une part par un développement structuré d'applications fondé sur des "squelettes" tolérants aux pannes, et d'autre part, par l'utilisation de collaborations. Ces dernières existent entre le programmeur, le framework et l'intergiciel d'exécution/environnement sous-jacents. Couplées à l'approche structurée, les collaborations contribuent à obtenir des tailles de points de reprise réduites ainsi que des surcoûts de checkpointing et de reprise réduits. Ce rapport introduit les principaux principes de MoLOToF ainsi que la conception du framework FTG. Pour évaluer la facilité d'utilisation du framework et l'efficacité de la tolérance aux pannes, une série de tests ont été menés jusqu'à 128 noeuds d'une grappe de PCs multi-coeurs. Ces tests impliquaient une application financière existante de valorisation d'actifs de stockage de gaz, initialement développée en collaboration avec EDF, et qui a été adaptée pour pouvoir profiter de FTG et de ses fonctionnalités. Les résultats expérimentaux témoignent de surcoûts faibles par rapport à des solutions équivalentes de niveau système.

Mots-clés : tolérance aux pannes distribuée, points de reprise applicatifs, paradigme SPMD, framework, squelettes

Introduction

As a result of increased competition, many industries have strived for new, more accurate simulation models. Their design, assessment and exploitation often requires huge amounts of computational power which eventually led those industries to adopt high performance distributed systems (HPDS).

In many cases, computations may be broken into several independent parts which can straightforwardly exploit HPDS. However, not all computations fit efficiently the embarrassingly parallel computation model. Programmers with little parallel background face three major issues. Firstly, they need to learn how to make their applications run efficiently on HPDS. Secondly, as the probability of failure increases with the number of computational resources (or nodes) used, developped applications have to be robust and fault-tolerant. Finally, HPDS systems involve rapidly evolving hardware and software which compels the use of portable fault tolerance (FT) solutions. Efficiency of these solutions is not an option as industrial environments are subject to time constraints: for example, some long-running applications in financial institutions are run overnight, and their results are expected in the morning in order to decide on the daily strategy to follow. Hence, FT efficiency translates into little slowdown of applications during failure-free time intervals. In case of failure, it translates into little wasted work and short restart times.

Checkpointing is a widely spread FT technique which consists in periodically saving the application state to stable storage. Following a failure, application execution is interrupted and the most recent saved state is then used to resume execution. This paper presents FTG, a specialized framework derived from the MoLOToF programming model [START_REF] Makassikis | A Skeletal-Based Approach for the Development of Fault-Tolerant SPMD Applications[END_REF]. MoLOToF aims to facilitate the development of efficient parallel applications and their endowment with efficient application-level checkpointing.

After covering related works (Section 2), we describe the MoLOToF programming model and the software architecture of the FTG programming framework (Sections 3 and 4). Ease of development is assessed through the use of FTG on an existing industrial application for gas storage valuation from EDF company (Section 5). Finally, the efficiency of the approach is evaluated experimentally (Section 6) before concluding (Section 7).

Related Works

Existing FT research for message passing applications has focused a lot on transparency. By combining a sequential checkpointer (e.g.: BLCR [START_REF] Hargrove | Berkeley Lab Checkpoint/Restart (BLCR) for Linux Clusters[END_REF], MTCP [START_REF] Rieker | Transparent User-Level Checkpointing for the Native Posix Thread Library for Linux[END_REF] . . . ) and some rollback recovery (RR) protocol, it was possible to endow existing MPI libraries such as OpenMPI [START_REF] Hursey | The Design and Implementation of Checkpoint/Restart Process Fault Tolerance for Open MPI[END_REF] or MPICH [START_REF] Bouteiller | MPICH-V: a Multiprotocol Fault Tolerant MPI[END_REF] with transparent FT or to provide more general solutions such as DMTCP [START_REF] Ansel | DMTCP: Transparent checkpointing for cluster computations and the desktop[END_REF] for socket-based distributed applications. RR protocols ensure that individual checkpoints of MPI processes remain consistent despite interdependencies created by communications: they form a recovery line. Most available solutions implement a so called blocking RR protocol which "freezes" communications while a checkpoint is made. MPICH-V [START_REF] Bouteiller | MPICH-V: a Multiprotocol Fault Tolerant MPI[END_REF] is the exception, for it implements several other RR protocols. From our experience with the OpenMPI-BLCR pair or DMTCP, such solutions yield bulky checkpoint files, and are very sensitive to changes of the underlying runtime system. As a result of including too much system-dependent information, issues are often raised as the system is updated.

C 3 [START_REF] Bronevetsky | Recent Advances in Checkpoint/Recovery Systems[END_REF] and CPPC [START_REF] Rodríguez | CPPC: A compiler-assisted tool for portable checkpointing of message-passing applications[END_REF] strive for similar levels of transparency to the user, but at the application level. By leveraging source-to-source compilers, source code of C/Fortran MPI applications is automatically transformed such as the resulting application can checkpoint and restart itself. Besides unburdening the user from non-trivial source code transformations, this approach benefits from high portability for it works at the application level. Many legacy and recent applications written in C/Fortran can benefit from this approach. However, providing similar transformations is challenging for languages such as C++ which is commonly used in many applications.

On the other end, FT can be tackled manually. Though not clearly stated in the litterature, the manual approach to FT is not unusual. Admittedly, the approach is tedious even with the aid of third-party libraries (e.g.: TCS [START_REF] Stone | A Checkpoint and Recovery System for the Pittsburgh Supercomputing Center Terascale Computing System[END_REF], SCR [START_REF] Moody | Design, Modeling, and Evaluation of a Scalable Multi-level Checkpointing System[END_REF] checkpointing libraries). Nevertheless, resulting FT is very efficient and portable, and the user may improve it further based on his knowledge of the application. For production applications, the efficiency gained quickly outweighs the endeavour.

Finally, users can rely upon frameworks to achieve FT. PUL-RD [START_REF] Silva | Fault-Tolerance on Regular Decomposition Grid Applications[END_REF] and Cactus [START_REF] Allen | Cactus Tools for Grid Applications[END_REF] are examples of this approach where users are subject to some programming constraints in exchange of which they benefit among others from FT. Compared to the previous two approaches, frameworks also facilitate writing parallel applications. As a framework, FTG shares similarities with PUL-RD and Cactus as far as the parallelization model is concerned. All three involve iterative calculations with two array datastructures, swapped at the end of each iteration. FTG supports applications with none-trivial distribution of calculations among processes. It differs by relying on a specific programming model for fault tolerance named MoLOToF. The introduction of fault-tolerant skeletons yields an explicit structuration of the application, which in turn provides a simple, yet efficient, way to endow applications with FT. Hence, users are led to "actively" interact with the framework. Our fault-tolerant skeletons are inspired by the skeleton programming approach, and are designed at lower level. Hence they should be named "sub-skeletons" (or "low-level skeletons"). But for simplicity, in the rest of this paper we call them "skeletons".

3 MoLOToF programming model for FT MoLOToF (Model for Low-Overhead Tolerance of Faults) is a programming model geared towards easing the development of fault-tolerant parallel applications. MoLOToF relies (1) on a peculiar structuring of the application and (2) on establishing collaborations through interaction functionalities between the programmer, the framework and the underlying environment. To achieve these, MoLOToF introduces the concept of fault-tolerant skeletons. The latter are structured pieces of code endowed with fault-tolerant properties. Usually, they are made of fault-tolerant loops each of which has the ability to save and restore itself. As illustrated in fig. 1, a typical loop body contains calculation, and possibly communication phases (l. [START_REF] Bouteiller | MPICH-V: a Multiprotocol Fault Tolerant MPI[END_REF][START_REF] Ansel | DMTCP: Transparent checkpointing for cluster computations and the desktop[END_REF]. Checkpoint definition (l. 7) within the loop allows to save calculations and related application state.

Skeleton-based code structuring

Using fault-tolerant skeletons, the application is split into two broad types of code operations: heavy and light operations. Heavy operations correspond to time-consuming code. Such code is usually found within calculation and communication phases of a skeleton. Light operations designate every other piece of code which is really fast to reexecute, and hence not interesting to checkpoint. Such separation results in a straightforward save/restore (checkpointing) mechanism based on application reexecution. To illustrate the mechanics, let us assume an application made of a single skeleton as the one in fig. 1. As represented in fig. 2, the application comprises an initialization phase, a skeleton and a cleanup phase. During normal application execution (cf. Fig. 3), checkpointing occurs whenever a checkpoint location and a checkpoint condition are met. In our example, a checkpoint occurs on the second iteration. Among saved data is the iterator's value (i.e.: 2). Consistency of individual checkpoints is ensured by some adapted RR protocol. An application executes until it quits gracefully or until a failure occurs. In our example it fails somewhere during communications of the third iteration.

Source code Init Calc Comm Ckpt Cleanup [ [ i: 1..n

Runtime

Calc

Comm Ckpt After the detection of a failure (cf. Fig. 4), the application is restarted from its most recent recovery line: here, the checkpoint taken at iteration 2. The application enters recovery mode (cf. "restarting:true" in fig. 4). Namely, each process restarts from the very beginning as in its initial run til it reaches the checkpoint location where the checkpoint it is supposed to restart from was achieved. During this course, only light operations are reexecuted therefore resulting in fast restarts. When reaching the appropriate checkpoint location, application context is restored. To complete recovery, checkpoint file contents are loaded back into the application, which can then normally resume its execution: the application falls back into normal execution mode (cf. "restarting:false" in fig. 4) and resumes with the third iteration.

Collaborations

To ease parallel programming, skeletons enclose a parallelization paradigm and come usually as part of a specialized framework. Through their use, the programmer easily develops his parallel application and installs checkpoint/restart semantics in it. Further involvement of the user may be required when checkpoints of fault-tolerant skeletons are not "self-contained". This means that their default contents do not capture entirely the application state. The user has to register as part of the checkpoint contents missing data in order to capture the application state. When fault-tolerant skeletons are self-contained, such action on behalf of the user is unnecessary as far as correctness is concerned. However, it is desirable in order to reduce the resulting checkpoint size and hence improve checkpointing efficiency. Therefore, in order to either constitute a consistent checkpoint or improve efficiency, the programmer's intervention may be needed. It makes all the more sense as the programmer, knowing the application semantics, may come up with smart choices.

Moreover, checkpoint efficiency may depend on the execution environment characteristics over time. For example, platforms under heavy load or aged platforms are more prone to failures. Hence, it is important that the application is able to adapt its fault tolerance based on external information.

To allow such interventions, FT skeletons come up with additional functionality allowing the programmer to control data which is enclosed in a checkpoint as well as setting the checkpointing frequency. Such interaction is seen as collaboration between the framework and the programmer. The other collaboration envisionned by MoLOToF lies between the framework and the underlying middleware, and is fully compatible with the spirit of fault-tolerant ecosystems such as the Fault-Tolerant Backplane [START_REF] Gupta | CIFTS: A Coordinated infrastructure for Fault-Tolerant Systems[END_REF]. To date, two framework implementations following the MoLOToF programming model exist. The first implementation is a Javaspace-based Master-Worker framework [START_REF] Galtier | A Javaspace-based Framework for Efficient Fault-Tolerant Master-Worker Distributed Applications[END_REF] which considers self-contained skeletons. The second implementation is a framework named FT-SPMD [START_REF] Makassikis | A Skeletal-Based Approach for the Development of Fault-Tolerant SPMD Applications[END_REF] which targets a broad family of SPMD applications named GReLoSSS (cf. Section 4). FT-SPMD benefitted from several modifications in its architecture and API which resulted in a seriously improved version called FT-GReLoSSS (FTG for short) which is presented hereafter. Compared to [START_REF] Makassikis | A Skeletal-Based Approach for the Development of Fault-Tolerant SPMD Applications[END_REF] the framework architecture is presented more in-depth and it is validated using an industrial application (cf. Section 4).

FTG programming framework 4.1 GReLoSSS computation model

The GReLoSSS (Globally Relaxed, Locally Strict Synchronization SPMD) parallel computation model encompasses SPMD applications consisting in a main loop where each iteration (or superstep) contains a calculation and a communication phase. Such applications follow globally a classic BSP scheme [START_REF] Valiant | A bridging model for parallel computation[END_REF]. However, these applications differ from BSP in two ways:

1. to improve efficiency on big parallel architectures, algorithms do not have a global synchronization between supersteps: each process starts its next superstep as soon as it has issued, but not necessarily completed, all its communications. Hence, global synchronization is relaxed as in the PRO model [START_REF] Gebremedhin | PRO: A Model for Parallel Resource-Optimal Computation[END_REF].

2. yet, to improve checkpointing efficiency, each process completes all its communications before starting the next superstep. Dependencies due to communications between processes disappear and consistent checkpointing is facilitated. Hence, compared to PRO, GReLoSSS has a locally strict synchronization which prevents overlapping of calculations with communications.

Application class supported by GReLoSSS and examples

The applications of the class supported by FTG (FT-GReLoSSS) involve two array datastructures one of which contains data used in the current superstep. The role of the second one depends on the application. A domain decomposition application such as Jacobi relaxations [START_REF] Gropp | Using MPI: Portable Parallel Programming with the Message-Passing Interface[END_REF] will use that datastructure to store new results at the current superstep. A data circulation application such as in some parallel matrix multiplication algorithms [START_REF] Makassikis | A Skeletal-Based Approach for the Development of Fault-Tolerant SPMD Applications[END_REF] will use it to receive data for the next superstep from neighbouring processes.

During communications, processes exchange data corresponding to initial data or to intermediary results. Communications between processes may be quite varied, yet, in most cases, they are sufficiently foreseeable to be specified by the programmer. For example, in a Jacobi relaxations application, it consists in exchanging borders (or shadow regions) between subdomains assigned to different processes. A subdomain designates the subset of the entire domain which was assigned to a process. Thus, a subdomain has the same datatype as array datastructures. At the end of communications, the two array datastructures are swapped.

Algorithms using two array datastructures (of 1 or more dimensions) such as the ones described may seem restrictive. But from our experience, they cover the needs of a fairly wide class of applications.

Software architecture and features

FTG is a C++ framework implementing programming principles described in MoLOToF. It is built on top of the MPI specification which makes it compatible with every MPI library. It provides a set of classes to ease the development of fault-tolerant GReLoSSS parallel applications.

Relations between main FTG classes are depicted in fig. 5. The ftg_Skel class represents a GReLoSSS skeleton and comprises:

• a calculation kernel (cf. ftg_Calc_Kern) which provides the calculation method and the skeleton's main loop iterator. It is also possible to specify a condition for communications achievement.

• a routing plan (cf. ftg_Routing_plan) which is responsible for determining and scheduling communication exchanges between processes.

• a checkpoint (cf. ftg_Checkpoint) which sets the location in the skeleton where checkpoints will be taken.

• two array datastructures (cf. ftg_Domain) as introduced previously in the GReLoSSS computation model (cf. Sections 4.1 and 4.2). The source code excerpt of ftg_Skel class (cf. Fig. 6) shows how these abstractions are layed out within the execute method. The resulting structure closely follows the example skeleton in fig. 1. It differs in the existence of conditional communications (l. 50-51) and the presence of a swap method (l. 61) where the two datastructures exchange their data according to the GReLoSSS computation model. Condition has to be verified globally: all processes enter communication phase or none. Notice the name prefixes: ftgu_ are user-provided (l. 29-30, 41, 44), ftg_ are framework provided (l. [START_REF] Gebremedhin | PRO: A Model for Parallel Resource-Optimal Computation[END_REF]), and ftgf_ are framework-provided based on user-provided methods (l. 38).

Concerning FT, the GReLoSSS skeleton defines a checkpoint right after the swap operation. By default, it contains data internal to the skeleton such as the two datastructures and skeleton iterators. ftg_Skel exposes to the programmer an interface allowing him to control data to include in the checkpoint. Furthermore, data which can be included in checkpoint files has to be either a primitive C++ type or derived from the framework's ftg_Serializable class (cf. Fig. 5), in which case the programmer has to provide a pair of save and load methods. This approach was used within FTG to make ftg_Domain and ftg_Domain_desc classes serializable. The latter is a helper class which describes the extent of a domain (or subdomain) along each dimension.

Upon instantiation of the skeleton, the user has to inform the routing plan whether application communication consists in borders exchanges or in circulating data. Aside from this, he may choose among different communication schemes. Currently, the routing plan integrates two schemes based on the MPI_Issend-MPI_Irecv pair of MPI primitives. The first scheme performs all communications in parallel, while the second tries to "pace" them so that there is only a limited amount of ongoing concurrent communications at a time [START_REF] Makassikis | Distribution of a Stochastic Control Algorithm Applied to Gas Storage Valuation[END_REF]. Concerning data movement, the routing plan groups sparse data into contiguous large messages and also minimizes copies from/to communication buffers by detecting ranges of contiguous elements. This last feature proved useful in applications where data to transfer is contiguous as it avoids using intermediate buffers.

The calculation kernel is merely made of pure virtual methods which the programmer has to define.

The array datastructure provides the programmer with an interface to N-dimensional arrays enclosed in ftg_Domain. When interfacing his own array datastructure (by inheriting from ftg_Domain), the programmer provides information regarding (1) the way the domain is split among processes as well as (2) the way to access an element given its coordinates (x 1 , x 2 , • • • , x n ). Moreover, the programmer provides information regarding the storage order. In 2D, C storage order consists in storing the array line-wise, while a Fortran storage order consists in storing the array column-wise. Storage order generalizes in higher dimensions and can lead to more different storages depending on the order of storing each dimension.

Since the interface proposed by ftg_Domain targets N -dimensional domains, it manipulates a vector of coordinates which can be inconvenient to the programmer in small dimensions. It can also be less efficient. Therefore, FTG proposes specific interfaces for dimensions 2 and 3 (resp. ftg_Domain2D and ftg_Domain3D in fig. 5).

Finally, after splitting a domain among processes (cf. Fig. 7), any given element can be located either relative to the first element of the entire domain or relative to the first element of the subdomain it belongs to (resp. Global ref and Local ref in fig. 7). The first are called global coordinates while the second are called local coordinates.

FTG allows the programmer to use either of them without any further involvement. Access through local coordinates is preferable since it is faster, but it is not always the most convenient. Using global coordinates incurs a slight overhead due to an additional offset computation.

From a design standpoint, the use of dynamic polymorphism (through classic inheritance) is sufficient. However, since array access methods are often used, the resulting overhead is tremendous. That is why FTG relies also on static polymorphism, and more specifically on the Curiously Recurring Template Pattern (CRTP) [START_REF] Coplien | Curiously recurring template patterns[END_REF].

In the following section, we decribe the Swing application and the steps to interface it efficiently with FTG. 5 Swing financial application migration

t g _ C h e c k p o i n t _ t c ; / / S k e l e t o n ' s c h e c k p o i n t b o o l i s _ c i r c u l a t i o n ; / / C o m mu n i c at i o n t y p e p u b l i c : / / C o n s t r u c t o r , D e s t r u c t o r and o t h e r m e t h

Due to house heating, demand in gas is higher in winter than in summer. Gas is mainly provided by pipes that cannot deliver more than a specific amount of gas per day and it results in far higher gas prices in winter. In order to provide energy to their customers, gas market actors have to own some storages that can help them smooth peak demand.

In order to assess a project's rentability, gas companies can use a financial real option approach that can be implemented in a software.

Swing application goal and implementation

The Swing application is used at EDF company to valuate a gas storage facing the energy market. It aims at giving the average cash flow generated by the asset depending on some prices scenario. It also gives the management and hedging strategies [START_REF] Warin | Gas storage hedging[END_REF]. These calculated strategies are used in a second application simulating the way the storage is used. This second application gives some cash flow distributions obtained by the storage management. While the time needed for the second application is short, the swing valuation can be very costly depending on the price models used to generate scenarios. The price models used at EDF for this software are a gaussian one-factor model (g), a normal inverse gaussian model (nig), and a two-factor gaussian model (g2d) [START_REF] Warin | Electricity asset management with future hedging[END_REF].

The resolution method for the Swing application is the dynamic programming method that has been written in C++ and uses MPI for parallelization following the methodology in [START_REF] Makassikis | Distribution of a Stochastic Control Algorithm Applied to Gas Storage Valuation[END_REF]. It also makes extensive use of the Blitz library for its convenient array manipulation facilities [START_REF]The Blitz++ library[END_REF], and is about 18380 logical lines of code.

From an algorithmic standpoint, the Swing application fits perfectly the GReLoSSS computation model. While being a domain decomposition application, it is more complex than classic Jacobi relaxations. Indeed, over supersteps:

• calculations involve a subdomain of the entire calculation domain; that subdomain may change;

• shadow regions between processes have no fixed extent and may change as well.

To interface the existing application with FTG, we follow some typical development steps which are described hereafter and summarized in fig. 8.

Development workflow

As part of the typical development steps, we have to inherit from ftg_Calc_Kern to define a calculation kernel: t e m p l a t e < c l a s s T_numtype , i n t N_rank , c l a s s T _ i t e r > c l a s s Swing :

p u b l i c f t g _ C a l c _ K e r n e l <T_numtype , N_rank , T _ i t e r , Swing_Domain > { } and also to inherit from one of the available ftg_Domain classes to define the calculation domain: t e m p l a t e < c l a s s T_numtype , i n t N_rank , c l a s s T _ i t e r > c l a s s Swing_Domain :

p u b l i c ftg_Domain2D <T_numtype , N_rank , T _ i t e r , When inheriting from ftg_Calc_Kern we provide the calculation method:

Swing_Domain > / / i m p l e m e n t s C u r i o u s l y R e c u r r i n g T e m p l a t e P a t t e r n { } RR n°7797

v o i d f t g u _ c a l c u l a t e ( f t g u _ D o m a i n _ t * d1 , f t g u _ D o m a i n _ t * d2 , T _ i t e r s t e p ) ;
the main loop calculation iterator:

T _ i t e r f t g u _ b e g ( v o i d ) ; T _ i t e r f t g u _ e n d ( v o i d ) ; T _ i t e r f t g u _ n x t ( T _ i t e r ) ;

and a method which tells when to achieve communications: When inheriting from one of the available ftg_Domain classes, we have to define two so called partition methods: which tell what data is owned and what is needed by each process at each superstep.

f t g _ D o m a i n _ d e s c <N_rank > f t g u _
In an application with borders exchange, such as the Swing application, data owned is a subset of the data needed. Moreover, we have to define a method telling the framework how to resize a domain:

v o i d
f t g u _ r e s i z e ( T i n y V e c t o r < i n t , N_rank > &e x t e n t )

and some methods on how to access any element or retrieve its memory address, given its local coordinates:

d o u b l e f t g u _ l g e t ( i n t l x , i n t l y ) v o i d f t g u _ l s e t ( i n t l x , i n t l y , d o u b l e e ) d o u b l e * f t g u _ l g e t A d d r ( T i n y V e c t o r < i n t , N_rank > &l c o o r d )
The resize method is not relevant to every application. Yet, in the case of Swing, it is interesting as subdomains attributed to each process and subdomains sizes vary across the course of an execution. Leveraging this feature, FTG can optimize memory management. The resize, the partition and the local access methods which the user provides (ftgu_ prefix) are used to build the routing plan. The "ftg_Domain * " classes also use local access methods to implement the global access methods they provide the user with.

Finally, the application's main function consisted in instantiating a calculation kernel and using it to initialize and execute the GReLoSSS skeleton. The corresponding source code has to be enclosed between the initialization and finalization statements of the framework provided by the ftg_Cpr (Checkpoint recovery) singleton class (cf. Fig. 8). The last step consists in choosing relevant checkpoint data. In particular, we unregister from the checkpoint the Domain which is useless at the end of the iteration, and register some application specific data which is not restored upon recovery (by mere reexecution), and hence has to be saved.

Ease of development assessement

The definition of the calculation method within our Swing calculation kernel involves the reuse of calculation functions in Gas_storage_price, and access to data in Gas_storage_asset (cf. Fig. 8). Both classes stem from the original Swing application. Similarly, partition methods in our Swing_Domain class require access to data and code reuse from Gas_storage_asset.

The calculation and partition methods are usually the most time-consuming to write. Especially the partition methods, which are prone to error due to indexes. Since compatible partition functions already exist in the original parallel Swing application, it is a matter of adapting them to use ftg_Domain_desc type to describe subdomains. The same applies for the calculation method, but some pieces of code have to be rewritten to use the accessors provided by ftg_Domain instead of those provided by Blitz Arrays. In the process, the interface of Swing_Domain was enriched to allow efficient and convenient access to ranges of data.

Another challenge we encountered concerned efficiency. In particular, partition methods have to be really fast as they are often called by the framework either to access elements through global coordinates or to build the routing plan.

In the end, the full fault-tolerant Swing application with FTG is about 353 logical lines of code less than the original application. Provided we are careful, the resulting application displays a clean and simple design as shown in fig. 8. Moreover, the application is fault-tolerant and the programmer has less lines of code to write. As shown in the next section, performance is as good as without FTG, sometimes slightly better.

Performance experiments

The evaluation of FTG consists in a comparison with the popular OpenMPI (OMPI) library which implements a blocking checkpoint protocol in combination with BLCR. OMPI applications benefit from a system-level checkpoint/restart solution. Experiments led, assess FTG and OMPI performances without and with FT enabled. When FT is enabled, checkpoints are taken, and we consider application runtimes in the failure-free case and recovery times otherwise. We also consider checkpoint size reduction and their impact on runtime overhead.

In these experiments, we consider temporary crash failures of nodes which can be dealt with by having both FTG and OMPI store checkpoints locally to each node (i.e.: /tmp). Also, this setting avoids disrupting measurements due to additional communications resulting from significant data transfer across the network.

Testbed description

Runtime overheads comparison in absence of checkpoints

Fig. 9 shows that FTG exhibits negligible to no overhead for the g2d price model, and scales well up to 256 nodes where it performs better than the application not using FTG. Similar results were observed for the g and nig price models.

Checkpoint size

RR n°7797

Checkpoint size is dominated by the size of the two array datastructures used for parallelization. Since at the end of a superstep, the input data is useless, it can be omitted from the checkpoint. This optimization results in checkpoint sizes per node twice as light compared to system-level checkpoints (cf. Fig. 10). The same can be observed for the g and nig models. These models require less memory and result in checkpoints size ranging between 20 MB and several kilobytes with FTG compared to 40 MB et 7 MB for OMPI BLCR. In what follows, we focus on the g2d model: its' long runtimes make it a good candidate for FT.

Fault tolerance performance without failures

For this experiment, we have run both versions of our application with the g2d model and different number of checkpoints: the lengthier the benchmark the more checkpoints were achieved in order to minimize the wasted amount of time in case of failure. In the present experiments we set the maximum allowed wasted amount of time to 4 minutes. Fig. 11 reports the runtimes on 128, 64 and 32 nodes. The plots compare OMPI with BLCR and FTG. In the latter's case we measured runtimes without and with checkpoint size optimization in an attempt to quantify the impact of checkpoint size reduction. The difference remains marginal (0 -3%) and is the highest with the biggest checkpoint sizes as the ones involved on 32 nodes. The difference is expected to grow further with bulkier checkpoint files. The remaining overhead observed with OMPI-BLCR ranges between 6% and 40%. Overhead incurred by FTG is always lower and does not exceed 8% in all cases. OMPI-BLCR's overhead is mainly attributed to its blocking checkpoint protocol [START_REF] Hursey | The Design and Implementation of Checkpoint/Restart Process Fault Tolerance for Open MPI[END_REF], which, unlike FTG's protocol, involves communications and "freezes" the application execution in order to checkpoint.

Recovery overhead

Recovery from a checkpoint comprises (1) a negotiation phase where processes decide from which recovery line to recover, followed by a (2) context recovery phase, and finally (3) a recovery from checkpoint file phase (i.e.: time to load data). Measured negotiation phase time is negligible (< 10ms). Context recovery phase is small as well. Overall recovery time from a checkpoint is rather small as it does not exceed 1s for the g2d model on 128, 64 and 32 nodes. Thus, fast restarts combined with the low-overhead checkpointing pointed out earlier make FTG suited for applications with time constraints.

Conclusion and Perspectives

Endowing parallel applications with efficient checkpoint-based FT at the application level can be a tedious task which adds up to the existing difficulties of parallelization. Our approach based on the MoLOToF programming model introduces fault-tolerant skeletons and results in a tractable way for users to endow efficient fault tolerance into their applications. Moreover, combined with a specialized framework, MoLOToF eases parallel programming, and encourages a synergy between the user, the framework and the runtime environment to improve FT efficiency. MoLOToF is applied to the GReLoSSS family of applications which we characterized in this paper. The application of the resulting FTG framework to an industrial application of EDF company showed that initial development with FTG involved simple steps and points out some elements to watch in order to minimize runtime overhead. Finally, experiments show the effectiveness of the overall approach and especially the efficiency of FT. For the same number of checkpoints achieved, FTG yields smaller checkpoint sizes than OMPI-BLCR and incurs at most 8% runtime overhead. Recovery from achieved checkpoints exhibited negligible overheads. These results confirm previous ones we achieved on more elementary but varied benchmark applications [START_REF] Makassikis | A Skeletal-Based Approach for the Development of Fault-Tolerant SPMD Applications[END_REF]. Many principles of MoLOToF were used in the design of FTG. But some of them, such as the integration with fault-tolerant ecosystems, still have to be integrated and experimented with. Moreover, due to its inherent portability, the approach appears viable for hybrid GP-GPU applications. FTG might further be extended to support other parallelization models such as asynchronous distributed iterative algorithms. Future works are planned along the aforementioned lines.

5

  model for FT 3.1 Skeleton-based code structuring . . . . . . . . . . . . . . . . . . . . 3.2 Collaborations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .4 FTG programming framework 4.1 GReLoSSS computation model . . . . . . . . . . . . . . . . . . . . . 4.2 Application class supported by GReLoSSS and examples . . . . . . . 4.3 Software architecture and features . . . . . . . . . . . . . . . . . . . Swing financial application migration 5.1 Swing application goal and implementation . . . . . . . . . . . . . . 5.2 Development workflow . . . . . . . . . . . . . . . . . . . . . . . . . 5.3 Ease of development assessement . . . . . . . . . . . . . . . . . . . 6 Performance experiments 6.1 Testbed description . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.2 Runtime overheads comparison in absence of checkpoints . . . . . . . 6.3 Checkpoint size . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.4 Fault tolerance performance without failures . . . . . . . . . . . . . . 6.5 Recovery overhead . . . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 1 :

 1 Figure 1: MoLOToF fault-tolerant skeleton example.

Figure 2 :

 2 Figure 2: Application source code representation.

[FailureFigure 3 :

 3 Figure 3: Normal execution with checkpoint and failure.

Figure 4 :

 4 Figure 4: Recovery execution.

Figure 5 :

 5 Figure 5: FTG UML class diagram

Figure 6 :Figure 7 :

 67 Figure 6: FTG's fault-tolerant skeleton.

Figure 8 :

 8 Figure 8: Swing application within FTG.

b o o l f

  t g u _ d o _ e x e c u t e _ r o u t i n g _ p l a n ( T _ i t e r s t e p ) ;

Figure 9 :

 9 Figure 9: Runtime comparison under the g2d price model.

Figure 10 :

 10 Figure 10: Checkpoint sizes comparison under the g2d model.

Figure 11 :

 11 Figure 11: Runtimes comparison under the g2d model when increasing number of achieved checkpoints.

RR n°7797

Acknowledgment

This research is partially supported by Region Lorraine.

This research is partially supported by Region Lorraine Project Team § EDF R&D ¶ OSIRIS group