Ahmad Almhdie-Imjabber

Christophe Léger
email: christophe.leger@univ-orleans.fr

Local Comparison of Multimodal Medical Surfaces using a New Variational Splines Fitting Algorithm

Keywords: Local comparison, multimodal, surface fitting, surface reconstruction

In this paper, we propose a modified variational splines fitting (MVSF) algorithm for the reconstruction and comparison of medical surfaces. The MVSF algorithm presents the advantage to use high precise and low complex approximation of the equations of energy in the discrete temporal domain, compared to previously presented methods. It also takes into account the periodicity constraints encountered when reconstructing sphere-like shape closed surfaces. Moreover, it gives straightly the mean square error (MSE) between the original data and the reconstructed data, which is useful to quantify the approximation introduced by the reconstruction. The developed model has been successfully applied for real biomedical data; in particular for the reconstruction and comparison of the left ventricle of human heart, acquired by SPECT and ultrasound imaging modalities.

Introduction

Surface reconstruction is a step generally carried out before comparing or representing/ visualizing discrete data sets. It is essential when the two data sets to be compared present dissimilar numbers of samples, different sampling processes, or missing samples (sparse data). These cases are commonly encountered in medical imagery, and particularly in the comparison of data obtained from several imaging modalities [START_REF] Garcia | Left ventricle's surface reconstruction and volume estimation[END_REF]. Hence, it is necessary to carry out resamplings, using interpolation and/or smoothing methods.

Surface reconstruction and smoothing methods are widely used in practice to best estimate the original surface represented by a scattered noisy point. For instance, the number and the distribution of the initial samples of an organ surface lead to incomplete meshes. Stacks of two-dimensional contours can also be used to reconstruct three-dimensional surfaces [START_REF] Hoppe | Surface reconstruction from unorganized points[END_REF].

The problem of surface reconstruction can be solved using a variety of techniques. The finite element method is one of these techniques. An example of this method is the use of Voronoi diagram and Delaunay triangulation [START_REF] Gao | A Fast Algorithm for Delaunay based Surface Reconstruction[END_REF] to find the topological connection of the sample points. Fourier-based interpolation methods [START_REF] Bonciu | 4D reconstruction of the left ventricle during a single heart beat from ultrasound imaging[END_REF], to mention a few, are also used in the domain of surface reconstruction.

Another type of reconstructing smooth surfaces relies on the finite difference methods using deformable models [START_REF] Wang | Surface reconstruction using deformable models with interior and boundary constraints[END_REF] such as thin plate model [START_REF] Enciso | 3D Tooth Shape from Radiographs using Thin-Plate Splines[END_REF]. The main idea is to create an initial mesh and deform it to best match the range input data. This is usually expressed as an energy minimization problem. There are two forces that formulize the final shape of the reconstructed surface. One attracts the surface towards the input data and the other tries to keep the surface smooth. These methods are generally computationally inexpensive but comparably more sensitive to number of missing points. The former approach was selected for further study in this work as it is generally more appropriate for smooth surfaces such as surfaces of left ventricle of heart. In general, methods based on deformable models give nice smooth surfaces but are computationally expensive and complex. We present in this work a faster and less complex algorithm compared to a previous algorithm proposed by Richard Szeliski. The reconstruction step makes it easier to compare medical surfaces obtained by different modalities or at different time or position.

The next section concerns the description of the original VSF method and our proposed variant, required prior to local comparison of medical surfaces with incomplete samples. Data used in our experiments are described in section 3. Results are then presented and discussed. Finally, conclusions are given.

Variational Splines Fitting (VSF) Algorithm

In [START_REF] Szeliski | Fast surface interpolation using hierarchical basis functions[END_REF], the author proposed to use a deformable model (here we call it VSF) to estimate the missing points. In this model, the problem is formulated as an optimization one. The function to be minimized is written.

) x ( E ) x ( E ) x ( E s d    (1)
Where X = [X i,j ], (i = 0 : N-1, j = 0 : M-1) are the mesh regular points of the reconstructed surface, i and j indicate spatial positions. This function includes two constraints: the data compatibility constraint E d (X), and the smoothness constraint E s (X). The above formulation is usually expressed as an energy minimisation problem where an attracting force draws the mesh towards the sparse data « E d (x) » and a tension in the mesh keeps the surface smooth « E s (x) » [START_REF] Hoppe | Surface reconstruction from unorganized points[END_REF].  (>0) is the regularization parameter which is used to adjust the closeness of the fit between the surface and the sparse data set. Figure 1 shows an example of a 2D curve construction under the effect of the regularisation parameter. This parameter depends on the sparse data set and can be estimated using a generalised cross validation technique. In general, as  increases, the reconstructed surface becomes smoother. However, the probability of error between the original sparse samples and their corresponding estimated ones becomes higher. For very high values of  , (  ), the fitted surface tends towards a flat one. As  tends to zero, this probability becomes smaller, but the reconstructed surfaces might no longer be smooth. The data compatibility constraint measures the distance between the original sparse points and the interpolated smooth surface. The energy expression corresponding to the data compatibility constraint can be written as:

    i i i i i d d v u f w x 2 ) ) , ( ( 2
1 ) ( (2)
The discrete form is:

   ) , ( 2 , , , ) ( 2 1 )
, ( E

j i j i j i j i d d x w d x (3)
where x i,j stands for the discrete values of f(u i ,v i ) : x i,j = f(u i ,v i ), d i,j (d i,j = 0 at missing points) are the sparse samples of the original incomplete surface and the weights w i,j (w i,j = 0 at missing points) are inversely related to the variance of the measurements. The higher the weights, the better the reconstructed surface fits the original sparse samples. In this case, the data compatibility constraint over influences the overall energy function and the reconstructed surface become no longer smooth.

Reassembling all mesh points into a vector x, one can rewrite equations (3) in a matrix form. The energy corresponding to the data compatibility constraint becomes:

      d x A d x d x,    d T d 2 1 E (4)
where d is a zero-padded vector of data values and the diagonal matrix A d has entries w i at which the data coincide with the sparse data points and zeros elsewhere. In particular, this allows treating problems with missing or unknown data. Using the thin plate model, the energy function corresponding to the smoothness constraint can be written in continuous form as:

        dudv f f f f uv vv uu s 2 2 2 2 2 1 ) ( (5)
where f: f(u,v) is the smoothed continuous functional of the interpolated surface in u and v directions; and the subscripts uu, vv and uv indicate partial derivatives.

Original VSF algorithm

In the original VSF algorithm, the discrete form of the above energy function is derived using a classical finite-difference scheme.

The resulting discrete function is:

                                                  2 , , 1 1 , 1 , 1 2 2 1 , , 1 , 2 , 2 , 1 , , 1 2 2 2 2
) ( E v u j i j i j i j i v j i j i j i j i u j i j i j i v u s h h h h h h x x x x x x x x x x x (6)
where

u h u   and v h v  
are the step sizes of the regular mesh of the reconstructed surface in the u and v directions respectively.

Modified VSF algorithm (MVSF)

In this work, a more accurate derivation process is proposed to get the discrete energy corresponding to the above continuous function [START_REF] Wang | Surface reconstruction using deformable models with interior and boundary constraints[END_REF]. In particular, we propose to take into consideration higher element components of the corresponding two-variable Taylor formula. The discrete corresponding function is obtained by:

                                                      2 1 , 1 1 , 1 1 , 1 1 , 1 2 2 1 , , 1 , 2 , 2 , 1 , , 1 8 1
2 2 2 ) ( v u j i j i j i j i v j i j i j i j i u j i j i j i v u s h h h h h h x E x x x x x x x x x x (7)
We can see from the above expression that the first two terms are similar to the original formulation given in equation [START_REF] Enciso | 3D Tooth Shape from Radiographs using Thin-Plate Splines[END_REF]. They both use the so-called 5-starpoints scheme. The main advantage of the expression given in equation ( 7) is its accuracy in using the 4 diagonal points to approximate the crossed secondorder derivative at position (i,j), as shown in Figure 2a.

Obviously, Szeliski's approximation is a biased forward approximation, as shown in Figure 2b. The previous approximation is a biased approximation as presented in Figure 2a. Reassembling all mesh points into a vector X, we can rewrite equations ( 6) and [START_REF] Szeliski | Fast surface interpolation using hierarchical basis functions[END_REF] in a matrix form. The energy corresponding to the thin plate model can also be written in compact form as:

v u i j    → fuu,fvv  → fuv        v u i j       → fuu,fvv  → fuv     (
x A x x s T 2 1 ) ( E  (8)
where the stiffness matrix A s is a sparse matrix block diagonal matrix. The matrix A s has at most 13 non-zero entries per row, as discussed in [START_REF] Szeliski | Fast surface interpolation using hierarchical basis functions[END_REF]. The rows with the maximum number of entries are :

(n + 1)  M + (m + 2); n = 1 : N  4 et m = 1 : M  4; N et M > 4.
However, when periodicity is imposed upon the stiffness matrix A s , all rows would then have 13 non-zeros entries. Note that such periodicity constraint is appropriate when considering spherical coordinates used in describing human organs. The resulting stiffness matrix has a homogeneous structure and can be easily coded. It has been constructed following the method proposed by [START_REF] Sainsaulieu | Calcul Scientifique[END_REF].

If h u =h v =1, by imposing the continuity constraint, the stiffness matrix can be written in the following matrix form:

                           B C D 0 0 D C C B C D 0 0 D D C B C D 0 0 0 D C B C D 0 0 0 D C B C D D 0 0 D C B C C D 0 0 D C B A                8 1 s (9)
where B, C and D are M  N symmetric matrices. 0 is an MN zero matrix.

Discrete problem

Using equations ( 4), ( 6) and ( 7), the combined discrete energy expression can hence be written in matrix form as:

  c E T T    b x x A x x 2 1 (10)
where

A =   A s + A d , b = A d  d
and c is a constant that may be omitted in the minimization process. This energy function has a minimum at x = x * , which is the solution of the following linear system that is obtained via the Euler equation:

b x A *  (11)
since A is a strictly positive matrix. The above set of linear equations is, hence, a high dimension homogeneous positive definite system, that can be solved using the conjugate gradient method, thanks to the strict positivity of the matrix A. Our approximation is found to lead to better performance compared to that used in original VSF algorithm.

Data description

In this work, data consist of two sequences of surfaces of the left ventricle (LV) of the heart reconstructed within a cardiac cycle. Examinations were carried out on the same patient using successively two medical imaging modalities within a short period of time, in order to assume that the LV deformations are reproducible and hence medical comparisons would be applicable. The first sequence is composed of eight LV surfaces obtained after automatic segmentation from nuclear medicine imaging (NMI), known to be a "gold standard" examination for cardiac observation.

The second LV sequence is provided by a multidimensional ultrasound technique (US) called LV4D for Left Ventricle in 4 Dimensions [START_REF] Bonciu | 4D reconstruction of the left ventricle during a single heart beat from ultrasound imaging[END_REF] [9], previously proposed by the authors. The objective was to use the NMI examination to validate this ultrasound method. The closed incomplete LV surfaces are first developed to a rectangular domain of size (I,J), with I = [/2, +/2] and J = [0, 2]. The transformation between the two domains, spherical S(x,y,z) (Figure 3a) and rectangular D(, ) (Figure 3b), is obtained by using the following relations: x = cos() cos(), y = cos() sin(), z = sin(). The associated surface of the spherical representation is identical to that of the rectangular surface. The samples of the developed surface are duplicated in (,) and placed in a doubly periodic grid according to these two variables. When the number of samples is too high (or the sampling grid is of reduced size), several samples {ρi(i,i)} may correspond to the same element (i,j) in the grid. In this case, the average of these samples is taken at this position. By the same way, to reduce artefacts in =/2, the average of samples available for each of these values of  is assigned to all values of corresponding . In fact, the radius at the vertices of the closed surface are the same whatever the values of . In contrast, a sampling grid of large size, or a poor distribution of the initial samples of the reconstructed surface generates an incomplete grid.

Each missing value is then set to zero and its position is saved.

Results and discussion

Performance analysis

The results displayed in Figure 4 show that the proposed algorithm, MVSF, outperforms the original VSF algorithm, e.g. when the tolerance value of the iterative conjugate gradient algorithm is set to 10-4, the modified algorithm converges about 100 iterations faster than the original algorithm. One can also notice from Figure 4 that an inappropriate choice of the values of  and w can induce a higher number of iterations to get the same value of MSE between the reference surface and the reconstructed one. Author names and affiliations are -or three-column format, with their affiliations italicized and centered below their respective names. Include e-mail addresses if possible. Author information should to be centered beneath the title and printed in Times 12-point, nonboldface type. Multiple authors may be shown in a two be followed by two 12-point blank lines.

Comparison of LV surfaces

In order for the given LV surfaces to be locally compared, the irregularity of the data, the number of the available data, or even the absence of certain samples requires the use of an interpolation step. The resampling of data was implemented to visualize the local differences of the compared surfaces. It was produced by the MVSF algorithm, Figure 5c

Global comparison

The evolution of the NMI and US LV volumes as a function of time provides a global comparison of the reconstructed surfaces, as shown on Figure 6. To facilitate the comparisons, each volume is normalized relative to a maximum volume obtained in diastole. From these volumes, measurements are realised: estimation of the telediastolic (V D ) and telesystolic (V S ) moments and volumes, then evaluation of the ejection fraction:

D S D V V V FE   (12)
However, the lack of registration does not allow local comparisons, since the two surface absolute orientations are not known. The registration of the two sequences was performed by estimating the parameters of translation and rotation for each pair of corresponding surfaces of the two data sets using the CICP algorithm [START_REF] Almhdie | 3D Registration Using a New Implementation of the ICP Algorithm Based on A Comprehensive Lookup Matrix: Application to Medical Imaging[END_REF]. The sequence of surfaces obtained by scintigraphy served as a « reference », while the sequence of surfaces generated by ultrasound constructed the « target » data. It can be seen from this figure that the variation of evolution curves is comparable, even if a volume difference during diastolic phase perceptible. This slight difference is probably explained by the unequal precise information obtained during surface estimation. In fact, the acquisition and reconstruction modes are very different for both techniques. However, the determination of the fraction of ejection, current measurement performed during the cardiac examination, gives results of the same order of magnitude, and sufficient precision to make a diagnosis, according to cardiologists: F ej (US) = 0,624  F ej (MN) = 0,706.

In this case, the reported fractions of ejection indicate normal cardiac contractions (0,5  F ej 0,7). From the medical point of view, the results presented in this study show that the surface provided by the LV4D method are close to the reference gated SPECT surface. Actually, for all considered surfaces, the coding of colors remains in average around the homogeneous light blue -green -yellow colors, indicating low radial distance between the two LV surfaces. These results validate the ultrasound LV4D method as an alternative to the reference gated SPECT examination.

Local comparison

Conclusions and perspectives

The MVSF algorithm can be used to reconstruct closed surfaces of an organ of the human body, such as the left ventricle of the heart. It gives more precise surface reconstructions and is faster compared to its original VSF algorithm. However, the major drawback of both methods is the size of the rigidity matrices which increases in proportional to the square of the desired mesh. In fact, moving from a mesh of samples of 3232 to a mesh of 6464 results in moving from a dimension of stiffness matrix of 10241024 (~10 +6 ) elements to 40964096 (~1610 +6 ). The size of the handled matrices, for a relatively small number of samples (4096) compared to the possibilities offered today by the medical imaging equipment, is an obstacle to the use of these variational splines based methods.

There exist already multi-scale methods to address this problem. Concerning the structure matrix presented in this work, it seems interesting to investigate whether a direct method, which would not require intermediate excessive storage, could not be developed to solve the linear system of equation ( 11).

Another research perspective could be to develop a 3D version of the MVSF algorithm in order to eliminate the step of converting Cartesian coordinates into spherical coordinates to convert closed surfaces into developed surfaces. Finally, it would be interesting to implement automatic methods to estimate the parameters of considered algorithms.

Using a reliable extension of the classical splines smoothing SVF method, the work presented in this paper describes a comparison technique based on color representation of the radial distances between two corresponding surfaces, projected on the surface of reference. It has been successfully applied to highlight differences between corresponding LV surfaces obtained starting from ultrasound and gated SPECT examinations. Finally, it was shown that the original LV4D method provides reconstructions of the LV surfaces equivalent to the surfaces obtained from gated SPECT, "gold standard" examination in the field of dynamic cardiac reconstruction. This result has a great medical interest, since the new method could be used as a promising alternative to the gated SPECT examination, taking all the benefits of the ultrasound modality: non invasive method, no radioisotope exposition, no contraindication, cost less expensive, mobile system.

The study presented can be extended to measure the differences between closed or developed corresponding surfaces (images). Then, it can be used to compare different reconstructions of a unique surface, or to compare a unique reconstruction of different surfaces. Additional work must be carried out to confirm the preliminary promising results. Particularly, focus should be done to improve the accuracy of samples that characterize the observed surfaces

Figure 1 :

 1 Smoothing degree of reconstructed surfaces in function of  :  1 (a) >  2 (b).

Figure 2 :

 2 Point estimation using VSF (a) and MVSF (b).

Figure 3 :

 3 Transformation of a point from spherical domain (a) to rectangular domain (b).

Figure 4 :

 4 Figure 4 : Number of iterations necessary for the convergence of the two algorithms. Rate of convergence of the conjugate gradient algorithm.

Figure 5 :

 5 Figure 5 : Samples of a closed surface of a LV of heart (a). Reconstruction of surfaces by MVSF algorithm (b). Corresponding closed surface (c).

Figure 6 :

 6 Figure 6 : Evolution of the NMI and US LV volumes.

Figure

  Figure 7c presents an NM reconstructed surface of the left ventricle of heart using the proposed MVSF algorithm. A color coding allows cartograhing the distances between the surfaces obtained by the two modalities by visually local differences between the two registered surfaces.From the medical point of view, the results presented in this study show that the surface provided by the LV4D method are close to the reference gated SPECT surface. Actually, for all considered surfaces, the coding of colors remains in average around the homogeneous light blue -green -yellow colors, indicating low radial distance between the two LV surfaces. These results validate the ultrasound LV4D method as an alternative to the reference gated SPECT examination.

Figure 7 :

 7 Figure 7 : Local comparison of two LV registered surfaces (c), acquired from gated SPECT (in red color) and ultrasound (in brown color before registration (a) and in blue color after registration (b)).


A Almhdie
email: ahmad.almhdie@univ-orleans.fr

R Jennane
email: rachid.jennane@univ-orleans.fr

O Rozenbaum

E Lespessailles
email: eric.lespessailles@chr-orleans.fr

C Leger

Image processing for porous media characterization

Keywords: Classification, thinning, segmentation

In digital image processing, skeletonization is a valuable technique for the characterization of complex 3D porous media, such as bone, stone and soils. 3D thinning algorithms are usually used to extract one-voxel wide skeleton from 3D porous objects while preserving the topological information. Models based on simplified skeletons have been shown to be efficient in retrieving morphological information from large scale disordered objects at a local level. In this paper, we present a series of 3D skeleton-based image processing techniques for evaluating the micro-architecture of large scale disordered porous media. The proposed hybrid skeleton method combines curve and surface thinning methods with the help of an enhanced shape classification algorithm. Results on two different porous objects demonstrate the ability of the hybrid skeleton method to provide significant topological and morphological information.

Introduction

Most variants of 3D binary skeletons are based on curve thinning or surface thinning. Since curve thinning does not preserve the geometry of noncylindrical shapes, and surface thinning cannot erode enough beam shapes, neither curve nor surface thinning is suitable for the analysis of porous media. To overcome this problem, shape information should be integrated into the skeleton models. Hybrid skeleton has recently been proposed [START_REF] Aufort | A new 3D shape-dependant skeletonization method. Application to porous media[END_REF] to create structural models that take into account the shape of the object which enables a precise description of the local geometry. In this paper, different processing tools are presented and applied to the analysis of real porous objects. Different parameters are mathematically defined to demonstrate the efficiency of the hybrid skeleton technique for the characterization of 3D porous media.

Hybrid skeleton method

The hybrid skeleton algorithm is inspired from the work of Pothuaud [START_REF] Pothuaud | A new method for three-dimensional skeleton graph analysis of porous media[END_REF], but takes more accurately into account the shape of the object, in order to improve the geometrical approximation. It consists in classifying all object's voxels before thinning. To highlight the performance of the proposed method, bone and stone data sets at 16-µm and 3-µm of resolution, respectively, are used, as seen in figure 1(a). Before thinning, three pre-processing operations are used. First Filtering with a median filter to reduce noise. Then, binarisation using a threshold determined as the local minimum between the modes of the histogram of each 3D image. Finally, volume correction for which the Hoshen-Kopelman (HK) clustering algorithm [START_REF] Hoshen | Percolation and cluster distribution. I. Cluster multiple labeling technique and critical concentration algorithm[END_REF] is used to remove non-connected solid voxel sets. The object has hence only two phases: a 26-connected solid phase and a 6-connected pore phase. The topology of a porous medium can be characterized using the Betti numbers. Considering a 3D space, there are 3 distinct Betti numbers that completely define the topology of an object. β 0 is the number of connected elements of the solid phase Ω. β 1 is the number of loops and closed paths of Ω. β 2 is the number of internal surfaces of Ω. The connectivity of an object is usually evaluated using the Euler-Poincare Characteristic, N3, which is linked to the Betti numbers by

N 3 = β 0 -β 1 + β 2 .
To overcome the limitations of curve and surface skeletons for modeling objects of mixed rod/plate shapes, we developed a new skeleton called Hybrid which is based on an original combination of three techniques: element classification as rods or plates, surface thinning and curve thinning. First, a method of shape classification is used to decompose the object into two sets of plates and rods [START_REF] Bonnassie | A new method for analyzing local shape in three-dimensional images based on medial axis transformation[END_REF]. This classification is performed at the voxel level as shown in where TV is the total volume in mm 3 , ρ is the voxel resolution in mm 3 , Card(y) is the number of elements in the set y, S H is the set of voxels of the hybrid skeleton, N (m) represents the neighborhood of voxel m, and S H P and S H are the rod and plate parts of S H , respectively. The results of table 1 represent the values of the previously described parameters extracted from the hybrid skeleton. The SV /T V parameter reflects the presence of plates in the porous samples which increases its density. The number of line-ends is an indicator of the number of broken branches. The hybrid skeleton distinguishes plates and rods, so as to provide two efficient and complementary values.

Feature

Table 1 :

 1 Conn.D (mm -1 ) SV/TV No.N/TV (mm -1 ) Le.N/TV (mm -1 ) Values of different parameters characterizing the porous objects of figure 1.

	Rd/Pl

figures 1(b). After the classification step, surface [START_REF] Mokherjee | On connectivity issues of ESPTA[END_REF] and curve [START_REF] Morgenthaler | Three dimensional simple points: serial erosion, parallel thinning and skeletonization[END_REF], thinning techniques are applied respectively to the plate and rod subsets. The joint use of these thinning techniques generates a skeleton composed of both 2D surfaces at plate shapes and 1D paths at rod shapes, as shown in figures 1(c).

Results and Discussion

Let Ω and Ω c be the set of solid and pore voxels, respectively. Interesting features characterizing the morphology and the topology of each sample can be extracted directly from the skeleton. These features include: skeleton Connectivity Density (Conn.D), Skeleton density (SV /T V ), normalized Number of Line-end voxels (Le.N ), normalized Number of Node voxels (N o.N ) and skeleton Rods to Plates ratio (Rd/P l).


Benedikt Bollig

Paul Gastin

Benjamin Monmege

Marc Zeitoun

Pebble weighted automata and transitive closure logics ⋆

We introduce new classes of weighted automata on words. Equipped with pebbles and a two-way mechanism, they go beyond the class of recognizable formal power series, but capture a weighted version of first-order logic with bounded transitive closure. In contrast to previous work, this logic allows for unrestricted use of universal quantification. Our main result states that pebble weighted automata, nested weighted automata, and this weighted logic are expressively equivalent. We also give new logical characterizations of the recognizable series.

Introduction

Connections between logical and state-based formalisms have always been a fascinating research area in theoretical computer science, which produced some fundamental theorems. The line of classical results started with the equivalence of MSO logic and finite automata [START_REF] Büchi | Weak second-order arithmetic and finite automata[END_REF][START_REF] Elgot | Decision problems of finite automata design and related arithmetics[END_REF][START_REF] Trakhtenbrot | Finite automata and logic of monadic predicates[END_REF]. Some extensions of finite automata are of quantitative nature and include timed automata, probabilistic systems, and transducers, which all come with more or less natural, specialized logical characterizations. A generic concept of adding weights to qualitative systems is provided by the theory of weighted automata [START_REF]Handbook of Weighted Automata[END_REF], first introduced by Schützenberger [START_REF] Schützenberger | On the definition of a family of automata[END_REF]. The output of a weighted automaton running on a word is no longer a Boolean value discriminating between accepted and rejected behaviors. A word is rather mapped to a weight from a semiring, summing over all possible run weights, each calculated as the product of its transition outcomes. Indeed, probabilistic automata and word transducers appear as instances of that framework (see [7, Part IV]).

A logical characterization of weighted automata, however, was established only recently [START_REF] Droste | Weighted automata and weighted logics[END_REF], in terms of a (restricted) weighted MSO logic capturing the recognizable formal power series (i.e., the behaviors of finite weighted automata). The key idea is to interpret existential and universal quantification as sum and product from a semiring. To make this definition work, however, one has to restrict the universal first-order quantification, which, otherwise, appears to be too powerful and goes beyond the class of recognizable series. In this paper, we follow a different approach. Instead of restricting the logic, we define an extended automata model that naturally reflects it. Indeed, it turns out that universal quantification is essentially captured by a pebble (two-way) mechanism in the automata-theoretic counterpart. Inspired by the theory of two-way and pebble automata on words and trees [START_REF] Globerman | Complexity results for two-way and multi-pebble automata and their logics[END_REF][START_REF] Engelfriet | Tree-walking pebble automata[END_REF][START_REF] Bojańczyk | Tree-walking automata[END_REF], we actually define weighted generalizations that preserve their natural connections with logic.

More precisely, we introduce pebble weighted automata on words and establish expressive equivalence to weighted first-order logic with bounded transitive closure and unrestricted use of quantification, extending the classical Boolean case for words [START_REF] Engelfriet | Automata with nested pebbles capture firstorder logic with transitive closure[END_REF]. Our equivalence proof makes a detour via another natural concept, named nested weighted automata, which resembles the nested treewalking automata of [START_REF] Cate | Transitive closure logic, nested tree walking automata, and XPath[END_REF]. The transitive closure logic also yields alternative characterizations of the recognizable formal power series.

Proofs omitted due to lack of space are available in [START_REF] Bollig | Pebble weighted automata and transitive closure logics[END_REF].

Notation and background

In this section we set up the notation and we recall some basic results on weighted automata and weighted logics. We refer the reader to [START_REF] Droste | Weighted automata and weighted logics[END_REF][START_REF]Handbook of Weighted Automata[END_REF] for details. Throughout the paper, Σ denotes a finite alphabet and Σ + is the free semigroup over Σ, i.e., the set of nonempty words. The length of u

∈ Σ + is de- noted |u|. If |u| = n ≥ 1, we usually write u = u 1 • • • u n with u i ∈ Σ and we let Pos(u) = {1, . . . , n}. For 1 ≤ i ≤ j ≤ n, we denote by u[i..j] the factor u i u i+1 • • • u j of u. Finally, we let Σ ≤k = 1≤i≤k Σ i .

Formal power series.

A semiring is a structure K = (K, +, •, 0, 1) where (K, +, 0) is a commutative monoid, (K, •, 1) is a monoid, • distributes over +, and 0 is absorbing for •. We say that K is commutative if so is (K, •, 1). We shall refer in the examples to the usual Boolean semiring B = ({0, 1}, ∨, ∧, 0, 1) and to the semiring (N, +, • , 0, 1) of natural numbers, denoted N. A formal power series (or series, for short) is a mapping f : Σ + → K. The set of series is denoted K Σ + . We denote again by + and • the pointwise addition and multiplication (called the Hadamard product) on K Σ + , and by 0 and 1 the constant series with values 0 and 1, respectively. Then (K Σ + , +, •, 0, 1) is itself a semiring.

Weighted automata. All automata we consider are finite. A weighted automaton (wA) over K = (K, +, •, 0, 1) and Σ is a tuple A = (Q, µ, λ, γ), where Q is the set of states, µ : Σ → K Q×Q is the transition weight function and λ, γ : Q → K are weight functions for entering and leaving a state. The function µ gives, for each a ∈ Σ and p, q ∈ Q, the weight µ(a) p,q of the transition p a -→ q. It extends uniquely to a homomorphism µ :

Σ + → K Q×Q . Viewing µ as a mapping µ : Q × Σ + × Q → K, we sometimes write µ(p, u, q) instead of µ(u) p,q . A run on a word u = u 1 • • • u n is a sequence of transitions ρ = p 0 u1 -→ p 1 u2 -→ • • • un --→ p n . The weight of the run ρ is weight(ρ) def = λ(p 0 ) • [ n i=1 µ(p i-1 , u i , p i )] • γ(p n ),
and the weight A (u) of u is the sum of all weights of runs on u, which can also be computed as A (u) = λ • µ(u) • γ, viewing λ, µ(u), γ as matrices of dimension 1 × |Q|, |Q| × |Q| and |Q| × 1, respectively. We call A ∈ K Σ + the behavior, or semantics of A. A series f ∈ K Σ + is recognizable if it is the behavior of some wA. We let K rec Σ + be the collection of all recognizable series. Example 1. Consider (N, +, • , 0, 1) and let A be the automaton with a single state, µ(a) = 2 for all a ∈ Σ, and λ = γ = 1. Then, A (u) = 2 |u| for all u ∈ Σ + .

It is well-known that K rec Σ + is stable under + and, if K is commutative, also under •, making (K rec Σ + , +, •, 0, 1) a subsemiring of (K Σ + , +, •, 0, 1).

Weighted logics. We fix infinite supplies Var = {x, y, z, t, . . .} of first-order variables, and VAR = {X, Y, . . .} of second-order variables. The class of weighted monadic second-order formulas over K and Σ, denoted MSO(K, Σ) (shortly MSO), is given by the following grammar, with k ∈ K, a ∈ Σ, x, y ∈ Var and X ∈ VAR:

ϕ ::= k | P a (x) | x ≤ y | x ∈ X | ¬ϕ | ϕ ∨ ϕ | ϕ ∧ ϕ | ∃xϕ | ∀xϕ | ∃Xϕ | ∀Xϕ .
For ϕ ∈ MSO(K, Σ), let Free(ϕ) denote the set of free variables of ϕ. If Free(ϕ) = ∅, then ϕ is called a sentence. For a finite set V ⊆ Var ∪ VAR and a word u ∈ Σ + , a (V, u)-assignment is a function σ that maps a first-order variable in V to an element of Pos(u) and a second-order variable in V to a subset of Pos(u). For x ∈ Var and i ∈ Pos(u), σ[x → i] denotes the (V ∪ {x}, u)assignment that maps x to i and, otherwise, coincides with σ. For X ∈ VAR and I ⊆ Pos(u), the (V ∪ {X}, u)-assignment σ[X → I] is defined similarly.

A pair (u, σ), where σ is a (V, u)-assignment, can be encoded as a word over the extended alphabet

Σ V def = Σ × {0, 1} V . We write a word (u 1 , σ 1 ) • • • (u n , σ n ) ∈ Σ + V as (u, σ) where u = u 1 • • • u n and σ = σ 1 • • • σ n .
We call (u, σ) valid if, for each first-order variable x ∈ V, the x-row of σ contains exactly one 1. If (u, σ) is valid, then σ can be considered as the (V, u)-assignment that maps a first-order variable x ∈ V to the unique position carrying 1 in the x-row, and a second-order variable X ∈ V to the set of positions carrying 1 in the X-row.

Fix a finite set V of variables such that Free(ϕ) ⊆ V. The semantics ϕ V ∈ K Σ + V of ϕ wrt. V is given as follows: if (u, σ) is not valid, we set ϕ V (u, σ) = 0, otherwise ϕ V is given by Figure 1. Hereby, the product follows the natural order on Pos(u) and some fixed order on the power set of Pos(u). We simply write ϕ for ϕ Free(ϕ) and say that ϕ is recognizable if so is ϕ . We note K MSO Σ + the class of all series definable by a sentence of MSO(K, Σ).

Example 2. For K = B, recognizable and MSO(K, Σ)-definable languages coincide. In contrast, for K = (N, +, • , 0, 1), the very definition yields ∀x∀y 2 (u) = 2 |u| 2 , which is not recognizable [START_REF] Droste | Weighted automata and weighted logics[END_REF]. Indeed, the function computed by a wA A satisfies A (u) = 2 O(|u|) . Also observe that the behavior of the automaton of Example 1 is ∀y 2 . Therefore, recognizable series are not stable under universal first-order quantification.

Let bMSO(K, Σ) be the syntactic Boolean fragment of MSO(K, Σ) given by

ϕ ::= 0 | 1 | P a (x) | x ≤ y | x ∈ X | ¬ϕ | ϕ ∧ ϕ | ∀xϕ | ∀Xϕ , ϕ1 ∨ ϕ2 V (u, σ) = ϕ1 V (u, σ) + ϕ2 V (u, σ) k V (u, σ) = k, for k ∈ K ϕ1 ∧ ϕ2 V (u, σ) = ϕ1 V (u, σ) • ϕ2 V (u, σ) Pa(x) V (u, σ) =  1 if u σ(x) = a 0 otherwise ∃xϕ V (u, σ) = X i∈Pos(u) ϕ V∪{x} (u, σ[x → i]) x ∈ X V (u, σ) =  1 if σ(x) ∈ σ(X) 0 otherwise ∃Xϕ V (u, σ) = X I⊆Pos(u) ϕ V∪{X} (u, σ[X → I]) x ≤ y V (u, σ) =  1 if σ(x) ≤ σ(y) 0 otherwise ∀xϕ V (u, σ) = Y i∈Pos(u) ϕ V∪{x} (u, σ[x → i]) ¬ϕ V (u, σ) =  1 if ϕ V (u, σ) = 0 0 otherwise ∀Xϕ V (u, σ) = Y I⊆Pos(u) ϕ V∪{X} (u, σ[X → I])
+ → ψ def = ¬ϕ∨(ϕ∧ψ) so that, if ϕ is a Boolean formula (i.e., ϕ (Σ + ) ⊆ {0, 1}), ϕ + → ψ (u, σ) = ψ (u, σ) if ϕ (u, σ) = 1, and ϕ + → ψ (u, σ) = 1 if ϕ (u, σ) = 0.
A common fragment of MSO(K, Σ) is the weighted first-order logic FO(K, Σ), where no second-order quantifier appears (note that second order variables may still appear free):

ϕ ::= k | P a (x) | x ≤ y | x ∈ X | ¬ϕ | ϕ ∨ ϕ | ϕ ∧ ϕ | ∃xϕ | ∀xϕ.
We define similarly bFO(K, Σ) as the fragment of bMSO(K, Σ) with no second-order quantifiers. We also let bFO+mod be the fragment of bMSO consisting of bFO augmented with modulo constraints x ≡ ℓ m for constants 1 ≤ m ≤ ℓ (since the positions of words start with 1, it is more convenient to compute modulo as a value between 1 and ℓ). The semantics is given by x ≡ ℓ m (u, σ) = 1 if σ(x) ≡ m mod ℓ and 0 otherwise: it can be defined in bMSO by

x ≡ ℓ m def = ∀X (x ∈ X) ∧ ∀y(y ∈ X ∧ y > ℓ) + → y -ℓ ∈ X + → m ∈ X .
For L ⊆ bMSO closed under ∨, ∧ and ¬, an L-step formula is a formula obtained from the grammar ϕ ::

= k | α | ¬ϕ | ϕ ∨ ϕ | ϕ ∧ ϕ, with k ∈ K and α ∈ L.
In particular, quantifications are only allowed in formulas α ∈ L. The following lemma shows in particular that an L-step formula assumes a finite number of values, each of which corresponds to an L-definable language. Lemma 3. For every L-step formula ϕ, one can construct an equivalent formula ψ = i (ϕ i ∧ k i ) with ϕ i ∈ L and k i ∈ K, with the same set of free variables.

From now on, we freely use Lemma 3, using the special form it provides for L-step formulas. All bMSO-step formulas are clearly recognizable. By [START_REF] Droste | Weighted automata and weighted logics[END_REF], ∀x ϕ is recognizable for any bMSO-step formula ϕ. The fragment RMSO(K, Σ) of MSO is defined by restricting universal second-order quantification to bMSO formulas and universal first-order quantification to bMSO-step formulas.

Theorem 4 ([6]

). A series is recognizable iff it is definable in RMSO(K, Σ).

Transitive closure logic and weighted automata

To ease notation, we write ϕ(x, y) (u, i, j) instead of ϕ(x, y) (u, [x → i, y → j]). We allow constants, modulo constraints and comparisons, like e.g. x ≤ y + 2. We use first and last as abbreviations for the first and last positions of a word. All of these shortcuts can be replaced by suitable bFO-formulas, except "x ≡ ℓ m" with 1 ≤ m ≤ ℓ, which is bMSO-definable.

Bounded transitive closure. For a formula ϕ(x, y) with at least two free variables x and y, and an integer N > 0, we let ϕ 1,N (x, y) def = (x ≤ y ≤ x + N ) ∧ ϕ(x, y) and for n ≥ 2, we define the formula ϕ n,N (x, y) as

∃z 0 • • • ∃z n x = z 0 ∧ y = z n ∧ 1≤ℓ≤n (z ℓ-1 < z ℓ ≤ z ℓ-1 + N ) ∧ ϕ(z ℓ-1 , z ℓ ) . (1)
We define for each N > 0 the N -TC < xy operator by N -TC < xy ϕ = n≥1 ϕ n,N . This infinite disjunction is well-defined: ϕ n,N (x, y) (u, σ) = 0 if n ≥ max(2, |u|), i.e., on each pair (u, σ), only finitely many disjuncts assume a nonzero value. Intuitively, the N -TC < xy operator generalizes the forward transitive closure operator from the Boolean case, but limiting it to intermediate forward steps of length ≤ N . The fragment FO+BTC < (K, Σ) is then defined by the grammar

ϕ ::= k | P a (x) | x ≤ y | ¬ϕ | ϕ ∨ ϕ | ϕ ∧ ϕ | ∃xϕ | ∀xϕ | N -TC < xy ϕ ,
with N ≥ 1 and the restriction that one can apply negation only over bFOformulas. We denote by K FO+BTC < Σ + the class of all FO+BTC < (K, Σ)definable series.

Example 5. Let ϕ(x, y)

def = (y = x+1)∧∀z 2∧(x = 1 + → ∀z 2) over K = N. Let u = u 1 • • • u n . For any N ≥ 1, we have N -TC < xy ϕ (u, first, last) = n-1 i=1 ϕ (u, i, i+1) due to the constraint y = x+1 in ϕ. Now, ϕ (u, 1, 2) = 2 2|u| and ϕ (u, i, i+1) = 2 |u| if i > 1, so N -TC < xy ϕ (u, first, last) = 2 |u| 2 .
This example shows in particular that the class of recognizable series is not closed under BTC < . Example 6. It is well-known that modulo can be expressed in bFO+BTC < by

x ≡ ℓ m def = (x = m) ∨ [ℓ-TC < yz (z = y + ℓ)](m, x) .
We now consider syntactical restrictions of FO+BTC < , inspired by normal form formulas of [START_REF] Neven | On the power of tree-walking automata[END_REF] where only one "external" transitive closure is allowed.

For L ⊆ bMSO, BTC < step (L) consists of formulas of the form N -TC < xy ϕ, where ϕ(x, y) is an L-step formula with two free variables x, y. We say that f ∈ K Σ + is BTC < step (L)-definable if there exists an L-step formula ϕ(x, y) such that, for all u ∈ Σ + , f (u) = N -TC < xy ϕ (u, first, last). For L ⊆ bMSO, let ∃∀ step (L) consists of all MSO-formulas of the form ∃X∀x ϕ(x, X) with ϕ an L-step formula: this defines a fragment of the logic RMSO(K, Σ) introduced in [START_REF] Droste | Weighted automata and weighted logics[END_REF]. The following result characterizes the expressive power of weighted automata. Theorem 7. Let K be a (possibly noncommutative) semiring and f ∈ K Σ + . The following assertions are equivalent over K and Σ:

(1) f is recognizable. (2) f is BTC < step (bFO+mod)-definable. (3) f is BTC < step (bMSO)-definable.
(4) f is ∃∀ step (bFO)-definable.

(5) f is ∃∀ step (bMSO)-definable.

Proof. Fix a weighted automaton A = (Q, µ, λ, γ) with Q = {1, . . . , n}. For d ≥ 1 and p, q ∈ Q, we use a formula ψ d p,q (x, y) to compute the weight of the factor of length d located between positions x and y, when A goes from p to q:

ψ d p,q (x, y) def = (y = x + d -1) ∧ v=v1•••v d µ(v) p,q ∧ 1≤i≤d P vi (x + i -1) .
We construct a formula ϕ(x, y) of bFO+mod allowing to define the semantics of A using a BTC < : A (u) = 2n-TC < xy ϕ (u, first, last). The idea, inspired by [START_REF] Thomas | Classifying regular events in symbolic logic[END_REF], consists in making the transitive closure pick positions z ℓ = ℓn + q ℓ , with 1 ≤ q ℓ ≤ n, for successive values of ℓ, to encode runs of A going through state q ℓ just before reading the letter at position ℓn + 1. To make this still work for ℓ = 0, one can assume wlog. that λ(1) = 1 and λ(q) = 0 for q = 1, i.e., the only initial state yielding a nonzero value is q 0 = 1. Consider slices [ℓn + 1, (ℓ + 1)n] of positions in the word where we evaluate the formula (the last slice might be incomplete). Each position y is located in exactly one such slice. We write y = ℓn + 1 for the first position of that slice, as well as [y] def = y + 1y ∈ Q for the corresponding "offset". Notice that, for q ∈ Q, [y] = q can be expressed in bFO+mod simply by y ≡ n q. Hence, we will freely use [y] as well as y = y + 1 -[y] as macros in formulas. Our BTC < -formula picks positions x and y marked • in Figure 2, and computes the weight of the factor of length n between the positions x and y -1, assuming states [x] and [y] just before and after these positions. The formula ϕ distinguishes the cases where x is far or near from the last position:

ϕ(x, y) = ( x + 2n ≤ last) ∧ ( y = x + n) ∧ ψ n [x],[y] ( x , y -1) ∨ ( x + 2n > last) ∧ (y = last) ∧ q∈Q ψ y-x +1 [x],q ( x , y) ∧ γ(q) . (ℓ -1)n + 1 ℓn ℓn + 1 (ℓ + 1)n x y • x • y [x]
[y]

ψ n [x],[y]
Fig. 2. Positions picked by the BTC < -formula

Weighted nested automata

Example 2 shows that weighted automata lack closure properties to capture FO+BTC < . We introduce a notion of nested automata making up for this gap. For r ≥ 0, the class r-nwA(Σ) (r-nwA if Σ is understood) of r-nested weighted automata over Σ (and K) consists of all tuples (Q, µ, λ, γ) where Q is the set of states, λ, γ : Q → K, and µ :

Q × Σ × Q → (r -1)-nwA(Σ × {0, 1}).
Here, we agree that (-1)-nwA = K. In particular, a 0-nwA(Σ) is a weighted automaton over Σ. Intuitively, the weight of a transition is computed by an automaton of the preceding level running on the whole word, where the additional {0, 1} component marks the letter of the transition whose weight is to be computed.

Let us formally define the behavior

A ∈ K Σ + of A = (Q, µ, λ, γ) ∈ r-nwA(Σ). If r = 0, then A is the behavior of A considered as wA over Σ. For r ≥ 1, the weight of a run ρ = q 0 u1 -→ q 1 u2 -→ • • • un --→ q n of A on u = u 1 • • • u n is weight(ρ) def = λ(q 0 ) • n i=1 µ(q i-1 , u i , q i ) (u, i) • γ(q n ) , where (u, i) ∈ (Σ × {0, 1}) + is the word v = v 1 • • • v n with v i = (u i , 1
) and v j = (u j , 0) if j = i. As usual, A (u) is the sum of the weights of all runs of A on u. Note that, unlike the nested automata of [START_REF] Cate | Transitive closure logic, nested tree walking automata, and XPath[END_REF], the values given by lower automata do not explicitly influence the possible transitions.

A series f ∈ K Σ + is r-nwA-recognizable if f = A for some r-nwA A. It is nwA-recognizable if it is r-nwA-recognizable for some r. We let K r-nwA Σ + (resp., K nwA Σ + ) be the class of r-nwA-recognizable (resp., nwA-recognizable) series over K and Σ.

Example 8. A 1-nwA recognizing the series u → 2 |u| 2 over N is A = ({p}, µ, 1, 1)
where, for every a ∈ Σ, µ(p, a, p) is the weighted automaton of Example 1.

We can generalize the proof of (1) ⇒ (2) in Theorem 7 in order to get the following result. The converse will be obtained in Section 5.

Proposition 9. Every nwA-recognizable series is FO+BTC < -definable.

Pebble weighted automata

We now consider pebble weighted automata (pwA). A pwA has a read-only tape. At each step, it can move its head one position to the left or to the right (within the boundaries of the input tape), or either drop or lift a pebble at the current head position. Applicable transitions and weights depend on the current letter, current state, and the pebbles carried by the current position. Pebbles are handled using a stack policy: if the automaton has r pebbles and pebbles ℓ, . . . , r have already been dropped, it can either lift pebble ℓ (if ℓ ≤ r), drop pebble ℓ -1 (if ℓ ≥ 2), or move. As these automata can go in either direction, we add two fresh symbols ⊲ and ⊳ to mark the beginning and the end of an input word. Let Σ = Σ ⊎ {⊲, ⊳}. To compute the value of w = w 1 • • • w n ∈ Σ + , a pwA will work on a tape holding w = ⊲w⊳. For convenience, we number the letters of w from 0, setting w 0 = ⊲, w n+1 = ⊳, and

w i = w i for 1 ≤ i ≤ n.
Let r ≥ 0. Formally, an r-pebble weighted automaton (r-pwA) over K and Σ is a pair A = (Q, µ) where Q is a finite set of states and µ : Q× Σ ×2 r ×D×Q → K is the transition weight function, with D = {←, →, drop, lift}.

A configuration of a r-pwA A on a word w ∈ Σ + of length n is a triple (p, i, ζ) ∈ Q×{0, . . . , n+2}×{1, . . . , n} ≤r . The word w itself will be understood. Informally, p denotes the current state of A and i is the head position in w, i.e positions 0 and n + 1 point to ⊲ and ⊳, respectively, position 1 ≤ i ≤ n points to w i ∈ Σ, and position n + 2 is outside w.

Finally, ζ = ζ ℓ • • • ζ r with 1 ≤ ℓ ≤ r + 1
encodes the locations of pebbles ℓ, . . . , r (ζ m ∈ {1, . . . , n} is the position of pebble m) while pebbles 1, . . . , ℓ-1 are currently not on the tape. For i ∈ {0, . . . , n+1},

we set ζ -1 (i) = m ∈ {ℓ, . . . , r} | ζ m = i (viewing ζ as a partial function from {1, . . . , r} to {0, . . . , n + 1}). Note that ζ -1 (0) = ζ -1 (n + 1) = ∅.
There is a step of weight k from configuration (p, i, ζ) to configuration (q, j, η)

if i ≤ n + 1, k = µ(p, w i , ζ -1 (i), d, q), and      j = i -1 if d = ← j = i + 1 if d = → j = i otherwise and      η = iζ if d = drop ζ = iη if d = lift η = ζ otherwise.
A run ρ of A is a sequence of steps from a configuration (p, 0, ε) to a configuration (q, n + 2, ε) (at the end, no pebble is left on the tape). We denote by weight(ρ) the product of weights of the steps of run ρ (from left to right, but we will mainly work with a commutative semiring in this section). The run ρ is simple if whenever two configurations α and β appear in ρ, we have α = β. The series A ∈ K Σ + is defined by A (w) = ρ simple run on w weight(ρ). We denote by K r-pwA Σ + the collection of formal power series definable by a r-pwA, and we let K pwA Σ + = r≥0 K r-pwA Σ + . Note that a 0-pwA is in fact a 2-way weighted automaton. It follows from Theorem 11 that 2-way wA have the same expressive power as classical (1-way) wA.

Example 10. Let us sketch a 1-pwA A recognizing the series u → 2 |u| 2 over N. The idea is that A drops its pebble successively on every position of the input word. Transitions for reallocating the pebble have weight 1. When a pebble is dropped, A scans the whole word from left to right where every transition has weight 2. As this scan happens |u| times, we obtain A (u) = 2 |u| 2 . Theorem 11. For every commutative semiring K, and every r ≥ 0, we have

(1) K r-pwA Σ + ⊆ K r-nwA Σ + (2) K FO+BTC < Σ + = K pwA Σ + = K nwA Σ + .
(3) K FO+BTC < Σ + ⊆ K pwA Σ + holds even for noncommutative semirings.

Proof. [START_REF] Berstel | Recognizable formal power series on trees[END_REF] We provide a translation of a generalized version of r-pwA to r-nwA. That generalized notion equips an r-pwA A = (P, µ) with an equivalence relation ∼ ⊆ P × P , which is canonically extended to configurations of A: we write (p, i, u) fined by replacing equality of configurations with ∼ in the definition of simple run. To stress this fact, we henceforth say that a run is ∼simple. So let r ≥ 0, A = (P, µ) be an r-pwA over K and Σ, and ∼ ⊆ P × P be an equivalence relation. We assume wlog. that all runs in which a pebble is dropped and immediately lifted have weight 0. We build an r-nwA

∼ (p ′ , i ′ , u ′ ) if p ∼ p ′ , i = i ′ ,
A ∼ = (Q, ν, λ, γ) over Σ such that A ∼ = A ∼ .
The construction of A ∼ proceeds inductively on the number of pebbles r. It involves two alternating transformations, as illustrated in Figure 3. The left-hand side depicts a ∼-simple run of A on some word w with factor ab. To simulate such a run, A ∼ scans w from left to right and guesses, at each position i, the sequence of those states and directions that are encountered at i while pebble r has not, or just, been dropped. The state of A ∼ taken before reading the a at position i is q = p 0 → p 3 ← p 4 drop p 5 p5 lift p 6 ← p 7 drop p 8 p8 lift p 9 → (which is enriched by the input letter a, as will be explained below). As the micro-states p 0 , p 3 , . . . form a segment of a ∼-simple run, p 0 , p 3 , p 4 , p 6 , p 7 , p 9 are pairwise distinct (wrt. ∼) and so are p 5 , p5 , p 8 , p8 . Segments when pebble r is dropped on position i are defered to a (r -1)-nwA B q , which is called at position i and computes the run segments from p 5 to p5 and from p 8 to p8 that both start in i. To this aim, B q works on an extension of w where position i is marked, indicating that pebble r is considered to be at i.

We define the r-nwA A ∼ . Let Π = (P {→, ←} ∪ P {drop}P P {lift}) * P {→}. Sequences from Π keep track of states and directions that are taken at one given position. As aforementioned, they must meet the requirements of ∼-simple runs so that only some of them can be considered as states. Formally, given the dashed part of the run of A, we use phase numbers, so that, during the simulation two different states (p, j) and (p, j ′ ) of A q may correspond to the same original state p of A. Now, only simple runs are considered to compute A . Therefore, for the simulation to be faithful, we want to rule out runs of A q containing two configurations which only differ by the phase number, that is, containing two ∼ q -equivalent configurations, which is done by keeping only ∼ q -simple runs. A run of A q is illustrated in Figure 3(b) (with N = 2). Note that, if N = 0, then A q simply scans the word from left to right, outputting weight 1. Finally, we sketch the proof of (3).

We proceed by induction on the structure of the formula and suppose that a valuation of free variables is given in terms of pebbles that are already placed on the word and cannot be lifted. Disjunction, conjunction, and first-order quantifications are easy to simulate. To evaluate [N -TC < xy ϕ](z, t) for some formula ϕ(x, y), we either evaluate ϕ(x, y), with pebbles being placed on z and t such that z ≤ t ≤ z + N , or choose non-deterministically (and with weight 1) positions z = z 0 < z 1 < • • • < z n-1 < z n = t with n ≥ 2, using two additional pebbles, 2 and 1. We drop pebble 2 on position z 0 and pebble 1 on some guessed position z 1 with z 0 < z 1 ≤ min(t, z 0 + N ). We then run the subroutine to evaluate ϕ(z 0 , z 1 ). Next we move to position z 1 and lift pebble 1. We move left to position z 0 remembering the distance z 1z 0 . We lift pebble 2 and move right to z 1 using the stored distance z 1z 0 . We drop pebble 2 on z 1 and iterate this procedure until t is reached.

⊓ ⊔

Conclusion and perspectives

We have introduced pebble weighted automata and characterized their expressive power in terms of first-order logic with a bounded transitive closure operator. It follows that satisfiability is decidable over commutative positive semirings. Here, a sentence ϕ ∈ FO+BTC < is said satisfiable if there is a word w ∈ Σ + such that ϕ (w) = 0. From Theorem 11, satisfiability reduces to non-emptiness of the support of a series recognized by a pwA. For positive semirings, the latter problem, in turn, can be reduced to the decidable emptiness problem for classical pebble automata over the Boolean semiring. We leave it as an open problem to determine for which semirings the satisfiability problem is decidable.

Unbounded transitive closure. We do not know if allowing unbounded steps in the transitive closure leads beyond the power of (weak) pebble automata. We already know that allowing unbounded steps is harmless for bMSO-step formulas. It is also easy to show that such an unbounded transitive closure can be captured with strong pebble automata, i.e., that can lift the last dropped pebble even when not scanning its position. Therefore, we aim at studying the expressive power of strong pebble automata and unbounded transitive closure.

Tree-walking automata. Our results are not only of theoretical interest. They also lay the basis for quantitative extensions of database query languages such as XPath, and may provide tracks to evaluate quantitative aspects of the structure of XML documents. The framework of weighted tree (walking) automata [START_REF] Berstel | Recognizable formal power series on trees[END_REF][START_REF] Fülöp | Weighted Tree-Walking Automata[END_REF] is natural for answering questions such as "How many nodes are selected by a request?", or "How difficult is it to answer a query?". The navigational mechanism of pebble tree-walking automata [START_REF] Bojańczyk | Expressive power of pebble automata[END_REF][START_REF] Samuelides | Complexity of pebble tree-walking automata[END_REF][START_REF] Bojańczyk | Tree-walking automata[END_REF] is also well-suited in this context. For these reasons, we would like to adapt our results to tree languages.

We thank the anonymous referees for their valuable comments.

Fig. 1 .

 1 Fig. 1. Semantics of weighted MSO

Fig. 3 .

 3 Fig. 3. (a) Runs of r-pwA A and r-nwA A ∼; (b) run of (r -1)-pwA Aq

⋆ Supported by fp7 Quasimodo, anr-06-seti-003 dots, arcus Île de France-Inde.

π ∈ Π, we define projections pr 1 (π) ∈ P + and pr 2 (π) ∈ (P P ) * inductively by pr 1 (ε) = pr 2 (ε) = ε and pr 1 (p → π) = pr 1 (p ← π) = p pr 1 (π) pr 1 (p drop p 1 p1 lift π) = p pr 1 (π) pr 2 (p → π) = pr 2 (p ← π) = pr 2 (π) pr 2 (p drop p 1 p1 lift π) = p 1 p1 pr 2 (π).

Let Π ∼ denote the set of sequences π ∈ Π such that pr 1 (π) consists of pairwise distinct states wrt. ∼ and so does pr 2 (π) (there might be states that occur in both pr 1 (π) and pr 2 (π)). With this, we set Q = ( Σ ⊎ {✷}) × Π ∼ . The letter a ∈ Σ ⊎ {✷} of a state (a, π) ∈ Q will denote the symbol that is to be read next. Symbol ✷ means that there is no letter left so that the automaton is beyond the scope of ⊲w⊳ when w is the input word.

Next, we explain how the weight of the run segments of A with lifted pebble r is computed in A ∼ . Two neighboring states of A ∼ need to match each other, which can be checked locally by means of transitions. To determine a corresponding weight, we first count weights of those transitions that move from the current position i to i + 1 or from i + 1 to i. This is the reason why a state of A ∼ also maintains the letter that is to be read next. In Figure 3(a), the 7 micro-transitions that we count in the step from q to q ′ are highlighted in gray. Assuming q = (a, π 0 ) and q ′ = (b, π 1 ), we obtain a value weight a,b (π 0 | π 1 ) as the product µ(p 0 , a, ∅, →,

We are now prepared to define the components ν, λ, γ of A ∼ . For q 0 = (a 0 , π 0 ) and q 1 = (a 1 , π 1 ) states in Q, we set

Here, B q is the constant 1 if r = 0. Otherwise, B q = A q ∼q is an (r -1)-nwA over ∆ = Σ × {0, 1} obtained inductively from the (r -1)-pwA A q which is defined below together with its equivalence relation ∼ q . Notice that k • B q is obtained from B q by multiplying its input weights by k.

Let us define the (r -1)-pwA A q = (P ′ , µ ′ ) over ∆ as well as ∼ q ⊆ P ′ × P ′ . Suppose q = (a, π) and pr 2 (π) = p 1 p1 • • • p N pN . The behavior of A q is split into N + 2 phases. In phase 0, it scans the input word from left to right until it finds the (unique) letter of the form (a, 1) with a ∈ Σ. At that position, call it i, A q enters state (p 1 , 1) (the second component indicating the current phase). All these transitions are performed with weight 1. Then, A q starts simulating A, considering pebble r at position i. Back at position i in state (p 1 , 1), weight-1 transitions will allow A q to enter the next phase, starting in (p 2 , 2) and again considering pebble r at position i. The simulation of A ends when (p N , N ) is reached in position i. In the final phase (N +1) weight-1 transitions guides A q to the end of the tape where it stops. The relation ∼ q will consider states (p, j) and (p ′ , j ′ ) equivalent iff p ∼ p ′ , i.e., it ignores the phase number. This explains the purpose of the equivalence relation: in order for the automaton A q to simulate


Nicolas Jay
email: jay@loria.fr

François Kohler

Amedeo Napoli

Analysis of Social Communities with Iceberg and Stability-Based Concept Lattices

In this paper, we presents a research work based on formal concept analysis and interest measures associated with formal concepts. This work focuses on the ability of concept lattices to discover and represent special groups of individuals, called social communities. Concept lattices are very useful for the task of knowledge discovery in databases, but they are hard to analyze when their size become too large. We rely on concept stability and support measures to reduce the size of large concept lattices. We propose an example from real medical use cases and we discuss the meaning and the interest of concept stability for extracting and explaining social communities within a healthcare network.

Introduction

Knowledge Discovery in Databases (KDD) is an iterative and interactive process for identifying valid, novel, and potentially useful patterns in data [START_REF] Fayyad | The kdd process for extracting useful knowledge from volumes of data[END_REF]. It is usually divided into three main steps: data preparation, data mining, and interpretation of the extracted units. Data mining is often considered as the central step in the KDD process. However, interpretation of data-mining results is also an important step within the KDD process. Indeed, one of the the success keys in KDD practice relies on the ability of easily producing units understandable as knowledge units. One way of achieving such a goal is to provide an adapted organization and representation of the extracted units, especially when the KDD system has to be used by novice users.

In parallel, Formal Concept Analysis (FCA) is a theory of data analysis introduced in [START_REF] Wille | Restructuring lattice theory: an approach based on hierarchies of concepts[END_REF], that is tightly connected with KDD [START_REF] Wille | Why can concept lattices support knowledge discovery in databases?[END_REF][START_REF] Valtchev | Formal concept analysis for knowledge discovery and data mining: The new challenges[END_REF], particularly regarding the search of frequent itemsets and the extraction of association rules [START_REF] Agrawal | Mining association rules between sets of items in large databases[END_REF]. Many algorithms relying on FCA central property of closure have been proposed to extract frequent closed itemsets: e.g. CLOSE [START_REF] Pasquier | Efficient mining of association rules using closed itemset lattices[END_REF], CLOSET [START_REF] Pei | Closet: An efficient algorithm for mining frequent closed itemsets[END_REF], CHARM [START_REF] Zaki | Charm: An efficient algorithm for closed itemset mining[END_REF], TI-TANIC [START_REF] Stumme | Efficient Data Mining Based on Formal Concept Analysis[END_REF], and ZART [START_REF] Szathmary | ZART: A Multifunctional Itemset Mining Algorithm[END_REF]. The set of frequent closed itemsets may be used to determine the set of all frequent itemsets: closed itemsets are a loss less representation of frequent itemsets, while the set of closed itemsets can be orders of magnitude smaller than the set of all frequent itemsets.

FCA organizes information into a concept lattice representing inherent structures existing in data. A concept lattice can be visualized for analysis using graphical tools, e.g. TOSCANA [START_REF] Vogt | Toscana -a graphical tool for analyzing and exploring data[END_REF], GALICIA [START_REF] Valtchev | Galicia: an open platform for lattices. in using conceptual structures[END_REF]. FCA forms also the basis of a knowledge processing paradigm known as "Conceptual Landscapes" [START_REF] Wille | Conceptual landscapes of knowledge: a pragmatic paradigm for knowledge processing[END_REF]. Furthermore, Stumme [START_REF] Stumme | Efficient Data Mining Based on Formal Concept Analysis[END_REF] has introduced the so-called iceberg lattices, which are concept lattices of frequent closed itemsets. Iceberg lattices serve as a support for visualization of association rules mined in large database. They can help analysts in selecting interesting patterns and organizing these patterns into understandable and reusable knowledge units. However, iceberg lattices may hide non frequent but still relevant concepts.

Following the idea of compact, reduced (loss less), and concise representation of extracted units (i.e. itemsets, association rules, or concepts), a number of numerical measures used for pruning itemsets, association rules, and in a certain sense concepts, have been proposed [START_REF] Guillet | Quality Measures in Data Mining (Studies in Computational Intelligence)[END_REF]. In this way, Kuznetsov has introduced stability as a new interest measure for concepts [START_REF] Kuznetsov | Stability as an estimate of the degree of substantiation of hypotheses derived on the basis of operational similarity[END_REF][START_REF] Kuznetsov | On stability of a formal concept[END_REF]. Stability has been successfully used for pruning concept lattices, e.g. in the field of social networks [START_REF] Roth | Co-evolution in epistemic networks -reconstructing social complex systems[END_REF][START_REF] Roth | Towards concise representation for taxonomies of epistemic communities[END_REF][START_REF] Kuznetsov | Reducing the representation complexity of lattice-based taxonomies[END_REF]. Accordingly, in this article, we address the problem of exploring "social communities". By "social communities", we intend sets of agents or organizations whose members are linked by a common interest or objective [START_REF] Ahuja | Network structure in virtual organizations[END_REF]. One of our goals is to study the basis and to design a decision support system for assisting experts identifying social communities. The selection, organization, and discrimination of relevant units of knowledge, help to understand how agents interact in a social community and how they gather on specific topics. Moreover, we show in this paper that combining concept frequency together with concept stability provides a very efficient means for discovering and analyzing social communities.

The paper is organized as follows. Following the present first section, the second section introduces the definitions and the properties of FCA, of support and stability measures. The third section presents a qualitative discussion on stability and shows how stability enlighten concept with a high internal cohesion, i.e. stable and without exceptional individuals. Then, the fourth section gives details on an example of social community discovery within a healthcare network. A discussion on the example and on the knowledge units that can be extracted is proposed and precedes the conclusion of the paper.

Support and Stability: Interest Measures of Formal Concepts

Formal Concept Analysis

We describe here the FCA basics. FCA starts with a formal context K = (G, M, I) where G is a set of objects, M is a set of attributes, and the binary relation I = G × M specifies which objects have which attributes. Two operators, both denoted by , connect the power sets of objects 2 G and attributes 2 M as follows:

:

2 G → 2 M , X = {m ∈ M|∀g ∈ X, gIm}
The operator is dually defined on attributes. The pair of operators induces a Galois connection between 2 G and 2 M . The composition operators are closure operators: they are idempotent, extensive and monotonous. For any A ⊆ G and B ⊆ M, A and B are closed sets whenever A = A and B = B .

A formal concept of the context K = (G, M, I) is a pair (A, B) ⊆ G × M where A = B and B = A. A is called the extent and B is called the intent. A concept

(A 1 , B 1 ) is a subconcept of a concept (A 2 , B 2 ) if A 1 ⊆ A 2 (which is equivalent to B 2 ⊆ B 1
) and we write (A 1 , B 1 ) ≤ (A 2 , B 2 ). The set B of all concepts of a formal context K together with the partial order relation ≤ forms a lattice and is called concept lattice of K.

Iceberg Concept Lattices

This paragraph is based on [START_REF] Stumme | Efficient Data Mining Based on Formal Concept Analysis[END_REF] and introduces basics of iceberg lattices.

Definition 1. Let B ⊆ M. The support count of the attribute set B in K is σ(B) = |B | |G| (1)
Let minsupp be a threshold ∈ [0, 1], then B is said to be a frequent itemset if σ(B) minsupp.

A concept is called frequent concept if its intent is frequent.

Definition 2.

The set of all frequent concepts of a context K is called iceberg lattice of the context K.

The support function is monotonously decreasing: given two attribute sets B 1 and B 2 ,

B 1 ⊆ B 2 ⇒ σ(B 1 ) ≥ σ(B 2 )
. Thus an iceberg lattice is an order filter of the whole concept lattice and in general only a join-semi-lattice. Meanwhile, adding a bottom element makes it a lattice again. Iceberg Lattices can be used to discover and visualize association rules. Within a formal context K = (G, M, I), the task of mining association rules is to determine all pairs

X → Y of M such that σ(X → Y) = σ(X ∪ Y) minsupp, and the confidence conf(X → Y) = σ(X∪Y) σ(X) is above a given threshold minconf ∈ [0, 1]
. Mining associations rules with FCA has two major advantages [START_REF] Lakhal | Efficient mining of association rules based on formal concept analysis[END_REF]. First, frequent closed itemsets are sufficient to deduce all frequent itemsets. Thus, algorithms can benefit from this property to reduce the search space. Second, iceberg lattices offer a reduced and lossless representation of association rules. They allow to directly read Luxenbourger basis for approximate association rules [START_REF] Luxenbouger | Implications partielles dans un contexte[END_REF] from a line diagram.

Stability

Stability has been introduced (probably for the first time) in [START_REF] Kuznetsov | Stability as an estimate of the degree of substantiation of hypotheses derived on the basis of operational similarity[END_REF] and then revisited [START_REF] Kuznetsov | On stability of a formal concept[END_REF][START_REF] Kuznetsov | Reducing the representation complexity of lattice-based taxonomies[END_REF]. Here, we rely on the definition given in [START_REF] Kuznetsov | Reducing the representation complexity of lattice-based taxonomies[END_REF].

Definition 3. Let (A, B) a formal concept of B(K). Stability of (A, B) is γ(A, B) = |{C ⊆ A|C = A = B}| 2 |A| (2)
The stability index of a concept indicates how much the concept intent depends on particular objects of the extent. Given a concept (A, B), the stability index measures the number of elements of G that are in the same equivalence class of A, where an equivalence class is defined as follows.

Definition 4. Let X ⊆ G, we denote by X the equivalence class of X where:

X = {Y ⊆ G|Y = X } (3)
Note that when X is closed, any Y in X is a subset of X. Thus, considering a formal concept (A, B), definition 3 can be rewritten as:

γ(A, B) = | A | 2 |A| (4)
Then, the larger the equivalence class of an extent is (wrt to extent size), the more stable the concept is. The idea behind stability is that a stable concept is likely to have a real world interpretation even if the description of some its objects (i.e. elements in the extent) is "noisy". Figure 1 shows an example of stability in a concept lattice. Each concept is labelled by its extent, intent and stability. For example, for the concept ({1, 5, 6}, {a}), we have:

∅ = {a, b, c, d} = {a} {1} = {a} {5} = {a} {6} = {a, b, c} = {a} {1, 5} = {a} {1, 6} = {a} {5, 6} = {a} {1, 5, 6} = {a} Thus γ({1, 5, 6}, {a}) = 6 8 = 0.75.
It can be noticed that stability is (by definition) always between 0 and 1. It can be still noticed that γ(⊥) = 1. This is always true, since for any subset X from the extent of ⊥, X is included in the intent of ⊥.

Computing stability has been shown to be a #P-complete problem [START_REF] Kuznetsov | On stability of a formal concept[END_REF]. Meanwhile, once the concept lattice has been computed, a bottom-up traversal algorithm can efficiently compute stability [START_REF] Roth | Towards concise representation for taxonomies of epistemic communities[END_REF]. Actually, a concept stability depends on the stability of its subconcepts. This can be shown as follows:

γ(A, B) = 1 - X⊂A,X=X γ(X, X )2 |X|-|A| (5)
Proof. For a formal concept (A, B), from (4), we have:

γ (A,B) = | A | 2 |A| Let I A be the set of subintents of A: I A = {X ⊆ A|X = X }. The set of equivalent classes { X |X ∈ I A } forms a partition of 2 A . Thus |2 A | = X∈IA| X | , which gives: | A | = |2 A | - X∈IA,X =A | X | Dividing by |2 A | we obtain: | A | |2 A | = 1 - X∈IA,X =A | X | |2 A | γ (A,B) = 1 - X⊂A,X=X γ(X, X )2 |X|-|A|

A Qualitative Analysis of Stability

Stability and Cohesion

As stated in [START_REF] Kuznetsov | Reducing the representation complexity of lattice-based taxonomies[END_REF], a concept is stable if its intent does not depend much Âon each particular object of the extent. Stability is aimed at measuring how much a concept extent depends on some of its individual members. This may be useful in analyzing a dataset with a concept lattice, having a special attention to social communities. Here, a social community can be thought as a group of agents -human, software, or resource agents-sharing the same interests, or ideas, or needs [START_REF] Roth | Co-evolution in epistemic networks -reconstructing social complex systems[END_REF]. For example, patients visiting the same hospitals with similar medical problems can be identified as a special social community. In the associated formal context, objects correspond to patients and hospital stays correspond to attributes (see hereafter). Always following the line of [START_REF] Kuznetsov | Reducing the representation complexity of lattice-based taxonomies[END_REF], an actual community has to be "internally cohesive" enough: a stable concept continues to be a concept even if a few members stop being members. This means also that a stable concept is resistant to noise and will not collapse when some members will be removed from its extent.

In this way, a stable concept is a meaningful concept, in the sense that it covers a group of objects, that considered together, have a high internal cohesion. The most stable concepts determine the most interesting groups of objects, that constitute their extents.

Stable Concepts Are of High Interest

FCA provides a powerful framework for identifying social communities [START_REF] Roth | Co-evolution in epistemic networks -reconstructing social complex systems[END_REF][START_REF] Freeman | Using galois lattices to represent network data[END_REF][START_REF] Duquenne | Latticial structures in data analysis[END_REF]. Relations between agents and common interests can be modeled within a formal context. The associated concept lattice will allow to discover and identify which agents do share common interests and what are these interests. However, as the size of a formal context increases, the number of formal concepts in the lattice may grow dramatically. In this case, interest measures such as stability and support can reduce the complexity of the analysis of the concept lattice. Filtering concepts by support relies on the assumption that useful knowledge is represented by frequent patterns. But, this is not always true as pointed out in studies on rare itemsets, as in e.g. in [START_REF] Szathmary | Towards rare itemset mining[END_REF], where it is shown that association rules with a low support but a high confidence may be of high interest.

Stability gives an alternative point of view on formal concepts. It indicates the probability of preserving a concept intent while removing some objects of its extent. Considering social communities, stability helps to identify groups of commons interest that dot not entirely depend on some specific agents. As stability is somewhat independent from support, it can be used to discriminate low-support concepts and detect small communities of strongly related agents. Moreover, stability also detects frequent concepts only depending on a small number of objects. For example, considering a lattice composed of the two following concepts

C 1 = ({g 1 , ..., g n-1 , g n }, {m 1 }) and C 2 = ({g 1 , ..., g n-1 }, {m 1 , m 2 })
with n high, then it can be noticed that C 1 depends solely on the object g n . Although C 1 and C 2 have both a support close to 1, stability of C 2 is 1 while stability of C 1 is 1 2 . In terms of social communities, the group of individuals {g 1 , ..., g n } has not a sufficient "internal cohesion" or has not a "real existence".

Hence, stability, together with support, are a convenient means for identifying two types of concepts:

rare stable concepts: concepts with a low support and a high stability, -frequent unstable concepts: concepts with a high support and low stability.

In section 4, this point of view is discussed and illustrated within a real-world application aimed at detecting communities of patients, i.e. groups of patients being treated in the same groups of hospitals.

Stable Concepts Are Monothetic Rather Than Polythetic

As introduced above, stability can also be linked with exceptions, and, furthermore, with the so-called monothetic and polythetic characters of a class of individuals [START_REF] Sokal | Phenetic taxonomy: Theory and methods[END_REF][START_REF] Bailey | Monothetic and polythetic typologies and their relation to conceptualization, measurement and scaling[END_REF][START_REF] Lerman | Les bases de la classification automatique[END_REF]. When building a concept lattice and analyzing groups of individuals through the extents of concepts, one problem is to recognize and explain exceptions. A subsequent question is to understand whether exceptions are linked to monothetic or polythetic classes.

A class of individuals C is said to be monothetic if and only if there exists a set of attributes Att that determines the membership of an individual to the class C (Att is a set of necessary and sufficient membership conditions). By contrast, given a set of attributes Att = {a 1 , ..., a n }, a class of individuals C is said to be polythetic if and only if:

-Every object that is an instance of the class C has an "important" -not necessarily fixed-number of attributes of Att. -Every attribute of Att belongs to an "important" -not necessarily fixednumber of instances of C. -There is not necessarily an attribute of Att belonging to every instance of C.

Relying on the fact that a stable concept has a high "internal cohesion", is resistant to noise, and does not collapse when some members stop being members of its extent, the more a concept is stable, the more it does not represent exceptional individuals, and, accordingly, the more the concept is able to represent cohesive groups of individuals, such as social communities. This means that stable concepts are rather monothetic and that unstable concepts are rather "exceptional" or polythetic, i.e. they include some exceptional character, shared by only a few individuals. For illustrating this view, let us consider the following example.

A B C D E 1 x x x x 2 x x 3 x x x x 4 x x x 5 x x x x
Here, attribute D can be considered as a "necessary and sufficient condition" for the membership of an individual to the class including individual 1. Indeed, the concept lattice includes the concept (A, B) = (1, abde). However, this concept has a very low stability ( 4 16 ) and also a low support for its intent ( 1 5 when abde is considered as an itemset). This is in agreement with the view that a stable concept appears to be monothetic while an unstable concept tends to be polythetic.

4 Social Communities in a Regional Healthcare System

Motivations

Healthcare management and planning play a key role for improving the overall health level of the population. From a population point of view, even the best and state-of-the-art therapy is not effective if it cannot be delivered in the right conditions. Actually, many determinants influence the effective delivery of healthcare services: availability of trained personnel, availability of equipments, security constraints, costs, proximity. . . All of these should meet economics, demographics, and epidemiological needs in a given area. This issue is especially acute in the field of cancer care where many institutions and professionals must cooperate to deliver high level, long term, and costly care. Therefore, it is crucial for healthcare managers and decisions makers to be assisted by decision support systems that give strategic insights about the intrinsic behavior of the healthcare system.

On the one hand, healthcare systems can be considered as "data rich" as they produce massive amounts of data such as electronic medical records, clinical trial data, hospital records, administrative data, and so on. On the other hand, they can be regarded as "knowledge poor" as these data are rarely embedded into a strategic decision-support resource [START_REF] Abidi | Knowledge management in healthcare: towards 'knowledge-driven' decision-support services[END_REF]. In France, the PMSI database is a national information system used to describe hospital activity with both an economical and medical point of view. In a previous work, we used this system together with iceberg lattices to discover how several institutions organize themselves into an implicit network to provide coordinated care at a regional level [START_REF] Jay | Using formal concept analysis for mining and interpreting patient flows within a healthcare network[END_REF]. Our method has been used in real world by healthcare managers. It appeared that support based pruning had some limits, for example in analysing small institutions interactions.

The Difficulty of Choosing the Good Support Threshold

In this section we present an example of an iceberg lattice showing cooperations between hospitals in the field of cancer. We then discuss the choice of minsupp by studying concept support distribution. In our approach, we build a context in which objects are patients suffering from cancer and attributes are hospitals. A patient and a hospital are related if the hospital has delivered cancer care to this patient. In our experiment, the resulting context has 6036 patients and 170 hospitals. While the whole concept lattice holds 865 concepts, an iceberg lattice built with a minsupp of 0.0033 (20 patients) gives 93 frequent concepts. A small excerpt of this iceberg is shown in figure 2.

Fig. 2. Iceberg of cancer treatment cooperations

For the sake of clarity, ⊥ was removed. Although its right and leftmost parts are not drawn, the iceberg lattice is much more wide than deep because the context is sparse and data are poorly correlated. This means that cooperations are most of the time tightly partitioned, and that patients are rarely hospitalized in more than two hospitals. The intent of co-atoms, i.e. immediate descendants of , is always a singleton, indicating that a hospital never shares all of its patients with another one, or if it is so, less than 20 patients are involved in the interaction. The intent of atoms, i.e. the immediate ascendant of ⊥, is always a pair. The extent of atoms gives an idea of the strength of the cooperation between the two hospitals lying in the intent: the larger is the cardinal of the extent, the higher is the strength of the cooperation (i.e. the more patients are shared between the two hospitals). The examination of the iceberg brings different types of knowledge:

some concepts that are both atoms and co-atoms (for example : CL-SELI-THIO). They represent institutions that share a few patients with others. This is that either they treat a few patients, or they work in a relative autonomy, or cooperation is split with many other hospitals. -other concepts have at least a sub-concept (different from ⊥). They represent a hospital receiving a significant number of patients, and having collaborations with at least another establishment. For example MAT-REG and CLCC-AV-VAN share 28 patients. -The concept representing the CLCC-AV-VAN hospital has a high support and many sub-concepts. This hospital is a specialized anti-cancer center. It employs highly skilled and specialized personnel. Treatments given there rely on state-of-art technology. Furthermore, It actively participates in anti-cancer research programs and thus can be considered as a reference institution.

The choice of the minimum support strongly influences the interpretation of the iceberg. It must be sufficiently low to convey meaningful knowledge and sufficiently high to keep this knowledge readable for a human expert. Here, the whole lattice holds 865 concepts. Table 1 shows the distribution of their support according to the size of their intent. By choosing a minsupp of corresponding to 20 patients, we see we can miss interesting knowledge:

Table 1. Concept support count and intent size h h h h h h h h h h h h h h

-13 hospitals treat at least 10 patients, -52 cooperations of two hospitals involve at least 10 patients, -7 concepts represent cooperations between three hospitals and involve at least 10 patients.

In our case, support is a weak mean for discriminating "not-so-frequent" concepts. Some concepts not appearing in the iceberg may be of interest for different reasons:

-They illustrate a cooperation of one hospital sharing almost all of its patients with another one. -They concern a hospital treating few patients but not sharing them with any other one. -They concern a 3-hospitals interaction.

Lowering the support threshold can let these concepts appear in the iceberg, but at the expense of readability. Moreover support measure is not specific enough to discriminate concepts with the above characteristics.

Stability Analysis

In this section, we study stability of concepts within the whole concept lattice. Figure 3 shows the histogram of concept stability.

Most of the concepts hold only one object in their extent. As ⊥ extent is empty, they have a stability of 0.5. The next important group consists of concepts having stability very close to 1. It corresponds generally to the most frequent concepts. Indeed, 93% of the concepts having support greater than 10 have stability greater than 0.99. Figure 4 is a scatter-plot of stability and support, also featuring intent size.

Values are presented on a log-scale for better visualization of both support and stability ranges. Moreover, we prefer to use the stability odds defined as follows:

o γ (A, B) = γ(A, B) 1 -γ(A, B) (6)
Stability odds are the ratio between the number of subsets of an extent A which belong to the equivalence class of A, and of those which do not.

Stability odds illustrate more clearly the distribution of stability for values close to 1 and show that stability has a better discriminant power than support. The whole set of concepts has 87 distinct support values and 152 distinct stability values. For concept having support between 10 and 19, we observe 57 distinct stability values. Figure 4 displays thresholds for the 10% most frequent concepts and 10% most stable concepts. In the next sections we discuss the differences between those two thresholds.

Frequent Unstable Concepts

The upper left quarter of figure 4 shows the 9 concepts having support greater than the 10% most frequent threshold and stability less than the 10% most stable threshold, i.e. frequent unstable concepts. They represent institutions playing a significant role in cancer care as secondary or tertiary care1 centers. They share most of their patients with other hospitals. This gives rise to cooperations of 2 kinds.

Figure 5 shows the order ideals of two frequent unstable concepts. Two hospitals have tight interaction each with another tertiary care center, i.e. MAT-REG-NANCY with CRLCC-AV-VAN, CL-ARCENCIEL with CH-EPINAL. These two institutions are secondary care centers usually delivering surgery and referring patients to a tertiary care center for radiotherapy. Their activity can be almost entirely explained by an exclusive cooperation. This induces a form of dependency with the tertiary care center. Patient not concerned by the cooperation can be considered as exceptions, i.e. patients not following the usual care pathway for some reason (for example because they suffer from a very specific pathology). After pruning according to stability, only the concepts in grey on figure 5 will remain.

One center, as shown on Figure 6, is a highly specialized anti-cancer institution. IGR-PARIS is an international well-known anti-cancer center located in Paris. Many Lorraine local hospitals refer patients to IGR for rare tumors. If we apply stability pruning here, the whole sub-lattice will disappear, which may be desirable for readability. Meanwhile, its suggests that searching for this type of unstable frequent concept can also be an interesting knowledge mining task.

Rare Stable Concepts

Three concepts are located in the right lower corner of figure 4.

Two have an intent of size 2 and illustrate thus cooperations between two hospitals sharing 18 patients. These cooperations differ from others of the same support in that they are not split themselves in cooperations involving a third institution. Thus, similarity of concerned patients is entirely explained by those cooperations. The last concept has a size 3 intent and illustrates the cooperation between three large specialized hospitals located in the same city of Nancy. This is the most frequent and most stable concept of that type.

Stability allows to distinguish rare concepts that cannot be separated from others by support. Rare stable concepts differ from other rare concepts in that their attributes suffice to explain the similarity of their objects.

Synthesis

Our system is used in real world by healthcare managers. First, at the individual level, it helps healthcare professionals to assess their activity in the regional landscape. While physicians are able to cite the names of people they are used to cooperate with, they cannot measure the strength of these cooperations. And it is even harder for hospital managers to count patients shared with other institutions due to the gaps and lacks of adapted processes in information systems. Second, at the regional level, it provides for the administrative staff a decision support to reorganize care resources according to the implicit behavior of the healthcare system. Actually, French law establishes activity thresholds in the field of cancer: e.g. an hospital must treat at least 30 patients a year to be authorized for digestive cancer surgery. Our system has allowed to enlighten and accordingly to promote cooperations between institutions that could not reach the thresholds alone. It has also demonstrated how administrative decisions could impact the social healthcare network, given the existence of many dependencies between structures.

Conclusion

Our main objective is to build a system allowing for visualization of a social healthcare network in the field of cancer care. This system is used today by healthcare managers. Things must be kept simple while conveying enough information for assisting strategic decisions. The use of concept interest measures has a strong impact both on readability and semantics of discovered knowledge. Together with support, stability can successfully identify two kind of concepts: frequent unstable concepts and rare stable concepts. In our experiment, the formal context is sparse and we need to mine concepts with very low support. Stability brings additional knowledge that helps to discover interesting rare concepts that can not be discriminated by support. We believe that it could have significant implications in the field of rare itemsets mining [START_REF] Szathmary | Towards rare itemset mining[END_REF]. Furthermore, stability enhances lattice visualization when pruning frequent unstable concepts. Besides, frequent unstable concepts may also be a subject of interest.

Fig. 1 .

 1 Fig. 1. Stability example

Fig. 3 .Fig. 4 . 5 .

 345 Fig. 3. Concept stability histogram

Fig. 5 .Fig. 6 .

 56 Fig. 5. Frequent unstable concepts: exclusive cooperations

Secondary care is delivered by a broadly skilled specialist (e.g. a general surgeon, a general internist, or an obstetrician). Tertiary care is provided by a sub-specialist(e.g. an orthopedic surgeon, a neurologist, or neonatologist).