Samuel Kadoury

Nikos Paragios

Towards Shape Constellation Inference through Higher-Order MRF Optimization in Nonlinear Embeddings

Keywords: Image inference, statistical modeling, non-linear manifold, articulated model, Markov Random Field Inférence d'images, modélisation statistique, domaines non-linéaires, modèle articulé, Markov Random Field

This paper introduces a novel approach for inferring articulated shape models from images. A low-dimensional manifold embedding is created from a training set of prior mesh models to establish the patterns of global shape variations. Local appearance is captured from neighborhoods in the manifold once the overall representation converges. Inference with respect to the manifold and shape parameters is performed using a Markov Random Field (MRF). Singleton and pairwise potentials measure the support from the data and shape coherence from neighboring models respectively, while higher-order cliques encode geometrical modes of variation in localized shape models. Optimization of model parameters is achieved using efficient linear programming and duality. The resulting model is geometrically intuitive, captures the statistical distribution of the underlying manifold and respects image support in the spatial domain. Experimental results on articulated bodies such as spinal column geometry estimation demonstrate the potentials of our approach.

Introduction

Statistical models of shape variability have been successful in addressing fundamental vision tasks such as segmentation and registration in computer vision and medical imaging analysis. Such models help to understand the distribution in appearance of a group of shapes and offer efficient parametrization of the geometric variability in a given cluster. In medical imaging, one has the ability to compare information of a subject over a time period, between different subjects or between groups of pathological cases by using these types of representation [START_REF] Rueckert | Automatic construction of 3-D statistical deformation models of the brain using nonrigid registration[END_REF]. These models have been used extensively for localized structures. On the other hand object constellations have been sparely investigated in computer vision (mainly on body [START_REF] Andriluka | Pictorial structures revisited: People detection and articulated pose estimation[END_REF] and hand-pose estimation [START_REF] De La Gorce | Model-based hand tracking with texture, shading and self-occlusions[END_REF]), which is mostly due to the challenges in interpreting the images.

Various statistical approaches for three-dimensional modeling of structures have been proposed based on shape analysis. Active Shape (ASM) [START_REF] Cootes | Active shape models -their training and application[END_REF] and Appearance Models (AAM) [START_REF] Cootes | Active appearance models[END_REF] have been successful in recovering object geometries obtained from dense collection of data points. Implicit representations is an alternative formulation [START_REF] Rousson | Prior knowledge, level set representations and visual grouping[END_REF] to address model-based segmentation while more recently numerous methods based on point distribution models and embedding on various geometric spaces (spherical [START_REF] Nain | Multiscale 3-D shape representation and segmentation using spherical wavelets[END_REF]) have been proposed. However, model-based segmentation of single objects typically leads to fitting errors when no clear object boundary is present, similar structures in close vicinity, or pathologies. The result is thereby sensitive to model initialization and therefore limited to the capture range.

Simultaneous multi-object inference is often beneficial compared to the separate segmentation of individual objects. In [START_REF] De Bruijne | Quantitative vertebral morphometry using neighbor-conditional shape models[END_REF], an extension of point distribution models (PDMs) was considered for modeling relations between shape constellations using conditional probabilities between 2D-contours. Furthermore, rigid transformations were considered as statistics between different parts to evaluate the variability from a given class of models [START_REF] Boisvert | Geometric variability of the scoliotic spine using statistics on articulated shape models[END_REF] and to achieve multi-modal inference [START_REF] Kadoury | Surface/volume-based articulated 3D spine inference through Markov Random Fields[END_REF].

Last, but not least in the context of medical imaging, traditional body-pose estimation methods are prone fail. This is due to the fact that the aim of articulated models is to deduce the anatomical structure but most importantly to identify/segment the individual structures precisely, which is not the focus of vision methods. In vision where mostly 2D models are used, classification techniques are often considered to provide candidate positions for the different parts [START_REF] Felzenszwalb | Pictorial structures for object recognition[END_REF]. This information along with prior constraints are used to determine the most appropriate pose configuration.

The main limitation of the above mentioned methods is that incorporating imagedriven information is not straightforward. Furthermore, the high dimensionality and complex non-linear underlying structure unfortunately makes the commonly used linear statistics [START_REF] Cootes | A unified framework for atlas matching using active appearance models[END_REF] inapplicable for articulated structures. Given the important and often un-correlated variations in the orientation and translation of the global shape, such statistics can only be performed along a linear manifold that appropriately parameterizes the nonlinear space. Recent concepts on Riemannian manifolds in tensor spaces is particulary suited for this application [START_REF] Khurd | Manifold Learning Techniques in Image Analysis of High-dimensional Diffusion Tensor Magnetic Resonance Images[END_REF].

A manifold learning algorithm of particular interest to this work is locally linear embedding [START_REF] Roweis | Nonlinear dimensionality reduction by locally linear embedding[END_REF]. It maps high-dimensional observation data that are presumed to lie on a nonlinear manifold, onto a single global coordinate system of lower dimensionality. Therefore, it helps to preserve neighborhood relationships of similar object geometries, thereby revealing the underlying structure of the data which can be used for statistical modeling. Inferring a model from the underlying manifold is a novel concept but far INRIA from being trivial, and relies on a cost function which includes visual support and prior constraints.

In this paper, we introduce a deformable articulated body instantiation through a statistical modeling of inter-object transformations. We use nonlinear manifold embeddings which can handle both small and large deformations of the object constellations in a given dataset. Our principal contribution lies with regards to the representation of the model. Rather than directly optimizing a set of transformation parameters, we perform a statistical embedding on the transformation vectors which describe the articulations, thus significantly reducing the degrees of freedom. An explicit representation of the nonlinear manifold given by smooth functions enables the joint mapping between ambient and embedded space. Furthermore, because of the limited dimensionality of the search space, the inference with respect to the low-dimensional parameters is performed in an integrated and interconnected high-order MRF. This graph involves costs related both to the data, prior geometrical dependencies and global higher-order cliques. We consider a framework which applies recent advances in MRF labeling schemes generalizing energy function classes with degrees higher than two which can be decomposed into pair-wise potentials [START_REF] Rother | Minimizing sparse higher order energy functions of discrete variables[END_REF]. This demonstrates the capability to model complex interactions of random variables. Recent advances in the area of discrete optimization which explore duality theorem of linear programming [START_REF] Komodakis | Performance vs computational efficiency for optimizing single and dynamic MRFs: Setting the state of the art with primal-dual strategies[END_REF] are explored to obtain the lowest potential of the objective function.

This paper is organized as follows. Section 2 presents the theoretical methodology for the manifold embedding algorithm of articulated mesh models while in Section 3 we propose the MRF-based inference framework minimizing P n Potts models energy terms. In Section 4 we present our evaluation results applied to medical images of the spinal column and the last section concludes the paper.

Manifold Embedding of Articulated Models

The input to our method is a sample of articulated models which comprises a set of learning shapes. These shapes are a constellation of triangular meshes, each annotated with landmarks defined as characteristic points uniquely localized across a set of objects. We first build an articulated shape model from a training database by embedding the data into a low-dimensional sub-space which dimensionality corresponds to the domain of admissible variations. Local shape appearance is determined via analysis of variations within a sub-patch of the manifold.

Representation of Articulated Models

Our shape model S = {s 1 , . . . , s L } consists of an interconnection of L objects. For each local shape s i , we recover a triangular mesh with vertices {v i j |j = 1, . . . , V }, where the j th vertex corresponds to approximately the same surface location from one shape to another. Additionally, every s i is annotated with landmarks on each personalized model to rigidly register each object to its upper neighbor. The resulting rigid transforms are stored for each inter-object link. These transforms can also be determined via an ICP-like algorithm to recover the extrinsic parameters. Hence, an articulated deformable model (ADM) is represented by a vector of local inter-object rigid transformations A = [T 1 , T 2 , . . . , T L ] as illustrated in [Fig. 1]. To perform global shape modeling of the shape S, we convert A into an absolute representation

A abs = [T 1 , T 1 • T 2 , . . . , T 1 • T 2 • . . . • T L ] (1)
using recursive compositions. The transformations are expressed in the local coordinate system (LCS) of the lower object, which can be defined proper to the object's geometrical representation capturing main axes of deformation. Center of transformation is located at the midpoint of the triangular mesh. The rigid transformations described in this paper are the combination of a rotation matrix R, a translation vector t and the scaling s. We formulate the rigid transformation T = {s, R, t} of a triangular mesh model as y = sRx + t where x, y, t ∈ 3 . Composition is given by

T 1 • T 2 = {R 1 R 2 , R 1 t 2 + t 1 }, while inversion as T -1 = {R T , -R T t}.

Nonlinear Manifold Embedding of ADMs

Let us consider N articulated shape models expressed by the absolute vector representation

A i abs , of dimensionality D. The aim is to create a low-dimensional manifold consisting of N points Y i , Y i ∈ d , i ∈ [1, N ] where d D.
In such a framework, if an adequate number of data points is available, then the underlying manifold M is considered to be "well-sampled". Therefore, it can represent the underlying population structure. In the sub-cluster corresponding to a pathological population, each individual point of the training set and its neighbours would lie within a locally linear patch on the manifold.

Nearest neighbor selection

The main limitation of embedding algorithms is the assumption of Euclidean metrics in the ambient space to evaluate similarity between sample points. In our approach, we adopt the intrinsic nature of the Riemannian manifold geometry allowing us to discern between articulated shape deformations in a topological invariant framework. The K closest neighbours are selected for each point using a distortion metric particularly suited for geodesic metrics, defined as d M (A i abs , A j abs ) which estimates the distance of articulated models i, j. As such, A i abs and A j abs are represented by the feature vectors described in [START_REF] Rueckert | Automatic construction of 3-D statistical deformation models of the brain using nonrigid registration[END_REF]. The distance measure can therefore be expressed as a sum of articulation deviations:

d M (A i abs , A j abs ) = L k=1 d M (T i k , T j k ) (2)
= L k=1 d M R i k c i k ... 1 , R j k c j k ... 1
where the canonical representation encodes the intrinsic (c) and orientation (R) parameters. The difference between analogous articulations is computed within the geodesic framework:

d M (A i abs , A j abs ) = L k=1 c i k -c j k + L k=1 d Q (R i k , R j k ). (3)
The first term evaluates intrinsic distances in the L 2 norm. Using the geodesics, it is possible to define a diffeomorphism between rotation neighborhoods in M and a tangent plane T x M. The exponential map at x ∈ M maps vectors of the tangent plane T x M to a point in the manifold which is reached by the geodesic γ x,v in a unit time. In other words, if γ (x,v) (1) = y, then the inverse mapping is known as log x (y) = v. The quaternion distances are therefore computed with the Frobenius

norm d Q (R i k , R j k ) = log((R i k ) -1 R j k )
F based on the geodesic distance on the 3D manifold. This is feasible since rotations R i k , R j k are nonsingular, invertible matrices. The log map can be thought as the unfolding operation connecting the tangent space at RT n°0376

x and the manifold M. One can now proceed to the manifold reconstruction using the local support in high-dimension data.

Embedding algorithm

The manifold reconstruction weights are estimated by assuming the local geometry of the patches can be described by linear coefficients that permit the reconstruction of every model point from its neighbours. In order to determine the value of the weights, the reconstruction errors are measured using the following objective function:

min W ε(W ) = N i=1 A i abs - K j=1 W ij A j abs 2 (4)
subject to

W ij = 0 if A i abs not neighbor A j abs j W ij = 1 for every i, (5)
Here, A i abs is the absolute vector describing an articulated model described above and ε(W ) sums the squared distances between all data points and their corresponding reconstructed points. The weights W ij represent the importance of the j th data point to the reconstruction of the i th element. It is easy to show that each weight can be calculated individually [START_REF] Roweis | Nonlinear dimensionality reduction by locally linear embedding[END_REF]. Each sample A i abs and weight W ij contributes to ε such that:

ε (i) (W ) = K j=1 K m=1 W ij W im Q i jm (6) with Q as a K ×K matrix s.t. Q i jm = (A i abs -A j abs ) T (A i abs -A m abs ). With R = (Q+ ρI) -1
, given a suitably chosen regularization constant ρ (see [START_REF] De Ridder | Locally linear embedding for classification[END_REF]), then W ij are solved by a least squares problem given the constraint in [START_REF] Cootes | Active appearance models[END_REF],

W ij = ( K m=1 R jm )/( K p,q=1 R pq ).
We propose to solve Q based on distance matrix D using articulated metric in [START_REF] De La Gorce | Model-based hand tracking with texture, shading and self-occlusions[END_REF]

such that Q = (D ij + D im -D jm )/2.
The algorithm then maps each high-dimensional A i abs to a low-dimensional Y i . These global internal coordinates are determined with a cost function minimizing the reconstruction error

Φ(Y ) = N i=1 Y i - K j=1 W ij Y j 2 (7) = N i=1 N j=1 M ij Y T i Y j
with M as a sparse and symmetric N × N matrix enclosing the reconstruction weights W ij such that M = (I -W) T (I -W), and Y spanning the Y i 's. By constraining Y to identity, the problem becomes a straightforward one with min tr(YMY T ) given the constraint 1 N YY T = I. Using Lagrange multipliers and setting the derivative to zero gives (M -Λ)Y T = 0, where Λ is the diagonal Lagrange multiplier matrix. The optimal embedding, up to a global rotation, is obtained from the bottom d + 1 eigenvectors of the and helps to minimize the cost function Φ(Y ) as a simple INRIA eigenvalue problem. The d eigenvectors form the d embedding coordinates. The coordinates Y i can be translated by a constant displacement without affecting the overall cost Φ(Y ). The eigenvector corresponding to the smallest eigenvalue corresponds to the mean value of the embedded data Y 0 = {y 1 , . . . , y d }, s.t. y i = 0, ∀i. This can be discarded with Y i = 0 to obtain an embedding centered at the origin. Hence, a new ADM denoted as Y new can be inferred in the embedded d-space as a low-dimensional data point by finding its optimal manifold coordinates y i .

Inverse mapping to the ambient space

To obtain the articulation vector for a new embedded point in the ambient space (image domain), one has to determine the representation in high-dimensional space based on its intrinsic coordinates. We first assume an explicit mapping f : M → D from manifold M to the ambient space D . The inverse mapping of Y i is then performed by estimating the relationship between D and M as a joint distribution, such there exists a smooth functional which belongs to a local neighborhood. Theoretically the manifold should follow the conditional expectation:

f (Y i ) ≡ E(A i abs |M(A i ) = Y i ) = A i f (Y i , A i ) f M(Ai) (Y i ) dA (8)
which captures the overall trend of the data. Here, both f M(Ai) (Y i ) (marginal density of M(A i )) and f (Y i , A i ) (joint density) are unknown. Using Nadaraya-Watson kernel regression theorem [START_REF] Nadaraya | On estimating regression[END_REF], densities are replaced by kernel functions based in a conditional expectation setting [START_REF] Davis | Population shape regression from random design data[END_REF]. We replace these variable with f M(Ai)

(Y i ) = 1 K j∈N (i) G h (Y i , Y j ) and f (Y i , A i ) = 1 K j∈N (i) G h (Y i , Y j )G g (A i , A j ).
The Gaussian regression kernels G require the neighbors A j abs of j ∈ N (i) to determine the bandwidths h, g so it includes all K data points. Plugging these estimates in [START_REF] De Bruijne | Quantitative vertebral morphometry using neighbor-conditional shape models[END_REF], this gives:

f NW (Y i ) = A i 1 K j∈N (i) G h (Y i , Y j )G g (A i , A j ) 1 K j∈N (i) G h (Y i , Y j ) dA. (9)
By assuming G is symmetric about the origin and generalize the expectation such that the observations Y are defined in terms of a metric d M in manifold space M, we obtain:

f NW (Y i ) = argmin A i abs j∈N (i) G(Y i , Y j )d M (A i abs , A j abs ) j∈N (i) G(Y i , Y j ) (10)
which integrates the distance metric d M (A i abs , A j abs ) defined in (3) and updates f NW (Y i ) using the closest neighbors of point Y i in the manifold space. This constrains the regression to be valid for similar data points in its vicinity since locality around Y i preserves locality in A i abs .

Local Shape Appearances in the Manifold

The key idea of capturing local shape appearance lies on the assumption that models, represented in a given neighborhood M k of the general manifold M s.t. M k ∈ M : d , will also manifest similar local geometries. We assume here that local shape appearances follow a linear distribution within the low-dimensional manifold. Hence, given a data point Y j and its K neighbors, the local shape model s i , representing the RT n°0376 i th element of the ADM, is obtained by building a particular class of shapes given the set of examples {s 1 i , ..., s K i }. This set belongs to the sub-patch M k , such that each shape s j i ∈ S j . We approximate the distribution of the shape using a parameterized linear model by computing the deformation vectors formed for the K -1 shape samples. We compute the n eigenvalues λ 1 , λ 2 , . . . , λ n and corresponding eigenvectors, v 1 , v 2 , . . . , v n of the i th covariance matrix:

C i = 1 K -1 K j=1 (s j i -si )(s j i -si ) T ( 11
)
where si is the mean shape of the K neighboring local objects si = 1 K K j=1 s j i . By arranging eigenvalues in increasing order, a new instance s new i of the shape model can be defined. The new shape is approximated by a linear combination of the deformation vector as follows:

s new i ≈ si + n i=1 ω i v i = si + [v 1 . . . v n ][ω 1 . . . ω n ] (12)
= si + Vw where V is a K × n matrix with eigenvectors in each column, and w = [ω 1 . . . ω n ] the weight vector of size n. In ( 12), w provides a compact parametrization of the transformation we seek for each object composing the model shape S based on the principal deviations. This final step warps individual instances from the training set to infer new local shape models.

Inference through MRF Optimization

Once an appropriate modeling of shape variations is determined, a successful inference between the image and manifold must be accomplished. We describe here how a new model is deformed, the similarity criterions as well as the adopted optimization procedure which infers the shape model to the sparse data. We search the optimal embedded manifold point Y = (y 1 , . . . , y d ) of the global ADM. Such a strategy offers an ideal compromise between the prior constraints, as well as the individual shape variation described by the main weight vectors W = (w 1 , . . . , w n ) in a localized sub-patch. Formally, the inference of the model S to the image I is given by:

({y 1 , . . . , y d }; {w 1 , . . . ,

w n }) = argmin δi,ωi E(S 0 , I, ∆, Ω). (13)
The energy E of inferring the model S in the image I is a function of the displacement vectors ∆ = ( δ 1 , . . . , δ d ) in the non-linear manifold space for global shape representation. This involves a data-related term V (Y 0 + ∆, I) expressing the image cost and a global prior term V (N, ∆) measuring deformation between low-dimensional vectors with shape models. Furthermore we introduce a higher-order term V (H, ∆, Ω) which is expressed by the reconstruction weights Ω = (ω 1 , . . . , ω n ) for local shape modeling in a linear space. The energy function E can therefore be decoupled to a global and local optimization scheme controlled by the weighting parameters α-β:

E S 0 , I, ∆, Ω = V Y 0 + ∆, I + α V N, ∆ (14)
+ β V H, ∆, Ω .
We explain in this section how we define the global and local energy terms.

Rigid Alignment of the ADM

The global alignment of the model with the target (i.e., image data) primarily drives the deformation of the ADM in the first phase of convergence controlled by the α-term such that it is minimized when the neighbors K in a manifold sub-patch stabilizes.

The purpose here is to estimate the set of articulations describing the global model shape by determining its optimal representation Y 0 in the embedded sub-space. This is performed by reformulating the global representation using the inverse mapping in ( 10)

so f NW (Y i + ∆) = f NW ({y 1 + δ 1 , . . . , y d + δ d }) = A i
abs + D with D as deformations in ambient space D . This allows to represent the model in image space based on its manifold space coordinates of M. The global cost is expressed as:

V Y 0 + ∆, I = V f NW ({y 1 + δ 1 , . . . , y d + δ d }), I) = V A 0 abs + D, I . ( 15
)
Since the transformations T i are implicitly modeled in the absolute representation A 0 abs , we can formally consider the singleton image-related term as a summation of costs RT n°0376 associated with each L objects of the ADM:

V A 0 abs + D, I = L i=1 V i s i * (T 0 i + d i ), I (16)
where

V i (s, I) = vi∈s n T i (v i )∇I(v i ),
is a modular term seeking to minimize the distance between the mesh vertices of the inferred ADM and image data I by a rigid transformation of the vectors. This globally modifies the mesh vertices in order to align them with the I. In ( 16), v i is the set of triangle centers, n i is the triangle normal pointing outwards, and ∇I(v i ) is the image gradient at location v i . This modular term effectively measures the strength of the edges over the triangles corresponding to the inferred model. The prior constraint for the rigid alignment is defined in the manifold space by modeling pairwise potentials between low-dimensional features y i represented by the second term in ( 14):

α V N, ∆ = α ij i∈G j∈N (i) V ij (y 0 i + δ i , y 0 j + δ j ). ( 17
)
This potential measures the distance between pairs y 0 i and y 0 j from the current data point coordinates to a prior distribution built with the K points y 1,..,K i , y 1,..,K j in the local neighborhood M k . This term represents the smoothness term of the global cost function to ensure that the deformation δ i applied to point coordinates is regular in the non-linear vicinity of variations. This is measured by the geodesic distance in the Riemannian manifold

V ij = P (y 0 i + δ i ) -(y 0 j + δ j ) F , N (i) . (18)
The Frobenius distance is measured to a density function determined from the the K samples in the neighborhood such that P (x, N ) =

1 σ N √ 2π (1 -e (-x-µ N ) 2 /2σ 2 N )
assigns a cost based on the probability that x belongs to a Gaussian distribution, with µ and σ calculated from the i and j vectors of N (i) ∈ M k .

Non-Rigid Adaptation of Local Shapes

Local shape geometry for each of the ADM's components is obtained by varying the weight parameters of the principal variations. We parameterize these potentials with clique variables w c taking on corresponding costs θ q if the cliques are given to the weight vectors ω c . Hence the third term of ( 14) is described as a high-order functional:

β V H, ∆, Ω = β c c∈C V c (w 0 c + ω c ) (19)
where independent clique variables c are treated as a graph minimization problem. Our prior term is represented by higher-order potentials of degree n, based on the eigenvalues in 2.3 of the L local objects from our model S. The potential functions are defined as:

V c (w 0 c ) = min{ min q∈{1,2,...,t} θ q + ∆ q (w 0 c ), θ max } (20)
INRIA with ∆ q (w 0 c ) = i∈c;l∈L ψ q il κ(ω i = l) is a deviation function to evaluate the assignment of w 0 c . The Kronecker delta function κ generates binary variables, while the weights are assigned such that ψ q il = {0, θ max } depending on whether the clique variable w c is given the appropriate label.

Our work is inspired from a mesh reconfiguration where the costs θ q depend on the positions of the morphed mesh vertices v i . It introduces a tradeoff between adherence to the image and coherence with regards to the prior principal modes of shape variations in the neighborhood M k . The strategy finds the point in the target image with an optimal boundary criterion. Hence for each mesh triangle, a search along the normal n i from v i with optimal compromise between boundary detection and distance jγ is performed. It measures the cost of driving the triangle centers towards detected features points in the target image [START_REF] Kaus | Automated 3-D PDM construction from segmented images using deformable models[END_REF]:

θ q = V i=1 w q ∇I(v i ) ∇I(v i ) γn i argmin j=-m,...,-m (j 2 γ 2 (21) -n T i ∇I(v i + jγn i )) 2 .
In (21), V is the number of triangles, m is the sample space to find the closest detected feature in ∇I and γ specifies the distance between two points on the normal. The weights w q controls the tradeoff of how far the directed gradient distance lies from the mean eigenvalue shape.

One can integrate the global data and prior terms along with local shape terms parameterized as the higher-order cliques, by combining [START_REF] Rother | Minimizing sparse higher order energy functions of discrete variables[END_REF], [START_REF] De Ridder | Locally linear embedding for classification[END_REF] and ( 19): E S 0 , I, ∆, Ω = V f NW ({y 1 + δ 1 , . . . , y d + δ d }), I)

+ α ij i∈G j∈N (i) V ij (y 0 i + δ i , y 0 j + δ j ) + β c c∈C V c (w 0 c + ω c ). (22)
RT n°0376

Energy Minimization

The optimization strategy of the resulting MRF (22) in the continuous domain is not a straightforward problem. The convexity of the solution domain is not guaranteed, while gradient-descent optimization approaches are prone to non-linearity and local minimums. We therefore considered results obtained from discrete optimization approaches [START_REF] Komodakis | Performance vs computational efficiency for optimizing single and dynamic MRFs: Setting the state of the art with primal-dual strategies[END_REF]. By assuming the initial zero-mean shape model in the given manifold space M falls in the image domain I and that the desired accuracy is specified so to define the quantization step, we can approximate the deformation of the shape model towards the optimal solution. We seek to assign the optimal labels L ∆ = {l 1 , . . . , l d } and L Ω = {l 1 , . . . , l n } which are associated to the quantized space ∆ of displacements and local weight parameters Ω. If we consider that displacing the coordinates of the sub-domain point y 0 i by δ li is equivalent to assigning label l i to y 0 i , we can formulate the energy as:

({l ∆ 1 , . . . , l ∆ d }; {l Ω 1 , . . . , l Ω n }) = argmin l ∆ i ,l Ω i ∈L ∆ ,L Ω E(S 0 , I, L ∆ , L Ω ). (23)
The quantization domain of the label set is crucial to achieve both inference accuracy and computational efficiency. To this end, we adopt a coarse-to-fine approach which continuously increases the number of displacements while decreasing the search space. This speeds up the convergence and helps to reduce the interaction between graph nodes. Furthermore to account for previously searched labels, an incremental approach is used where in each iteration t we look for the set of labels that improves the current solution s.t. y t i = y 0 i + t δ lit , which is a temporal minimization problem. Then (22) can be re-written as a labeling problem:

E t (L ∆ , L Ω ) = V f NW ({y t-1 1 , l ∆ 1 , . . . , y t-1 d , l ∆ d }), I) + α ij i∈G j∈N (i) V ij (y t-1 i , y t-1 j , l ∆ i , l ∆ j ) + β c c∈C V c (w t-1 c , l Ω c ). (24)
We solve the minimization of the higher-order cliques in (24) by transforming them into quadratic functions [START_REF] Rother | Minimizing sparse higher order energy functions of discrete variables[END_REF] using a (t + 1)-state switching variable which finds the deviation function which assigns the lowest cost to the labeling:

min V c (w 0 c ) = min w 0 c ,z∈{1,2,...,t+1} f (z) + i∈c g(z, w i ) (25)
where f (z) = {θ q , θ max } is a cost assigning function depending on the state variable z and g(z, w i ) = ψ q il when z = q and w i = l ∈ L Ω , while g(z, w i ) = 0 when z = t + 1. We apply a Primal-Dual algorithm called FastPD [START_REF] Komodakis | Performance vs computational efficiency for optimizing single and dynamic MRFs: Setting the state of the art with primal-dual strategies[END_REF] which can efficiently solve the inference problem in a discrete domain by formulating the duality theory in linear programming. The advantage of such an approach lies in its generality, efficient computational speed, and guarantees the global optimum without the condition of linearity.

Experimental Validation

In order to validate the performance of the method, we considered modeling pathological spinal columns represented as articulated shapes for computerized tomography (CT) inference. We used a database of 711 spine models reconstructed in 3D from X-ray images, exhibiting different types of deformations relative to global and local shape geometries. For each spine, 6 landmarks on each of the 17 vertebrae composing the spinal column where used to define the LCS in order to extrapolate the inter-object transformation as described in Section 2.1. Models were composed between 3831 and 6942 vertices depending on the vertebra level. Dimensionality reduction provides the embedding manifold shown in [Fig. 2] with local variation mode from a selected patch. Optimal neighborhood size was found at K = 10, while intrinsic dimensions was d = 7 and n = 5 as shown in [Fig. 3(a)] dictating the number of nodes in our global graph model. We tested the algorithm on a subset of 20 unseen cases from the database. Successful examples of the obtained results from sparse sample points and CT inferred data are shown in [Fig. 4]. Additional results are presented in [Fig. 5]. Towards checking the robustness of the method, we evaluated the performance based on the density of the input sample points. [Fig. 6(a)] demonstrates the model's accuracy with respect to the number of available points, visibly affecting the global shape when density is reduced to 20%. To further validate the robustness, we added gaussian distributed noise to the target data. While distortions begin to appear when σ = 3, the inferences obtained in [Fig. 6(b)] are promising by simply varying the prior weight complements α and β. The method also performs well when local objects are completed occluded or missing [Fig. 6(c)]. We quantitatively compared our method to an AAM modeling based on global PCA [Fig. 3(b)]. Dice scores and root-mean-square (RMS) landmark distances computed from the test cases show improvement of the proposed MRF approach via a non-linear shape analysis. To demonstrate the potential of RT n°0376 our method for other modeling applications, we produced results on deformable arm poses (GRAIL, Univ. Washington) via an articulated representation shown in [Fig. 7]. These validations prove how our method elegantly encodes prior knowledge with image constraints in an MRF framework, and efficiently minimizes the energy term to converge towards an optimal solution. One drawback remains the computational time due to the inverse regression mapping and higher-order clique potential minimization.

Experiments were conducted in C++ on a 2.8GHz Intel P4 processor and 4GB memory.

Discussion

We have proposed a novel method for inferring articulated shape models. Our main contribution consists in modeling complex, non-linear patterns of prior deformations in a Riemannian manifold embedding. New point-based models are created from statistical knowledge in terms of global and local variations. To this end, we introduced a conditional regression kernel to perform the inverse mapping to the ambient space from neighbors selected with an articulated distance metric based on intrinsic and orientation properties. One observation of the non-linear embedding is that it avoids creating shape distortions and eliminates the need to solve large-dynamic programming problems, thus saving computational time and memory space. The proposed framework based on higher-order MRFs is general, and can be extended for other applications in vision and medical imaging to accommodate for pose modeling. Real-time inference of articulated models, motion tracking and medical guidance based on optimal higherorder clique decomposition would be beneficial towards this end.

RT n°0376

INRIA

Figure 1 :

 1 Figure 1: Schematic representation of inter-object transformations.

Figure 2 :

 2 Figure 2: Manifold embedding result of the spine dataset with illustration of a local patch M k for individual shape warping of two first modes n = 2 with variation of ±3SD.

Figure 3 :

 3 Figure 3: (a) Dimensionality of the global manifold based on residual variance (top) and of local shape eigenvalues using the % of covered shape variance (bottom);(b) Boxplots of dice overlap coefficients and landmark errors comparing our method to AAM.

Figure 4 :

 4 Figure 4: Successful spine inference results of triangulated mesh ADMs to sparse (red points) and CT image with orthogonal views.

Figure 5 :

 5 Figure 5: (a)-(b) Examples of successful and erroneous inference results in various configurations.

Figure 6 :

 6 Figure 6: Triangulated mesh models represent inferred ADM; red points illustrate the target image data. (a) Decreasing number of available points to infer the ADM affecting global shape coherence. (b) Local shape distortions with significant noise level increase σ added to target points (error-coded models for ground-truth distances). (c) Robustness of method towards missing target components and occluded parts.

Figure 7 :

 7 Figure 7: Human body part deformation results for the arm.

RT n°0376

INRIA


Alexander Heußner

Jérôme Leroux

Anca Muscholl

Grégoire Sutre

Reachability Analysis of Communicating Pushdown Systems

The reachability analysis of recursive programs that communicate asynchronously over reliable Fifo channels calls for restrictions to ensure decidability. We extend here a model proposed by La Torre, Madhusudan and Parlato [LMP08], based on communicating pushdown systems that can dequeue with empty stack only. Our extension adds the dual modality, which allows to dequeue with non-empty stack, and thus models interrupts for working threads. We study (possibly cyclic) network architectures under a semantic assumption on communication that ensures the decidability of reachability for finite state systems. Subsequently, we determine precisely how pushdowns can be added to this setting while preserving the decidability; in the positive case we obtain exponential time as the exact complexity bound of reachability. A second result is a generalization of the doubly exponential time algorithm of [LMP08] for bounded context analysis to our symmetric queueing policy. We provide here a direct and simpler algorithm.

Introduction

The verification of safety properties for distributed programs, e.g., client/server environments, peer-to-peer networks, or Grid applications, relies on the decidability of the reachability problem. In this paper we reconsider recursive queueing concurrent processes (Rqcp), one possible model for such systems which was studied recently by La Torre, Madhusudan, and Parlato [START_REF] La Torre | Context-bounded analysis of concurrent queue systems[END_REF]. It is a natural idea to combine peer-to-peer asynchronous communication (via point-to-point, unbounded, reliable Fifo channels) with some automaton-based model for individual peers (e.g., pushdown automata or Petri nets). We call such combined models queueing concurrent processes (Qcp). Since communicating finite-state automata are the most elementary instantiation of Qcp, reachability is in general undecidable [START_REF] Brand | On communicating finite-state machines[END_REF]. Furthermore, adding recursion (i.e., replacing finite-state by pushdown automata) yields an additional source of undecidability. One of the main motivations in this paper is to separate these two sources of undecidability: we consider behavioral restrictions for which reachability for communicating finite-state machines is decidable, and then look under which conditions we can add recursion to the model. The challenging task is to derive conditions that conserve the simplicity and expressiveness of the model.

In general, there are three main directions to cope with the undecidability of communicating finite-state machines: restricting the communication architecture, assuming that channels are lossy, or adding semantic restrictions, e.g., that sending/receiving of messages can be scheduled in such a way that runs can be executed with channels of bounded size. In this paper we stick to the latter approach, as described in more detail below. Our point of departure is the work of La Torre et al. [START_REF] La Torre | Context-bounded analysis of concurrent queue systems[END_REF], which introduced Rqcp together with a behavioral restriction on the combined use of channels and pushdowns. Informally, a Rqcp is well-queueing if pushdown processes can only dequeue (read) messages when the stack is empty (they can enqueue messages without restriction). Wellqueueing expresses an event-based programming paradigm: tasks are executed by threads without interrupt, i.e., a thread accepts the next task only after it finished the current one. One of the results of [START_REF] La Torre | Context-bounded analysis of concurrent queue systems[END_REF] is that Rqcp have a decidable reachability problem if and only if their communication architecture is a directed forest; in the decidable case, the latter paper provides a doubly exponential upper bound by a reduction to bounded-phase multi-stack pushdown systems [START_REF] La Torre | A robust class of contextsensitive languages[END_REF].

Our contribution. We extend the work of La Torre et al. [START_REF] La Torre | Context-bounded analysis of concurrent queue systems[END_REF] in several directions. First, we add a dual notion to well-queueing: a pushdown process can enqueue (send) messages only with empty stack (but can dequeue messages without restriction). Oriented communication architectures, as presented here, combine these two notions, by fixing the behaviour of the two endpoints of each channel. This dual notion to well-queueing arises naturally if one wants to model interrupts: a server might need to accept tasks from high priority clients independently of the status of the running task.

Second, we exhibit a precise characterization of those oriented architectures for which the Rqcp model has a decidable reachability problem over so-called eager runs. Informally, a run is eager if the sending of a message is immediately followed by its reception, a notion closely related to existentially 1-bounded communication [START_REF] Lohrey | Bounded MSC communication[END_REF]. Eagerly communicating finite-state machines are a well-studied model, enjoying good expressiveness and decidability properties [START_REF] Genest | A Kleene theorem and model checking algorithms for existentially bounded communicating automata[END_REF]. Here, we use eager runs in order to rule out undecidability stemming from unbounded channels. We show that reachability of Rqcp over eager runs is ExpTimecomplete in the decidable case. This result generalizes and improves the doubly exponential time decision procedure of [START_REF] La Torre | Context-bounded analysis of concurrent queue systems[END_REF], which holds for architectures without undirected cycles (polyforest architectures).

Eagerness is a relatively strong requirement, hence, we show how it arises rather naturally, by imposing a semantic restriction on the communication flow: the mutex restriction demands that in every reachable configuration there is no more than one non-empty channel per cycle. In particular, Qcp over polyforest architectures are mutex. Actually mutex can be seen as a generalization of the half-duplex restriction studied in [START_REF] Cécé | Verification of programs with half-duplex communication[END_REF].

La Torre et al. [START_REF] La Torre | Context-bounded analysis of concurrent queue systems[END_REF] propose a second approach to solve the reachability problem for Rqcp, inspired by recent work on reachability with bounded contexts in the verification of concurrent programs [START_REF] Qadeer | Context-bounded model checking of concurrent software[END_REF]. They show that bounded-context reachability for well-queueing Rqcp is decidable in time doubly exponential in the number of contexts. Again, this result is obtained by a reduction to bounded-phase multi-stack pushdown systems [START_REF] La Torre | A robust class of contextsensitive languages[END_REF]. Our second main contribution is to extend the bounded-context result to Rqcp that allow for the two dual notions of well-queueing. Moreover, our algorithm is direct and simpler than the one involving bounded-phase multi-stack pushdown systems.

A long version of this paper that includes all proofs omitted due to space limitations can be found at http://hal.archives-ouvertes.fr/hal-00443529/.

Related work. In the context of thread programming, other notions of synchronization between pushdowns arise naturally. Earlier publications considered synchronization via shared memory, such as local/global memory in [START_REF] Bouajjani | Reachability analysis of multithreaded software with asynchronous communication[END_REF][START_REF] Bouajjani | Regular symbolic analysis of dynamic networks of pushdown systems[END_REF] or bags in [START_REF] Sen | Model checking multithreaded programs with synchronous atomic methods[END_REF][START_REF] Jhala | Interprocedural analysis of asynchronous programs[END_REF]. The paper [START_REF] Bouajjani | Reachability analysis of multithreaded software with asynchronous communication[END_REF] showed that boundedcontext reachability can be solved in exponential time, whereas [START_REF] Sen | Model checking multithreaded programs with synchronous atomic methods[END_REF] provided an exponential space lower bound for reachability (without context bounds). More recently, synchronization in the form of state observation was considered in [START_REF] Atig | On the reachability analysis of acyclic networks of pushdown systems[END_REF]. The latter model was shown to be decidable only for acyclic architectures, and is strongly related to lossy systems [START_REF] Abdulla | Verifying programs with unreliable channels[END_REF][START_REF] Finkel | Well-structured transition systems everywhere![END_REF].

Queueing Concurrent Processes

We write X • ∪X ′ for the disjoint union of X and X ′ . Let Σ denote an alphabet (i.e., a finite set of letters). We write Σ * for the set of all finite words (words for short) over Σ, and we let ε denote the empty word. Moreover, we use standard complexity classes such as polynomial space (PSpace), deterministic exponential time (ExpTime), and doubly exponential time (2ExpTime). For more detailed definitions the reader is referred to textbooks like [START_REF] Papadimitriou | Computational Complexity[END_REF].

A communication architecture T (or architecture for short) is a pair P, Ch with a finite non-empty set P of processes and a finite set of point-to-point channels Ch ⊆ (P × P) \ id P .

Remark 1.1. Our definition of architecture forbids self-loops, as well as two distinct channels in the same direction between a pair of processes (without further restriction, these settings are immediately Turing equivalent).

Definition 1.2. A system of queueing concurrent processes (Qcp) over a given architecture T = P, Ch is a tuple A = (S p ) p∈P , (Σ p ) p∈P , (∆ p ) p∈P , (s 0 p ) p∈P , M with M a finite message alphabet. For each process p ∈ P, the tuple S p , Σ p , ∆ p , s

Σ com p = {p!q(m) | (p, q) ∈ Ch and m ∈ M } ∪ {p?q(m) | (q, p) ∈ Ch and m ∈ M }.
Local transitions are given by the rules in ∆ p ⊆ S p × Σ p × S p , and the initial state of process p is s 0 p . The global state space is S = p∈P S p and the global initial state is s 0 = (s 0 p ) p∈P ∈ S. By Σ = • p∈P Σ p we denote the set of all possible actions in A.

The size of A is p∈P |∆ p |.
As usual, p!q(m) denotes the send of message m from process p to process q, whereas q?p(m) denotes the matching receive on process q.

Note also that the S p and hence S need not necessarily be finite. If S is finite, we will call A a finite Qcp (or communicating finite-state machine, Cfm). The local transition systems given by S p , Σ p , ∆ p , s 0 p could be, for example, finite automata, counter automata (including Petri nets), pushdown automata. As usual, we define the semantics of A as labeled infinite-state transition system:

Definition 1.3. A Qcp A represents an LTS A = C, Σ, →, c 0 with config- urations C = S × (M * )
Ch and the initial configuration c 0 = s 0 , (ε, . . . , ε) , i.e., all channels are initially empty. We write a configuration as c = s, w where s = (s p ) p∈P is the global state and w = (w p,q ) (p,q)∈Ch are the channel contents. Further, for any p ∈ P and a ∈ Σ p , s, w a -→ s ′ , w ′ is a transition in C ×Σ ×C with s ′ = (s ′ p ) p∈P , w ′ = (w ′ p,q ) (p,q)∈Ch , if (s p , a, s ′ p ) ∈ ∆ p and the following holds: (i) s q = s ′ q for all q = p, (ii) if a ∈ Σ loc p then w = w ′ , (iii) if a ∈ Σ com p with a = p!q(m) then w ′ p,q = w p,q m and w ′ s,t = w s,t for (s, t) ∈ Ch\{(p, q)}, (iv) if a ∈ Σ com p with a = p?q(m) then w q,p = mw ′ q,p and w ′ s,t = w s,t for (s, t) ∈ Ch\{(q, p)}.

A finite run ρ in the labeled transition system A from a configuration c 0 to

c n ∈ C is a sequence c 0 , a 1 , c 1 , a 2 , c 2 , • • • , a n , c n where c i-1 ai -→ c i for 1 ≤ i ≤ n.
The length of ρ is n, and a run of length 0 is defined as ρ = c 0 .

A configuration c ∈ C is reachable in the Qcp A if there exists a finite run ρ = c 0 , a 1 , c 1 , • • • , c n starting in the initial configuration c 0 = c 0 and ending in c n = c. We define the reachability set as Reach A = {c ∈ C | c is reachable in A}. The reachability question asks for a given A and a configuration c ∈ C whether c ∈ Reach A . Given a state s ∈ S, the control state reachability question asks whether one can reach a configuration with (control) state component s regardless of the channel content, i.e., whether {s} × (M * ) Ch ∩ Reach A = ∅. Both questions are undecidable for finite Qcp with at least two processes that are connected by two channels [START_REF] Brand | On communicating finite-state machines[END_REF].

The trace of a run

ρ = c 0 , a 1 , c 1 , a 2 , c 2 , • • • , a n , c n is the sequence of actions tr (ρ) = a 1 • • • a n ∈ Σ * .
Since channels are Fifo, we can speak about matching send/receive pairs: a i , a j form such a pair if (1) a i = p!q(m), a j = q?p(m), and

(2) |{ℓ | ℓ ≤ i, a ℓ = p!q(n), n ∈ M }| = |{ℓ | ℓ ≤ j, a ℓ = q?p(n), n ∈ M }|.
We call two runs ρ, ρ ′ order-equivalent if they can be transformed one into the other by iteratively commuting adjacent transitions labeled by a and b, resp., such that (i) a, b do not belong to the same process, and (ii) a, b are not a matching send/receive pair.

Lemma 1.4. If ρ, ρ ′ are order-equivalent runs of A starting in the same configuration, then ρ, ρ ′ end in configurations with the same control state. Definition 1.5. A run ρ with trace tr (ρ) = a 1 • • • a n is eager if the following holds: if a i = p!q(m) for some 0 ≤ i < n then either a i+1 = q?p(m) or no action a j with j > i is a receive on (p, q). A channel that does not permit further receives is in its "growing phase".

A Qcp A is eager if each c ∈ Reach A is reachable by some eager run. Eager runs, modulo the fact that Definition 1.5 allows for runs which end in a sequence of (unmatched) send actions, are closely related to globally 1-bounded runs, whereas eager Qcp are close to existentially globally 1-bounded Cfm [START_REF] Lohrey | Bounded MSC communication[END_REF][START_REF] Genest | On communicating automata with bounded channels[END_REF]. The (control-state) reachability question on eager runs asks whether one can reach a given (control state) configuration by some eager run.

Recursive QCPs and oriented communication architectures. In the following we introduce Rqcp together with a symmetric version of the "well-queueing" property from [START_REF] La Torre | Context-bounded analysis of concurrent queue systems[END_REF]. Informally speaking, Rqcp are Qcp where all basic processes are pushdown automata.

A recursive Qcp (Rqcp) is a Qcp given by (S p ) p∈P , Σ, (∆ p ) p∈P , s 0 , M, Γ where each process p is a pushdown process over the local states S p ⊆ Z p × Γ * with a finite set Z p of control states and the content of the pushdown stack represented by a word over the stack alphabet Γ ; further, the local actions Σ loc p contain push(γ) and pop(γ) for each γ ∈ Γ , and we assume that in the initial state all stacks are empty, i.e., s 0 ∈ p∈P (Z p × {ε}).

A well-queueing Rqcp in [START_REF] La Torre | Context-bounded analysis of concurrent queue systems[END_REF] is one where a process can only receive when its stack is empty. Here, we dualize this concept by also allowing channels where the sender (but not the receiver) must have an empty stack.

Definition 1.6. An architecture T = P, Ch together with a labeling of the channels Ch = Ch s ∪ Ch r as "send restricted" (Ch s ) and / or "receive restricted" (Ch r ), is called oriented.

For pushdown networks the previous definition translates, informally speaking, as follows: a process p can send on a channel (p, q) ∈ Ch s only with empty stack. Symmetrically, p can receive on a channel (q, p) ∈ Ch r only with empty stack. By definition, channels are restricted at least at one end.

The semantics of an Rqcp R is given by a labeled transition system A = C, Σ, →, c 0 analogously to Definition 1.3 except for transitions (s, w)

a -→(s ′ , w ′ ) that correspond to a (local) pushdown transition (s p , a, s ′ p ) ∈ ∆ p : (i)' if a ∈ Σ loc
p , then w = w ′ and further push and pop behave as local pushdown actions;

(ii)' for a = p!q(m) ∈ Σ com p , we demand additionally to (ii) that if (p, q) ∈ Ch s then s p ∈ Z p × {ε}; (iii)' for a = p?q(m) ∈ Σ com p , we demand additionally to (iii) that if (q, p) ∈ Ch r then s p ∈ Z p × {ε}.
Given an oriented architecture T = P, Ch we will use the following notation, that forgets about the direction of the channels and focuses on the (un)limited use of pushdowns: for two processes p, q ∈ P we write p

•-• q if (p, q) ∈ Ch \ Ch s or (q, p) ∈ Ch \ Ch r . Moreover, we write p •-• q if (p, q) ∈ Ch s ∩ Ch r or (q, p) ∈ Ch s ∩ Ch r .
Informally, p •-• q means that for at least one channel between p and q, process p can use its stack without restriction. Similarly, p •-• q means that neither p nor q can use their stacks when communicating. Finally, p •-• q is equivalent to q •-• p. As this notation refers implicitly to a given channel between p and q, we might have both p •-• q as well as p •-• q (or p •-• q) -since both channels (p, q) and (q, p) may exist.

Remark 1.7. A channel can be both send and receive restricted, but we exclude -per definition -channels that are unrestricted at both ends, as this leads immediately to undecidability: one can reduce right away the intersection of two context-free languages to the reachability question on a topology with a single channel that is unrestricted at both ends (and eager runs suffice).

Decidable Oriented Architectures

Several factors lead to the undecidability of the (control-state) reachability question for Rqcp. Particularly, the model is already undecidable even without pushdowns. Our actual motivation in this section is therefore to separate the undecidability which stems from unbounded Fifo queues from the undecidability originating from an unrestricted usage of pushdowns. Hence, we consider a restricted version of the control state reachability question, namely the one on eager runs. In the next section, we show how eager Qcp naturally arise from some natural (and decidable) restrictions on cyclic communication. The most simple example of an eager Qcp is an Rqcp on a polyforest architecture.

Definition 2.1. An oriented architecture T = P, Ch is called confluent if there exist distinct processes p = r 0 , r 1 , . . . , r k , r k+1 = q (k ≥ 1) in T , satisfying the following conditions: (1)

(r i , r i+1 ) ∈ Ch ∪ Ch -1 for all 0 ≤ i ≤ k, and (2) p •-• r 1 and r k •-• q. p r1 r2 r3 . . . rk q Fig. 1. Example of a confluent architecture (mixing •-• and → = Ch notation).
Theorem 2.2. An oriented architecture T = P, Ch admits a decidable Rqcp control-state reachability problem on eager runs if and only if it is non confluent. Moreover, the problem is ExpTime-complete in the latter case.

The only-if direction of the theorem above is not difficult to show. The ifdirection is based on two main ingredients: first, we show how to reorder runs such that we can identify subruns on subarchitectures that start and end with empty stacks; second, we use induction on subarchitectures. The subsequent lemma is the core of the remaining proof of Theorem 2.2. Lemma 2.3. Let R be an Rqcp over a non-confluent architecture T , and consider an eager run ρ of R starting with all stacks empty. Further assume that ρ = ρ 1 a ρ 2 b ρ 3 with a, b ∈ Σ such that for some process p ∈ P:

(i) the stack of p is continuously empty during both subruns ρ 1 and ρ 3 , and continuously non-empty during ρ 2 , respectively; (ii) there is no eager run ρ ′ that is order-equivalent to ρ and has the form

ρ ′ = ρ 1 c ρ ′ 2 d ρ 3 where c, d ∈ Σ with c = a or d = b.
Then, the stacks of all processes occurring in ρ 2 are empty after both ρ 1 and ρ 2 .

Proof. By (ii), each process q occurring in ρ 2 has a "proof" for its presence in ρ 2 . This proof consists of a simple, unoriented path in T between p and q.

We explain this more formally: notice first that a is necessarily a push action and b its matching pop action on p. Consider now the first action of some q occurring in ρ 2 : this is a communication either with p, or with some different process r. Thus, we can assume inductively that there is a simple unoriented path between p and r, which can be extended to a path between p and q.

Since p's stack is continuously nonempty during ρ 2 , we know that the communication between p and the first node p ′ on the path above is on a channel of type p •-• p ′ . Due to the non-confluent property, we cannot have r •-• q; hence, q's stack must be empty at the beginning of ρ 2 . A symmetric argument applies at the end of ρ 2 .

Proposition 2.4. The control-state reachability problem for Rqcp on eager runs over non-confluent architectures is ExpTime complete.

Proof. For the upper bound we show how to compute inductively in ExpTime all pairs of global states (s, s ′ ) ∈ S 2 such that there is an eager run over T = P, Ch , starting in s with stacks and queues empty, as well as ending in s ′ with possibly empty stacks and queues. Actually, we need to compute more, namely for which sets P ⊆ P of connected processes we can reach s ′ from s (with empty stacks and queues). W.l.o.g. we assume that T is connected. We (arbitrarily) order the processes in P and choose the first process p ∈ P that has at least one edge (channel) of type p •-• * in T and at least one of type either p •-• * or p •-• * (with * denoting an arbitrary process in P \ {p}). Since the architecture is non confluent, only two cases can occur if such a process does not exist: either (i) T contains no •-• edge at all, or (ii) T = P = r, r 1 , . . . , r k , Ch , where:

a channel of type •-• is one of (r, r i ), for some i, -a channel of type •-• is one of (r i , r j ), for some i, j.

r

r 1 r k r i Fig. 2. Case (ii).
Let us first consider case (i) where T has no channel of type •-•, i.e., all channels in T are of type •-•. In this case we can reorder any eager run into an orderequivalent eager run where messages alternate with local pushdown runs, each of them starting and ending with empty stack. This amounts to solving the controlstate reachability problem for one pushdown automaton of size exponential in |P|, which is possible in ExpTime.

We consider now case (ii). Here, we may assume w.l.o.g. that r •-• r i for every 1 ≤ i ≤ k (just add channels of type r •-• r i if there is no channel between r and r i ). Assume that there exists an eager run from s to s ′ , starting with all stacks empty. For simplicity we also assume that the run ends with all channels empty (the special case of unmatched sends can be handled similarly). This run can be reordered into an order-equivalent eager run ρ of the following form:

ρ = σ 0 α 0 m 1 σ 1 α 1 m 2 σ 2 α 2 • • • m n σ n α n
(1) where process r does not occur in

σ 0 • • • σ n ,
each m i is a message (i.e., send/receive pair) between r and some r j , -each α i is a sequence of local (pushdown) actions of r, -each σ i (i < n) starts and ends with all stacks of processes r j empty.

We obtain the previous reordering by scheduling the messages between r and the r j as late as possible. That is, the run starts first with actions not involving r (subrun σ 0 ), plus some pushdown actions of r (subrun α 0 ). Then the first message m 1 between r and some r j follows. All remaining actions are in the future of this message, and σ 1 is a run over the r j (which may synchronize among themselves by communication), etc.

We check that the stacks of the r j are always empty at the beginning/end of each σ i by considering the subruns σ 0 α 0 m 1 σ 1 , σ 0 α 0 m 1 σ 1 m 2 α 2 σ 2 , etc.: consider the first occurrence of some process q occurring e.g. in σ 1 . Either this occurrence is the message m 1 between q and r (so r •-• q), or it is a synchronization with some r i , so r i •-• q. In particular, q's stack must be empty at the end of σ 0 .

The existence of an eager run as above can be checked by first pre-computing the control-state reachability for T \ r, which corresponds to the first case previously considered in this proof (recall that we compute summaries, i.e., all pairs of global states that can be reached starting/ending with empty stacks and queues). Then the question for T reduces to the control-state reachability of a pushdown automaton (that of process p) of size exponential in |P|, thus showing the claim.

We now get back to the situation where T contains some process p ∈ P with at least one channel of type p •-• * in T , and at least one channel of type either p

•-• * , or p •-• * in T .
Consider an (eager) run over T , starting (and possibly ending) with all stacks empty. Again we assume, for convenience, that all channels are empty at the end. The run can be reordered such that we obtain an order-equivalent run of the form ρ 0 σ 0 ρ 1 • • • σ n-1 ρ n , where:

every subrun σ i (i < n) starts and ends with empty stacks for all processes q = p occurring in σ i , -p's stack is continuously empty during ρ i , and continuously non-empty during σ i , for each i.

We need to explain why we may assume that the subruns σ i start / end with empty stacks for each q = p. The reason is that we schedule the internal push / pop actions of p that start / end a phase with non-empty p-stack such that push actions have lowest priority, and pop ones have the highest one. Such push / pop pairs on p delimit the subruns σ i and Lemma 2.3 can be applied to ρ 0 σ 0 ρ 1 , (ρ 0 σ 0 ρ 1 )σ 1 ρ 2 , etc. The existence of suitable subruns σ i can be checked inductively: notice that channels of type p •-• * and p •-• * can be removed from T , since p's pushdown is non-empty during every σ i , hence such channels are not used. More formally, we check for each (connected) subset P ⊆ P with p ∈ P , and each pair of starting/ending states s = (s q ) q∈P , s ′ = (s ′ q ) q∈P whether there is a run of processes from P from s to s ′ in the modified architecture, starting/ending with empty stacks and queues.

Consider now a set P corresponding to processes occurring in some run σ i . Notice first that every process in P can be reached from p via messages involving only P ; symmetrically, every process in P can reach p via messages involving only P . For each such set P we can introduce new synchronization messages between processes in P such that we replace σ i by a sequence S i of new messages with the following properties: the first message in S i is sent by p, the last message is received by p, and every process in P occurs among the receivers in S i . Such a sequence S i can be used to enforce that processes in P i go (in a sort of metatransition) from state s = (s q ) q∈P to state s ′ = (s ′ q ) q∈P , thus replacing σ i . We can enforce that the new messages occur only in form of sequences S i , by encoding S i with message contents. In order to avoid an exponential blow-up in the size of the Rqcp we record the possible sequences S i separately. Notice that using these sequences in the base step of the induction does not affect the ExpTime upper bound.

We can now apply induction on the modified Rqcp in order to check whether there is some run of the form ρ ′ = ρ 0 S 0 ρ 1 S 1 . . . S n-1 ρ n . The induction is possible since we can transform the channels of type p •-• * into type •-•(since p does not use its pushdown in ρ ′ ).

To summarize, the induction is done on two parameters: either we decrease the overall number of channels, or we change at least one channel of type •-• into type •-•. We first check the existence of runs σ i inductively in exponential time, by computing reachability for every pair of global states and subset of processes (of which there are exponentially many). Then we modify the Rqcp according to the previous calls and check inductively the existence of a single run ρ ′ (again, this is done for each pair of global states and set of processes).

Finally, let us comment on the lower bound: It is known (and probably folklore) that the following problem is ExpTime-complete: checking the emptiness of the intersection of a pushdown with n finite automata. The hardness follows easily by a reduction from linearly bounded alternating Turing machines. Actually, a closely related problem is shown to be ExpTime-hard in [START_REF] Esparza | Model checking LTL with regular valuations for pushdown systems[END_REF], namely the reachability problem for pushdowns with checkpoints. Clearly, the intersec-tion between a pushdown and n finite automata can be simulated on a topology T = P, Ch with P = {r, r 1 , . . . , r n } and r •-• r i for each i.

3 From Mutex QCP to Eager QCP

The previous section showed how to decide the control state reachability for Rqcp (and therewith finite Qcp) on eager runs. Nevertheless, restricting the communication to eager runs seems rather strong at first glance. In the following, we will show how eagerness arises naturally on two practically relevant communication architectures: polytrees and cyclic architectures with mutex restriction. Further, we discuss the reduction of the control state reachability problem for (possibly infinite) mutex Qcp to their underlying local transition systems, like e.g., Petri nets.

Any communication architecture T can be regarded as directed graph P, Ch ; let UCycle(T ) be the set of its undirected simple cycles. A cycle is undirected if we ignore the direction of the channels, and simple if it has no subcycle of smaller length. Definition 3.1. A configuration c of a Qcp A is mutex with respect to a given architecture T if for every cycle α of UCycle(T ) at most one of the channels occurring in α is non-empty in c. A Qcp A is called mutex with respect to a given architecture T if every c ∈ Reach A is mutex.

Before discussing mutex Qcp in detail, we first recall two known results that are subsumed by our definition of mutex: Remark 3.2. A special case of mutex Qcp was considered in [START_REF] La Torre | Context-bounded analysis of concurrent queue systems[END_REF]: polyforest architectures over finite Qcp, as well as (well-queueing) root-to-leaf directed forests for Rqcp. Their decidability proof relied on the idea that, on any tree architecture, we can reorder runs such that first all actions of the root process are scheduled, and then, in breadth-first order, the actions of all others. Consequently, each run could be partitioned into a bounded number of contexts (bounded by |P|) where in each context only one process executes all its actions by reading on one unique incoming channel from its tree parent (and -in the case of Rqcp-solely when its local stack is empty). Hence, the decidability problem reduced to the control state reachability for a bounded-phase multistack pushdown system, which is known to be decidable in doubly exponential time [START_REF] La Torre | A robust class of contextsensitive languages[END_REF].

We will show in the following that mutex Qcp are eager, and, consequently, apply the results of the previous section to obtain the decidability of controlstate reachability via a direct proof. Moreover, recall that the complexity of the algorithm of the previous section is ExpTime, so one exponential less than the positive results of [START_REF] La Torre | Context-bounded analysis of concurrent queue systems[END_REF].

Remark 3.3. Runs over an architecture of two finite processes connected by two channels where each reachable configuration is mutex are known as "half-duplex communication". For these, it is known how to decide the (general) reachability question by computing a recognizable description of the channel contents [START_REF] Cécé | Verification of programs with half-duplex communication[END_REF]. Quasi-stable systems are a semantic ad-hoc extension of this idea to larger, cyclic architectures of finite Qcp [START_REF] Cécé | Verification of programs with half-duplex communication[END_REF], which is subsumed by our mutex condition.

Proposition 3.4. Given a Qcp A that is mutex with respect to a given architecture T , each of its runs has an order-equivalent eager run.

Proof. In the following, we will differentiate the occurrences of one action by referring to them as events. For each process p ∈ P occurring in ρ there is a first, initial event w.r.t. ρ which will be abbreviated f p ; further, each receive event has a preceding matching send in ρ. All events that belong to the same process are totally ordered. Consequently, we define the partial order "before" (denoted by <) between events as the transitive closure of the previous two cases. In the following we will focus only on matched communication events, by considering send actions on channels that already entered their growing phase as internal actions. We will inductively define a reordering for a run ρ whose first configuration c 0 = s, w fulfills w p,q = ε if there exists a receive event q?p in ρ.

Assume we have a run ρ from c to c ′ and c fulfills the previous property. First, we pick an initial send event f p on a channel from process p to q which either (i) has no matching receive in ρ (i.e., it is the first send in a growing phase), or (ii) its matching receive on process q is also initial, hence, equal to f q . In case (i) we schedule f p first and reorder inductively the remaining run starting from c ′′ with c fp -→ c ′′ towards c ′ . For case (ii) we first schedule f p and then f q before we inductively reorder the remaining run starting from c ′′ with c fp -→ c ′′′ fq -→ c ′′ . Note than in both cases c ′′ satisfies our requirement.

Next we have to show that it is always possible to apply cases (i) or (ii) above, to any run ρ of our mutex Qcp. The general idea is as follows. Suppose that we pick an initial send event f p0 on process p 0 that has a matching receive e p0 on process p 1 , but e p0 is not initial. Then we can restart our search for an initial event from p 1 on. If f p1 is a send, then we proceed as for p 0 ; else, if f p1 is a receive, we continue with its matching send on process p 2 . As we only have finitely many processes, an unsuccessful, repeated search leads to a cycle in UCycle(T ): p 0 , p 1 , p 2 , . . . , p i , p i+1 , . . . , p k , p k+1 with p k+1 = p i and all p i (i ≤ k) pairwise different. Moreover, we show the existence of at least two non-empty channels on this cycle.

In the following, we slightly abuse notation by writing f i and instead of f pi for the initial event of process p i . We focus on the initial events f i , . . . , f k and their matching events e i , . . . , e k on processes p i+1 , . . . , p k , p k+1 = p i . Obviously, f j+1 < e j for all i ≤ j ≤ k, since both f j+1 , e j occur on process p j+1 and f j+1 is initial. We distinguish the following cases:

(a) all initial events f i , . . . , f k are receives (cf. Fig. 3(a)), then e j+1 < f j+1 < e j for all i ≤ j ≤ k; hence, we arrive at the contradiction e i < e i ; (b) there are at least two sends among the initial actions, for example f j and

f l with i ≤ j < l ≤ k; consequently, c fj -→ c ′′ f l -→ c ′′′ leads
to a configuration which is not mutex (cf. Fig. 3(b) ) and, hence another contradiction; (c) there is only one send event among the initial events of the cycle, say f i .

Then, f i = f k+1 is before e k , and e k is a send event, too (since all f j with i < j ≤ k are receives). It is easy to see that

e k < f k < e k-1 < • • • < f i+1 < e i .
In particular, all events on each of p i+1 , . . . , p k are after e k .

Consider now an event e with f i < e < e k . Notice that e cannot belong to any of p i+1 , . . . , p k , as all events on these processes must take place after e k (cf. Fig. 3(c) for e on a process p that does not participate in the cycle); consequently, the configuration obtained after executing all events before e k is not mutex, as the channels (p i , p i+1 ) and (p i , p k ) are both non-empty.

Corollary 3.5. If a Qcp A is mutex with respect to a given architecture T then A is eager.

Proposition 3.6. The control state reachability for finite Qcp that are mutex with respect to the given architecture is Pspace-complete.

Remark 3.7. Control-state reachability is decidable for particular infinite-state mutex Qcp. For example, if all local transition systems are Petri nets, then the control state reachability question reduces to a Petri net reachability question which is known to be decidable [START_REF] Mayr | An algorithm for the general Petri net reachability problem[END_REF][START_REF] Kosaraju | Decidability of reachability in vector addition systems[END_REF].

The mutex property can be checked effectively for Qcp.

Proposition 3.8. It is PSpace-complete to check whether a finite Qcp is mutex with respect to a given architecture.

Bounded Phase Reachability

Besides their proven practical relevance in the verification of concurrent programs [START_REF] Qadeer | Context-bounded model checking of concurrent software[END_REF], bounded-context reachability allows to attack the (control-state) reachability problem on Qcp from a different angle. In this section, we neither restrict the communication architecture, nor constrain the runs to be eager (or mutex). The price we pay is a (strong) restriction on the form of the possible runs, by fixing the number of contexts. We present in this section a construction that subsumes the 2ExpTime algorithm for bounded-context reachability for well-queueing Rqcp described in [START_REF] La Torre | Context-bounded analysis of concurrent queue systems[END_REF]. Recall that the latter algorithm is based on a reduction to bounded-phase reachability for multi-stack systems. In contrast, our construction below is direct and simpler. We define a phase in an Rqcp run ρ over an oriented architecture as a (contiguous) subrun of ρ consisting of actions on a unique process, say p, that are subject to one of the two restrictions below (defining M -phases and N -phases, resp.):

(M ) Receives are from a unique process, say q, with (q, p) ∈ Ch of type p •-• q or p •-• q. Sends go to (arbitrarily many) processes r with (p, r) ∈ Ch of type p •-• r. (N ) Sends go to a unique process, say q, with (p, q) ∈ Ch of type p •-• q or p •-• q. Receives come from (arbitrarily many) processes r with (r, p) ∈ Ch of type p •-• r.

Notice that M -phases are precisely the phases (contexts) used in [START_REF] La Torre | Context-bounded analysis of concurrent queue systems[END_REF], whereas N -phases represent the dual notion. We will refer below to the channel (q, p) (M -phase) resp. (p, q) (N -phase) as the special channel of the phase. A run ρ of an Rqcp is K-bounded, if we can write ρ = ρ 1 • • • ρ K , with each ρ i a phase as above.

Theorem 4.1. Given an Rqcp A and an integer K, the K-bounded control state reachability problem for A can be solved in time doubly exponential in the number K of phases (but polynomial in the size of A).

Sketch of proof.

The basic idea is to decrease the number of phases in a particular order: M -phases are deleted for right to left (a sort of pre-computation), whereas N -phases are deleted from left to right (post-computation). Deleting a phase i belonging to some process p amounts to synchronizing a finite automaton obtained from A p and phase i with the current automaton A q of the process communicating with p on the special channel of phase i. We obtain this finite automaton by exploiting the fact that p's stack is empty while communicating in phase i on the special channel. In addition we must ensure that the phase i that we delete starts and ends with empty p-stack. Finally, for a single phase we need to solve a reachability problem for a single pushdown with doubly exponentially many states. The details can be found in the long version of this paper.

Remark 4.2. Adapting proof ideas from [LMP07, ABH08], we can show that the complexity bound in Theorem 4.1 is tight.

Conclusion

Applications. Qcp combine an automata-based local process model with pointto-point communication, which results in an intuitive and simple framework. Since we subsume well-queueing Rqcp, we also inherit their application domains, e.g., event-based programs. The dual restriction to well-queueing (i.e., that sending on a channel is only possible if the stack is empty) covers e.g. "interrupt based" programming models, i.e., threads that can receive messages while still in recursion, as well as extended sensor networks where peers can collect and send data while using their pushdown for computations.

Fig. 4 (b) shows an example for non-confluent architectures that are on the rise with the current focus on Grid computing. The topology depicts a hierarchical overlay network as implemented, for example, in a master-worker protocols. Here, mutual communication is restricted with respect to the hierarchy (in general: •-• top-down and •-• between siblings). Notice also the use of the dual notion to well-queueing, when sending information from lower to higher levels.

Proposition 3.4 allows for further applications, since it does not assume that the Qcp is finite: we can combine locally decidable models for multi-threaded programs (with or without local data), as well as local event-based programs together with eager (or mutex) communication architectures; natural candidates for local models would be Petri Nets, Wsts, or multi-set pushdown systems [START_REF] Sen | Model checking multithreaded programs with synchronous atomic methods[END_REF].

Outlook. We discussed in detail the class of eager Rqcp (as well as mutex Qcp) which both generalize the current lineup of decidable models for asynchronously communicating pushdown systems. Further, we presented an optimal decision procedure for eager Rqcp over non-confluent architectures in ExpTime, as well as a direct and simpler construction for bounded phase reachability for Rqcp.

This paper dealt with the most basic form of verification, namely controlstate reachability. More general reachability questions (w.r.t. configurations) may be interesting to consider. Further decision problems for Qcp, like boundedness or liveness, will be investigated in future work.

As in the previous proof for Proposition 3.6, we have to take care of growing channels by a flag construction.

PSpace-hardness is not hard to obtain by reducing from the non-emptiness test of the intersection of arbitrarily many finite automata. We can construct the finite Qcp in such a way that a "bad" cycle is reached precisely from the target (final) states of the product of the automata.

Proofs for Section: Bounded Phase Reachability Theorem 4.1. Given a Qcp A and an integer K, the K-bounded control state reachability problem for A can be solved in time doubly exponential in the number K of phases (but polynomial in the size of A).

Proof. Let ρ = ρ 1 • • • ρ k be a K-bounded run of A = (S p ) p∈P , M, Σ, (∆ p ) p∈P , s 0 , with each ρ i an M -or an N -phase. W.l.o.g. we suppose that all stacks are empty at the end of ρ. Let (t p ) p∈P ∈ Π p∈P S p be the target control state. We first modify the Rqcp A in such a way that we explicitly refer to phases in action names (actually we only consider a subset of A's behaviors, namely those corresponding to K-bounded runs): we add to each send p!q(m) the number i of the phase in which the send occurs, as well as the number j of the phase in which the matching receive occurs; in the special case j = K + 1, the matching receive does not occur within the first K phases. So p!q(m) is replaced (non-deterministically) by p!q(m | i, j), with i < j ≤ K + 1. We proceed similarly for receives, this time just replacing p?q(m) by p?q(m | i, j) with i < j ≤ K. In addition, we label also the local actions by their phase number. The transition relation ∆ p of each automaton A p = S p , M, Σ p , ∆ p , s 0 p can check that the phase labeling for p's actions is consistent. The set S p of local states increases polynomially in K.

At the beginning of the construction we fix for each phase 1 ≤ j ≤ K: (1) its type, (2) the process to which it belongs, and (3) its special channel. Moreover, we fix the "reading graph" Read , which is a DAG with nodes (p, j), where p is the process of phase j. If j is an M -phase then there is an edge to (p, j) from each (q, i) such that (a) (q, p) is the special channel of phase j and (b) q in phase i sends to p in phase j. Dually, if j is an N -phase then there is an edge from (p, j) to each (q, k) such that (a) (p, q) is the special channel of phase j and (b) p in phase j sends to q in phase k. Clearly there are some natural consistency requirements on Read . Note also that there are only exponentially many graphs Read , so that we can check them one by one.

In addition, we will fix at the beginning some state information for 4 distinguished points in each phase i (belonging to some p): the starting/ending point of the phase, together with the first and last point where p's stack is empty (if any). Of course, some of these points might coincide. For each of the 4 points we fix p's state at that point. Again, there are exponentially many cases to check one by one.

We now show how to reduce the K-bounded control state reachability question for A to a (K -1)-bounded control state reachability question for a modified Rqcp A ′ . Assume first that there is at least one phase of type M , and let j ≤ K be the index of the last M -phase in ρ. Suppose moreover that phase j belongs to process p, and that the special channel of this phase is (q, p) ∈ Ch. We will delete phase j, and adapt the two local automata A p , A q accordingly. First, we show that we may assume w.l.o.g. that ρ j starts and ends with empty p-stack. So let ρ k = αβγ, where β starts and ends with empty p-stack and such that the stack of p is never empty during α, nor during γ (α, β, γ might be empty each). By the definition of M -phases, the only possible actions during α and γ are sends on channels of type p •-• r. Since all phases j < ℓ ≤ K are N -phases, these sends are unmatched within ρ, hence they must be labeled as p!r(c | j, K +1). We now declare all sends p!r(m | j, K +1) as local p-actions, and move α to the end of the last phase i < j that belongs to p. The subrun γ can be either moved to the beginning of the next phase belonging to p, or deleted if there is no such phase. Recall that we know p's state after α and before γ, respectively. So we can adapt the information about the 4 distinguished points in the phases belonging to p that are modified.

So assume that ρ j starts in p-state s p and ends in s ′ p , both with empty pstack. From A p we can define a finite automaton B p that accepts precisely those sequences (m 1 , i 1 ) • • • (m n , i n ) ∈ (M × {1, . . . , j -1}) * such that A p can go from s p to s ′ p on a sequence w of actions of phase j only, (1) starting and ending with empty stack, (2) performing p?q actions only when the stack is empty, and such that (3) the projection of w on p?q actions equals p?q(m 1 | i 1 , j) • • • p?q(m n | i n , j). Note that B p can be constructed on the same state space as A p . We finally modify A p by disallowing actions from phase j.

Finally we change the automaton A q of process q into A ′ q in order to ensure consistency with ρ j . The new automaton A ′ q starts simulating B p in parallel from state s p on the subsequence consisting precisely of symbols q!p(m | i, j) (for all m ∈ M and i < j). Moreover, the simulation of B p must end with state s ′ p . We need now to update the state information for the 4 distinguished points in the phases that are synchronized by B p above. Assume that the incoming edges of (p, j) are (q, i 1 ), . . . , (q, i k ) with i 1 < • • • < i k . So a run of B p "splits" over (a suffix of) phase i 1 , followed by phases i 2 , . . . , i k . So for each of the phases i 1 , . . . , i k and each of the 4 distinguished points (except maybe for some of the points of i 1 ) we need to fix a state of B p . We then run the recursive algorithm for each choice of these states. We will see below what this means for the running time.

Now it remains to deal with the case where there is no M -phase in ρ. Here, it is quite easy to check that we may assume from the beginning that each phase starts and ends with empty stack (of all processes). So here we do not need the extra information about the 4 distinguished points for each phase.

Since there are only N -phases we work from left to right and we eliminate the first phase of the run. Suppose that the first phase belongs to process p. Since it is an N -phase, the only possible actions are sends to a unique process, say q, such that (p, q) ∈ Ch is of type p •-• q or p •-• q. In particular, p's stack is empty during this first phase. So from A p we obtain directly a finite automaton B p that accepts precisely those sequences (m 1 , j 1 ) • • • (m n , j n ) ∈ (M × {2, . . . , k}) * such that A p can go from s 0 p to some s p on a sequence w of actions of phase 1 only, (1) working without stack, (2) performing only p!q actions (or local actions), and such that (3) the projection of w on p!q actions equals p!q(m 1 | 1, j 1 ) • • • p!q(m n | 1, j n ). We also set the initial local state of A p to s p and disable actions from phase 1. We change the automaton A q by synchronizing with B p on q?p(m | 1, j) actions (for all m ∈ M and 1 < j). Moreover the synchronization starts with B p in state s 0 p . Let us now estimate the running time of the construction. For simplicity we focus on the elimination of M -phases. During the phase elimination, automata A q get synchronized with automata B p (with the same state space as A p at the current stage). So local states become tuples of local states of the original automata A r . If the length of the tuples gets at most F (j) after eliminating j < K phases, then we can show that the overall construction is of complexity n O(K 2 (K)) : denote by t(i) the time needed for i phases. Then t(1) ≤ n O(F (K)) and t(i) ≤ (n F (K-i) ) 4K • t(i -1) (the first factor stands for the states of the 4 distinguished points needed in the recursive step). So we have t(K) ≤ n 4K•(F (1)+•••+F (k-1)) t(1) ≤ n O(K 2 •F (K)) as claimed.

Finally, it is not hard to see that F (j) is exponential in j, which gives the doubly exponential upper bound (in K).

Fig. 3 .

 3 Fig. 3. Cycles in the proof of Lemma 3.4 (⋄ : events, arrows : messages).

Fig. 4 .

 4 Fig.4. Non-confluent architectures: (a) ring, and (b) hierarchical master-worker setting -tree-like architecture with •-•-channels between master and workers (distribute tasks and collect results while in computation, send result to own master when computation finished, i.e., stack empty) as well as •-•-channels between workers of the same master; note the •-•-cycle on the top level.

⋆ Work supported by the ANR project AVERISS.

Appendix

Proofs for Section: Decidable Oriented Architectures Proposition . The Rqcp control-state reachability problem is undecidable if the topology is confluent. This holds even if we only consider reachability on eager runs.

Proof. Assume that we have a simple and unoriented path p, r 1 , . . . , r k , q between p = q, with p •-• r 1 and r k •-• q. Then we simulate non-deterministically the intersection of two context-free languages on p, q, by conveying the information about the common input guessed by p, q on the path r 1 , . . . , r k . Since the path between p, q is unoriented, the intermediary processes r i (for 1 ≤ i ≤ k) may have to guess the input, and deadlock if the guess was wrong.

Proofs for Section: From Mutex QCP to Eager QCP Proposition 3.6. The control state reachability for finite Qcp that are mutex with respect to the given architecture is Pspace-complete.

Proof. We want to know for a given control state s = (s p ) p∈P whether a configuration c = (s, w) can be reached, for some w, from the initial configuration c 0 . By Lemma 3.4 we know that c can be reached by a eager run ρ. Consequently, we can guess the run ρ storing the current global state, but without recording more than one message, which gives the (usual) polynomial space bound for products of automata.

Note that we need to take extra care of a channel's growing phase by nondeterministically guessing, at each step of the previous computation, whether a channel will enter its growing phase. This additional information can be stored in polynomial space.

PSpace-hardness follows by reducing from the non-emptiness test of the intersection of arbitrarily many finite automata. Proposition 3.8. It is PSpace-complete to check whether a finite Qcp is mutex (eager, resp.) with respect to a given architecture.

Proof. Assume that A over an architecture T is not mutex and consider a run ρ of minimal length from c 0 to a configuration c = (s, w) that is not mutex. Thus, by minimality all configurations in ρ up to c are mutex.

So by an argument analogous to the construction in the proof of Lemma 3.4, we can reach c ′ -the predecessor of c -by a eager run ρ ′ , and test whether there exists a transition c ′ a -→ c with a send action a that violates the mutex condition for c. Hence, we will guess ρ ′ without storing more than one message at each step, leading to a configuration c ′ in which one channel of a cycle in UCycle(T ) is non-empty and the action a would write on another channel in this cycle.