Laurent Galluccio
email: laurent.galluccio@oca.eu

Olivier Michel
email: olivier.michel@gipsa-lab.inpg.fr

Pierre Comon
email: pcomon@unice.fr

Alfred O Hero

Graph Based K-Means Clustering

Keywords: Data partitioning, unsupervised classification, graph-theoretic methods, minimal spanning trees, similarity measures, information theoretic measures, multispectral imaging

An original approach to cluster multi-component data sets is proposed that includes an estimation of the number of clusters. Using Prim's algorithm to construct a minimal spanning tree (MST) we show that, under the assumption that the vertices are approximately distributed according to a spatial homogeneous Poisson process, the number of clusters can be accurately estimated by thresholding the sequence of edge lengths added to the MST by Prim's alorithm. This sequence, called the Prim trajectory, contains sufficient information to determine both the number of clusters and the approximate locations of the cluster centroids. The estimated number of clusters and cluster centroids are used to initialize the generalized Lloyd algorithm, also known as k-means, which circumvents its well known initialization problems. We evaluate the false positive rate of our cluster detection algorithm, using Poisson approximations in Euclidean spaces. Applications of this method in the multi/hyper-spectral imagery domain to a satellite view of Paris and to an image of Mars are also presented.

Introduction

Motivation and contribution

One of the recurrent problems in pattern recognition, machine learning or data mining, is data clustering [START_REF] Jain | Data clustering: a review[END_REF][START_REF] Xu | Survey of clustering algorithms[END_REF][START_REF] Theodoridis | Pattern Recognition[END_REF]. This is also a major subject of research in the remote sensing community with the emergence of hyper-spectral sensors, which generate a significant amount of data. Clustering consists of partitioning a dataset into groups of points that have high similarity. Similarity in this context means that data belonging to a cluster have similar features.

In this paper a new method is introduced for estimating of the number of clusters present and determining good centroid locations to initialize the k-means algorithm [START_REF] Macqueen | Some methods for classification and analysis of multivariate observations[END_REF], which is widely known to be extremely sensitive to centroid initialization. The presented method relies on properties of the minimal spanning tree (MST) as grown by Prim's algorithm [START_REF] Prim | Shortest connection networks and some generalizations[END_REF], and in particular the Prim trajectory defined as the rate of growth of the MST as Prim's algorithm successively adds edges. The Prim trajectory describes the evolution of the cost of adding edges to the MST as a function of the algorithm iteration index. Data are L dimensional feature vectors sampled from a unknown multimodal distribution. It is assumed that distinct groups or clusters manifest themselves as modes in the density. We show that Prim trajectory is as a one-dimensional unfolded representation of the underlying data probability density function, exhibiting peaks and valleys. Thus, we propose to estimate the number of clusters by thresholding Prim trajectory, in which sharply peaked modes are associated with deep valleys [START_REF] Stuetzle | Estimating the cluster tree of a density by analyzing the minimal spanning tree of a sample[END_REF]. A similar Prim-trajectory clustering approach was used for a bioinformatics application in Olman et al. [START_REF] Olman | Identification of regulatory binding sites using minimum spanning trees[END_REF].

The proposed method can be viewed as an algorithm for initializing an iterative partitioning algorithm like k-means, which is the focus of this paper. k-means is one of the simplest and most popular approaches for solving clustering problems despite its severe sensitivity to initialization. While k-means appears as a final step in the proposed algorithm, other partitioning algorithms could be used. Alternative approaches can be used to identify the number of clusters. Jordan & Ng or Shi & Malik's spectral clustering methods [START_REF] Bach | Learning spectral clustering, with application to speech separation[END_REF][START_REF] Ng | On spectral clustering: Analysis and an algorithm[END_REF][START_REF] Shi | Normalized cuts and image segmentation[END_REF] estimate the number of clusters by thresholding the eigen-spectrum of the graph Laplacian [START_REF] Chung | Spectral Graph Theory[END_REF]. This however remains a difficult task in the presence of outliers, as ordered eigenvalues usually do not exhibit an obvious threshold. As our method involves only those vertices in the graph that are close to the most concentrated regions (modal peaks) of the data this problem is circumvented.

The problem of setting the threshold is tackled by adopting a hypothesis testing framework. Under the null (no clusters) hypothesis we model the data as a realization of a spatial homogeneous Poisson spatial point process over feature space. A Neyman-Pearson hypothesis test procedure is then used to detect the modes.

The construction of the Prim trajectory generally requires the computation of a distance matrix between the data points, which may be impractical for large amounts of data. We propose an algorithm for calculating Prim's trajectory that avoids the difficult computation of huge similarity matrices. The algorithm is based on a data-driven hierarchical classification.

Our Prim-based algorithm can be applied to a variety of distance metrics or similarity measures in feature space. Often a non-symmetric information theoretic similarity measure can be used to obtain improved clustering results. We show that use of such information measures can significantly improve clustering performance on real-world data imaging sets.

Context and brief review of previous works

The general data clustering problem has long been studied and many approaches have been proposed (see [START_REF] Jain | Data clustering: a review[END_REF][START_REF] Xu | Survey of clustering algorithms[END_REF][START_REF] Theodoridis | Pattern Recognition[END_REF] for review). One can divide the existing methods into two classes: hierarchical and partitional clustering algorithms. Hierarchical algorithms evolve the number of clusters, e.g., as in agglomerative or divisive clustering, while partitional algorithms fix the number of clusters, e.g., as in modal or k-means clustering.

Most hierarchical clustering algorithms are based on popular single-link or complete-link algorithms. These methods often suffer from prohibitive computational time due to the need to construct a dendrogram on a large datasets. Their stability in the presence of outliers and their sensitivity to the applied dendogram thresholds are problematic. One of the most popular partitional clustering algorithm is the generalized Lloyd algorithm used to implement k-means [START_REF] Macqueen | Some methods for classification and analysis of multivariate observations[END_REF], which will be the focus of this paper. The present work is developed in the context of unsupervised clustering. Neither the underlying distributions, nor the number and locations of meaningful clusters are known. Hence no prior information is required.

In k-means clustering k points are randomly chosen as the initial cluster centers. The choice of the number k of clusters is usually based on some heuristic. Each data point is assigned to the group that has the closest center. The cluster centers of mass are then recomputed. The assignment and re-computation steps are iterated until the intra-cluster variance converges to a minimum (see section 2.1). This algorithm, as well as the classical ISODATA algorithm [START_REF] Ball | ISODATA, a novel method of data analysis and classification[END_REF], belongs to the class of squared error algorithms which aim at minimizing some objective function in order to cluster the data.

The k-means algorithm clusters the data into non-overlapping convex groups and it always converges (although not necessarily to the global optimum). As with any partitional clustering algorithm, it is highly sensitive to the initial parameters: the number k of clusters and their centroids, respectively. In practice, to obtain the best clustering results, the k-means algorithm is often applied many times to different random initializations and the minimizer of the intracluster variance is selected as the final clustering result. However, this is at the cost of increased computation.

The main focus of this paper is automatic initialization of k-means by estimating both the number of clusters present in the dataset and the positions of initial cluster centroids. Initialization is crucial since improper initialization could lead the algorithm to converge to some undesirable local minimum of the objective function. Initialization of the k-means algorithm has been extensively studied during the last decade (see [START_REF] Peña | An empirical comparison of four initialization methods for the k-means algorithm[END_REF] or [START_REF] Redmond | A method for initialising the k-means clustering algorithm using kd-trees[END_REF] for a short comparative study of initialization methods for k-means). This paper proposes an improved method to initialize the k-means algorithm.

The issue of estimating the right number of clusters has been widely addressed. One of the most widespread approaches is to use statistical model selection criteria such as Akaike information criterion (AIC), Bayesian information criterion (BIC), minimum description length (MDL) [START_REF] Rissanen | Stochastic Complexity in Statistical Inquiry[END_REF], Tibshirani's gap [START_REF] Tibshirani | Estimating the number of clusters in a dataset via the gap statistic[END_REF], Hartigan's index [START_REF] Hartigan | Clustering Algorithms[END_REF] or Krzanowski & Lai's index [START_REF] Krzanowski | A criterion for determining the number of groups in a data set using sum of squares clustering[END_REF]. [START_REF] Milligan | An examination of procedures for determining the number of clusters in a data set[END_REF] have studied and compared thirty different model selection criteria and they concluded that Calinski & Harabasz's index [START_REF] Calinski | A dendrite method for cluster analysis[END_REF] was the best. Most of the aforementioned techniques are based on a clustering step with k-means, followed by a merge or split procedure to optimize these criteria. We briefly review some of these methods below.

Milligan and Cooper

In the vector quantization community, Bischof et al. [START_REF] Bischof | MDL principle for robust vector quantisation[END_REF] proposed an algorithmic complexity called MDL, to select the number of clusters (starting with a large value of k). Pelleg and Moore [START_REF] Pelleg | X-means: Extending K-means with efficient estimation of the number of clusters[END_REF] developed a variant of k-means referred to as X-means, starting from a minimum number of clusters followed by application of k-means. Their agglomerative cluster splitting procedure is based on the BIC criterion. Hamerly and Elkan [START_REF] Hamerly | Learning the k in k-means[END_REF] proposed to learn the integer k required in k-means by assuming a Gaussian mixture model. Starting with a small value of k, the algorithm splits the clusters that fail a test of spherical Gaussianity. Between each statistical test, the k-means algorithm is applied to refine the solution. Nevertheless, this algorithm referred to as G-means, performs poorly in the presence of non-spherical or nonelliptical clusters. Tibshirani et al. [START_REF] Tibshirani | Estimating the number of clusters in a dataset via the gap statistic[END_REF] defined the Gap statistic to determine the optimum number of clusters. Their method is based on the output of one of two clustering algorithms: k-means or hierarchical classification. The authors propose to compare the logarithm of the within-cluster dispersion to its expectation under a reference null-distribution. This method performs better if clusters are well separated and if the number of clusters is small. The location of the initial cluster centroids has serious impact on the performance of partitional clustering algorithms. Forgy's approach [START_REF] Forgy | Cluster analysis of multivariate data: Efficiency vs. interpretability of classifications[END_REF] consists in choosing initial centers by picking them at random. This method often converges to local minima. A refinement was recently proposed by Arthur and Vassilvitskii [START_REF] Arthur | k-means++: the advantages of careful seeding[END_REF] called k-means ++. The first initial center is randomly chosen among the data points. The next seed is selected among the remaining data points with a certain probability, depending on the distance between data points and their closest center. This operation is repeated until k centers are chosen. [START_REF] Katsavounidis | A new initialization technique for generalized lloyd iteration[END_REF] starts by initializing the first cluster as the point with the maximal norm. Then, the distance of all data points from the first cluster center is calculated. The point with the largest distance is chosen as the second cluster center.

Katsavounidis et al.'s algorithm

After that, for every data point, the distance with its closest center is computed; and the point with the largest distance is chosen as the next cluster center. This process is repeated until k centers are chosen. Unfortunately, this method is sensitive to outliers. Kaufman and Rousseeuw's method [START_REF] Kaufman | Finding groups in data : an introduction to cluster analysis[END_REF] selects the first center as the most centrally located point. Then, a heuristic rule described in [START_REF] Kaufman | Finding groups in data : an introduction to cluster analysis[END_REF] is used to select other clusters centers until k centers are chosen.

A good way to properly initialize the k-means algorithm is to place every cluster centroid at the modes of the joint probability density of the data. Motivated by this idea, Bradley and Fayyad [START_REF] Bradley | Refining initial points for k-means clustering[END_REF] propose a refinement of initial conditions of k-means near the modes of the estimated distribution, by the use of a recursive procedure executing k-means on small random sub-samples of the data. Global k-means was introduced by Likas et al. in [START_REF] Likas | The global k-means clustering algorithm[END_REF], and is also motivated by the same idea that modes play an important role. Their method is based on a recursive partitioning of the data space into disjoint subspaces by using kd trees. They then define a cutting hyperplane as a linear space perpendicular to the highest variance axis. This is similar to the data partitioning that occurs in adaptive vector quantization and regression trees [START_REF] Badel | Tree structured non linear signal modeling and prediction[END_REF].

Our method is motivated similarly to the Bradley and Fayyad's motivation of their recursive algorithm [START_REF] Bradley | Refining initial points for k-means clustering[END_REF]: the modes of the underlying unknown multivariate density are the most pertinent features for initializing a partitional algorithm. The originality of our approach resides in (1) the use of Prim trajectory to map the multivariate density to a singly indexed 1D space, and (2) the formulation of mode detection as a binary hypothesis testing problem (section 2). This allows us to specify a fully unsupervised procedure, which exhibits both reliability and robustness in the presence of noise and outliers. An implementation of our approach to large data sets, relying on the use of a pre-conditioning data driven classification tree, is described in section 3. Finally, in the last section, experimental results are presented, and alternative information theoretic similarity measures are introduced. Information theoretic similarity measures turn out to be better adapted to feature vectors that behave like spectra, i.e., non-negative valued and normalized. Our new Prim-based clustering approach is tested in this context.

Automatic detection of the relevant number of clusters

In this section, the definition of MST and Prim's algorithm are briefly reviewed. Then an original approach is proposed for initializing k-means by jointly estimating the number of relevant clusters and the locations of their corresponding centroids.

Formulation

Let V be a set of N data points or feature vectors in R L . In our unsupervised framework, it is assumed that the points are sampled from an unknown density P (v) and that distinct groups manifest themselves as distinct modes of the density P (v).

Each data point (actually a vector in R L ) is interpreted as a vertex in a graph so that V is the graph vertex set. The goal is to partition V into k classes. Denote C : (C 1 , . . . , C k ) the set of clusters, and let M : (µ 1 , . . . , µ k ) be the set of corresponding centroids. Let us emphasize that no prior information is introduced: in the context of unsupervised clustering, neither the underlying distributions nor the number and locations of meaningful clusters are known.

The k-means algorithm aims at finding a partition of the data set such as the following error function (within-cluster summary distance to centroids) is minimized:

J(C, M ) = k j=1 v∈Cj ||v -µ j || 2 (1)
This minimization is well known in the pattern classification and vector quantization literature [START_REF] Gersho | Vector Quantization and Signal Compression[END_REF][START_REF] Duda | Patter Classification[END_REF]. In this paper the classical Lloyd algorithm is used which is summarized in the final loop of Algorithm 1 given in the sequel.

Minimum Spanning Tree and Prim's algorithm

Let T be an undirected acyclic graph (or tree) connecting all vertices in V . The graph is specified by its set of vertices and its set of edges E = (e(v i , v j ), (i, j) ∈ {1, . . . , N }). The length or weight |e(v i , v j )| of an edge measures a distance or dissimilarity between two vertices

(v i , v j ).
The total length of the tree T is the sum of all edge lengths: L(T ) = e∈T |e|. The minimal spanning tree T * is the one having minimal length among all spanning trees:

L(T * ) = min T e∈T |e| (2)
There exist many algorithms for constructing the MST (see [START_REF] Graham | On the history of the minimum spanning tree problem[END_REF] for state of art of MST's algorithms). This paper focuses on Prim's algorithm [START_REF] Prim | Shortest connection networks and some generalizations[END_REF] as some of its features will be used for estimating the presence of relevant clusters in the data set. Prim's algorithm is briefly described below.

Let T i be the partially connected graph at iteration i -1, hence i vertices are connected to the MST. At the ith iteration, one non-connected vertex, say v i , and one connected vertex of T i are selected, so that the dissimilarity measure between them is minimal. T i becomes T i+1 with this new vertex v i and the associated edge of minimal length. This operation is repeated until no unconnected vertex remains; N -1 such iterations are required to construct the complete MST from a set of N points. A consequence is that any vertex is connected to its nearest neighbor, and that any tree fragment is connected to its nearest neighbor by the shortest possible link.

A MST is unique if there are no ties in the pairwise distances. The choice of the initial vertex has no consequence on the resulting MST: the same set of edges will be obtained, although they will be constructed at different iteration indexes. The construction of the MST associated with a vertex set V only requires the knowledge of the dissimilarities between the pairs of vertices. These will be assumed to be symmetric quantities in the sequel, in agreement with the fact that only undirected graphs will be considered. The importance of the choice of the metric provided to the MST computing routine is discussed in Section 4.2.

Exploiting the MST and the Prim's Trajectory

The use of the MST for describing multivariate data in R L in order to exhibit clusters is not a new idea [START_REF] Theodoridis | Pattern Recognition[END_REF]. Existing methods share a common feature: they are based on graph cuts.

Single-cut clustering [START_REF] Gower | Minimum spanning trees and single linkage cluster analysis[END_REF] is easily performed using the MST, but it is known be highly unstable in the presence of noise and/or outliers in the data set. Slagle et al. [START_REF] Slagle | Experiments with some cluster analysis algorithms[END_REF] proposed to build a one dimensional function of the length of edges connecting vertices within a short subgraph, to detect the largest ones, subsequently applying graph cuts to form clusters. Stuetzle [START_REF] Stuetzle | Estimating the cluster tree of a density by analyzing the minimal spanning tree of a sample[END_REF] proposed to consider multiple cuts involving all edges larger than a threshold in order to preserve minimal cardinality (the 'runt size') of the resulting clusters (subgraphs). Tuning the adequate 'runt size' turns out to be difficult and critical for accurate clustering performance.

The Prim trajectory makes use of the MST construction scheme to propose candidate subsets of vertices that are representative of modes of the density P (v). Instead of trying to detect edges to cut the Prim algorithm decides which vertices to keep, which we use to form the seeds of clusters. Like the method presented in [START_REF] Slagle | Experiments with some cluster analysis algorithms[END_REF] the Prim trajectory projects the data into a one dimensional space.

Prim trajectory and the modes of P (v)

We formally define the Prim trajectory here. Prim's algorithm starts with a root vertex and then constructs edges by a sequence of nearest neighbor operations. Let g(i) = |e i | be the length of a new edge built by Prim's algorithm at iteration i. The ordered set of edge lengths [g(i), i = 1 . . . N -1] is defined as Prim trajectory. The edge function g allows us to "unfold" the MST built in an L dimensional space into a one-dimensional function (see Fig. 1).

A rough estimate of the L-variate density P (v) in the vicinity of v i and v j may be inferred from the weight of the edge connecting them: if v i and v j are connected by the MST then, for k = 1, the 1-NN density estimate is given by p(v i ) ∝ 1 minj |e(vi,vj)| L . This widely used density estimator motivates our approach: a series of consecutive connections with low weight (a valley in g) is the signature of a peak in the density estimate of P (v). Then, if enough vertices are involved in the valley, a mode of P (v) will be detected by thresholding g.

Discussion

• While the obtained MST does not depend on the vertex onto which it is rooted, g changes with the choice of the root. However, Prim's algorithm aggregates new edges of minimal weight at each iteration and the vertices sampled from the same mode will always be connected within a subset of consecutive iterations. Figure 1 illustrates this property. valleys, associated to series of consecutive iterations connecting close neighbors from the same mode.

• The choice of the threshold to apply to g(i) is crucial, since it defines the sensitivity of the mode detection to overlapping of clusters. As g is built from a finite size sample, a high threshold will lead to detection of spurious modes, whereas a low threshold may lead to missing some modes. The example in [START_REF] Michel | Unsupervised clustering with MST: Application to asteroïd data[END_REF] illustrates this behavior; the threshold here was arbitrarily set to ǫ t = σ(e i ), the standard deviation of the edge lengths. In Fig. 1(e) and (f), if a mode is associated to each group of successively connected points below the threshold, a lot of small spurious modes are introduced due to the high variability of the 1 NN density estimator for small N . This was addressed and an entropy based approach was discussed in [START_REF] Michel | Unsupervised clustering with MST: Application to asteroïd data[END_REF] for the case of Euclidean distances, but the resultant algorithm has high computation cost. We formulate a systematic approach to threshold selection using a binary hypothesis testing framework in the next subsection.

Setting the threshold

Consider the following question. What is the probability that a subset of vertices appear as a close connected neighbors in the absence of mode, thus leading to the erroneous detection of a mode? The closeness will depend on the given threshold, and the erroneous detection will be a false alarm. If we specifiy a null model for P (v) a Neyman-Pearson detection framework can be adopted. We specify the hypothesis testing procedure as follows. First, we propose a model for the spatial vertex distribution under the null hypothesis H 0 (no cluster). Second, we express the probability of finding a subset of k consecutive connections whose weights are below a given threshold ǫ t . Third, we relate the probability of false alarm to the parameters N , ǫ t and k.

Modeling the vertex distribution under H 0 .

Under H 0 we model the N vertices as a realization of a homogenous Poisson process over the bounded domain V ⊂ R L . This implies that under H 0 and given N , the vertices are samples of a multi-dimensional independent and identically distributed (iid) uniformly distributed over V. Note that modeling the distribution of vertices of a MST as a Poisson process was assumed by M.Steele in [START_REF] Steele | Probability theory and combinatorial optimization[END_REF].

Let v i be the vertex connected at iteration i in the Prim construction algorithm. Let p N,ǫ,i be the probability under H 0 that another vertex is found within B (v i ,ǫ) of radius ǫ, centered on v i :

p N,ǫ,i = B (v i ,ǫ) f (r)dr
where f (r) = 1/vol(B(v i , ǫ)) and C L is the volume of the L dimensional unit ball. Conditioned on the vertex v i , define k the number of vertices that fall in the ball B (v i ,ǫ) , excluding vertex v i . As the vertices are assumed to be iid samples conditioned on N , k is binomial distributed with parameters N and p N,ǫ,i . Therefore, the probability F N,ǫ,v i (k) = Pr(k ≤ k|N ) is the probability that k vertices are in B (v i ,ǫ) and is given by a binomial probability law with parameters N , p N,ǫ,i , whose cumulative distribution function is

F N,ǫ,vi (k) = k j=0   N -1 j   p j N,ǫ,i (1 -p N,ǫ,i ) N -1-j
It is well known that the convergence of the binomial converges to the Poisson probability with rate parameter β as N → ∞ and ǫ → 0, with β = lim N →∞ (N p N,ǫ,i ). Consequently, we obtain

F N,ǫ,vi (k) N →∞ -→ e -β k j=0 β j Γ(j + 1)
Note that 1 -F N,ǫ,vi (k) is simply an Erlang distribution with rate β.

Since f (v) = 1/|V| where |V| stands for the volume of V, the Poisson rate is

β ≃ N C L |V| ǫ L = λǫ L Parameter λ = N CL
|V| emphasizes the respective importance of the number of observations (N ) and the volume of the support (V).

Probability of detecting a mode under H 0

With the results above the probability that there exists at least one edge (or weight) of length smaller than ǫ is

F v i (ǫ) = 1 -e -β (3)
Then the probability that under H 0 , k consecutive connections of length smaller than ǫ and that the (k + 1)-th has a larger length is given by

P k,ǫ = (1 -F v i (ǫ)) k i=1 F v i (ǫ) = (1 -F v i (ǫ))(1 -e -λǫ L ) k (4) = 1 -e -λǫ L k e -λǫ L (5)
Note that the obtained expressions above satisfy the expected asymptotic equalities:

lim β→∞ P k,ǫ = 0 for any k < N value, lim β→0 P kǫ = 0 for k = 0.
Furthermore, it is easily verified that

∞ k=0 P k,ǫ = e -β ∞ k=0 (1 -e -β ) k = 1.
The expression relating β to ǫ needs to be reconsidered here, in the context of a MST construction with Prim's algorithm. Actually, one connects a new vertex to an existing growing subtree. In the limit of large N and large subtrees, this amounts to testing the presence of a vertex in a neighborhood which is no longer a sphere but a half sphere. This is illustrated on Fig. 2. Consequently the normalization constant C L (volume of the unit sphere in the L-dimensional Euclidean space) will be replaced by CL 2 . The above approximations were derived under the assumptions of N and large Prim subtree size. We validate that this asymptotic approximation is accurate for finite sample size by comparing the simulated average number of vertices that are consecutively connected by edges of length smaller than a threshold ǫ. In the previous derivations, k denotes the number of vertices in a neighborhood centered on a vertex v, but excluding v; the number of edges connecting successive neighbors and the number of vertices involved in a cluster therefore differ by one vertex. As a consequence, the mean number of these vertices, denoted by < k > is obtained as the average of (k + 1) where k is the number of edges whose length are less than ε. As a cluster will be considered only if there is at least one such edge, we obtain

< k > = ∞ k=1 (k + 1)P k,ǫ = e -λǫ L ∞ k=1 k(1 -e -λǫ L ) k + 1 -e -λǫ L = (1 -e -λǫ L ) + (1 -e -λǫ L ) e -λǫ L = 2 sinh(λǫ L ) = 2 sinh C L 2 N |V| ǫ L . (6)
Fig. 4 shows the average size of candidate modes detected as a function of the threshold value. All the detections are false alarms, as the vertices are drawn from a uniform distribution over [0, 1] 2 . This simulation suggests that there is good match between the asymptotic theory and the finite N simulations.

Threshold as a function of the probability of false alarm (PFA)

A false alarm mode occurs whenever some edges of length smaller than ǫ are created under the null hypothesis that the vertices are realizations from a homogeneous Poisson process.

Note that the null hypothesis then refers to the no-cluster hypothesis. To avoid these false alarms, a simple strategy is to apply additional criteria. In the sequel we will only declare a mode present if there are least k 'small' edges that are contiguous in the Prim trajectory. This requires relating k to the probability of false alarm, denoted PFA.

From (3) the expression of the false alarm probability, the probability of occurence of at least k consecutive connections of length less than ǫ, is given by:

P F A (k, ǫ) = 1 -e -λǫ L k . ( 7
)
If L-dimensional Euclidean space is considered, the volume of the half sphere of radius ǫ is

B L (ǫ) = 1 2 C L ǫ L .
Finally, since λ = N CL 2|V| :

P F A (k, ǫ) = 1 -e -C L 2 ǫ L N |V| k (8)
This formula determines the relationship between k m (minimum number of vertices required to detect a mode) and the threshold value ǫ in the framework of a Neyman-Pearson test (Fig. 3), for a fixed PFA.

k = log(P F A ) log 1 -exp -CL 2 ǫ L N |V| (9)
The ratio N |V| in equation ( 9) reveals the role of the density of vertices. By maintaining k and the PFA constant over the sample space, equation (9) specifies value of the threshold that assures a given false alarm rate PFA. In practice, for a given PFA one has to either fix the threshold ǫ and deduce the minimal k, or fix the minimal k and compute the threshold.

The proposed clustering algorithm is summarized in Algorithm 1. As illustrated in section 4, this algorithm performs well for identifying distinct groups in datasets. Assign v to nearest cluster C j where Index(v i ) = argmin j∈1,...,k ||v iµ j ||.

8:

Recompute mean µ j as µ j = v i ∈C j v i card(Cj) . 9: until convergence 10: return Index 13

Implementation of MST's for large datasets

The construction of the Prim-s trajectory or the MST in the previous section requires evaluation of a ranking of the set of all pairwise distances. If the sample set of N vertices is in an L-dimensional space, this will induce a computational load of O(L N 2 /2) flops for evaluating the similarity matrix, and an additional O(N log(N )) logical operations for the ranking step (for example by using the "quick-sort" algorithm [START_REF] Press | Numerical Recipes 3rd Edition: The Art of Scientific Computing[END_REF]). This clearly becomes prohibitive for large data sets.

However, the Prim's trajectory construction can be reformulated as a series of nearest neighbor searches. Finding nearest neighbors of a given point in space within an efficient computation time has been widely addressed. One of the most famous search algorithm is the space partitioning technique based on a k dimensional binary search tree, referred to as the kd tree [START_REF] Friedman | An algorithm for finding best matches in logarithmic expected time[END_REF]. The optimized kd tree method chooses the hyperplane passing through the median point perpendicular to the coordinate axis where the distribution exhibits the largest spread. Based on a search for nearest neighbors using the kd tree and priority queues to reduce these searches, Bentley and Friedman [START_REF] Bentley | Fast algorithms for constructing minimal spanning trees in coordinate spaces[END_REF] proposed several algorithms of MST construction. Guttman [START_REF] Guttman | R-trees: A dynamic index structure for spatial searching[END_REF] proposed the R-tree : Although kd tree methods use partitioning element that is a hyperplane, the R-tree uses hyper-rectangular regions. However these methods cannot easily handle high-dimensional data, since their complexity exponentially grows with the dimension. Other methods such as vp trees [START_REF] Yianilos | Data structures and algorithms for nearest neighbor search in general metric spaces[END_REF], quad-trees [START_REF] Finkel | Quad trees a data structure for retrieval on composite keys[END_REF] or hB-trees [START_REF] Lomet | The hB-tree: a multiattribute indexing method with good guaranteed performance[END_REF] have the same computational difficulty for large data sets.

Here, we reduce the computational load by using a preconditioning data-driven hierarchical classification tree, a procedure we denote as CT. The tree can be learned from a small randomly chosen subset of R data points involving only O(LR log R) logical operations. It is designed to easily identify the neighborhoods of each vertex by comparing the vertex coordinates along each of the L dimensions to a set of L thresholds. The size of the neighborhood can be set up in such a way that in the average, each neighborhood contains approximately M ≪ N vertices. Consequently, the number of pairwise distances that need to be evaluated is of the order of M 2 /2. Thus the computational load is lowered to the order of O(N L M 2 /2) flops and O(N M log M ) logical operations, giving substantial computational savings when M ≪ N .

The constructions of MST using the pre-conditioning classification tree (CT) described below will be referred to as "Nearest-Neighbor MSTs" (NN MST) in the sequel.

Assume that after l iterations of the CT growing procedure a partition Π was created. At the next iteration Π is refined by splitting each of its constitutive cells, called parent cells, into 2 L children cells. This step is conducted by employing a median split over each axis. In order to control the number of cells, only the cells with a minimum number (n min ) of vertices are split further, the others being declared terminal cells, or leaves in the CT. The advantage of such a procedure is that the cells tend to fall in areas of the sample space where the sample density is high. Furthermore, the CT procedure allows to maintain an approximately constant number of vertices in every leaf of the tree. The median splitting rule allows to define subcells with approximately equal marginal probabilities for arbitrary underlying distribution.

A criterion for setting n min using the variance of the median estimate was derived in [START_REF] Mood | Introduction to the Theory of Statistics[END_REF] and refined in [START_REF] Badel | Tree structured non linear signal modeling and prediction[END_REF]. The same idea is used here, and developed for each of the L axes. The L dimensional CT is obtained by merging all 1-dimensional CTs by a logical 'and'.

Consider a cell Π in the L-dimensional space. Then Π = × L j=1 π j where π j = [α j , β j ]. π j is a 1-D cell defined on axis j; let N πj be the number of vertices in π j . The coordinates are assumed to be iid with marginal continuous density function

f v j |π j on the interval [α j , β j ].
The sample median Tπj is asymptotically Gaussian distributed [START_REF] Mood | Introduction to the Theory of Statistics[END_REF]:

Tπj ∼ N T πj , 1 4N πj (f v|π j (T πj )) 2 (10)
where T πj is the theoretical median: T πj = (β j + α j )/2. A criterion to select n min is specified by requiring that N πj is large enough to ensure that the density of Tπj has maximum mass inside the interval [α, β]. It will be assumed that the sample medians Tπj , j ∈ {1, . . . , L} are statistically independent. Furthermore, it is assumed that f v|π j is a uniform distribution over π j (which will be practically true for cells whose size is small compared to the variations of the underlying distribution) we obtain 1 -

L j=1 P r(| Tπj -T πj | ≤ (β j -α j )/2) ≤ ǫ. ( 11
)
Under the Gaussian approximation [START_REF] Shi | Normalized cuts and image segmentation[END_REF], and integrating over all coordinate axes, we obtain

1 -P r(|Z| < N Π ) L ≤ ǫ
where Z is a standard normal random variable (zero mean and unit variance).

P r(|Z| < N

Π ) = erf ( N Π /2)
where erf is the error function:

erf (x) = 2 √ π x 0 e -t 2 dt.
This motivates the following stopping rule for the CT cell subdivision:

if N π ≥ 2 erf -1 ((1 -ǫ) 1/L ) 2 = n min then subdivide π else stop
The structure for the CT associated to the L-dimensional distribution is obtained by storing every 1-dimensional CT. Each CT is associated to an array describing the hierarchy and the thresholds. Thus, nearest neighbors of a vertex are found by performing a small number of logical tests to identify the terminal cell to which it belongs, thus sparing the costly flops needed to evaluate all N 2 dissimilarities.

Results and Discussions

In this section,we illustrate the proposed clustering algorithm for three applications. The first application is real-world data extracted from the UCI Machine Learning Repository [START_REF] Asuncion | UCI machine learning repository[END_REF].

The second application is the segmentation of multi-spectral satellite image. The third application is to classification of chemical species in hyper-spectral imaging of planet Mars.

Experiments on the estimation of the number of clusters

Estimation of the number k of clusters is performed on a set of simulated data, for which the actual k is known. All previous methods rely upon the analysis of the results of a clustering algorithm to select k. We used the k-means algorithm and values ranging from k = 1 to k = 10 were tested. The tests were conducted on 4 different data models described below, and 50 independent data samples were generated for each model.

Model 1: three spherical clusters in R 2 . The clusters were generated by simulating a mixture of three spherical normal variables with means (0, 0), (0, 5), (5, -3).

Model 2: four spherical clusters in R 3 . Each cluster is randomly chosen to contain 25 or 50 observations, with means randomly chosen from a N ([0, 0, 0], 5I 3 ) distribution, where I 3 is the identity matrix in R 3 .

Model 3: four spherical clusters in R 10 . Each cluster is randomly chosen to contain 25 or 50 observations, with means randomly chosen from a N (0, 3.6I 10 ) distribution.

Model 4: two elongated clusters in R 3 . Cluster 1 = {x + b/x 1 = x 2 = x 3 = t, t ∈ [-.5, .5]}
where b is a centered Gaussian noise in R 3 , with covariance matrix .1I 3 ; Cluster 2 is similarly obtained except that the noise is centered in m = [START_REF] Shi | Normalized cuts and image segmentation[END_REF][START_REF] Shi | Normalized cuts and image segmentation[END_REF][START_REF] Shi | Normalized cuts and image segmentation[END_REF]. The two clusters are thus stretched along the main diagonal.

For comparison we also simulated the methods of Calinski and Harabaz, Hartigan, Tibshirani, and Kaufman and Rousseeuw. Calinski and Harabasz [START_REF] Calinski | A dendrite method for cluster analysis[END_REF] proposed a test based on the optimization of a ratio between the inter and intra-cluster variances. Hartigan [START_REF] Hartigan | Clustering Algorithms[END_REF] derived a method that evaluates the relatives variation of the intra cluster sum of squares when k increases. Tibshirani et al. [START_REF] Tibshirani | Estimating the number of clusters in a dataset via the gap statistic[END_REF] constructed the Tibshirani's Gap: a criterion that relies on the intra cluster measure, but introduced a corrective term accounting for the statistical variation of the intra cluster measure under a null hypothesis. Kaufman and Rousseeuw [START_REF] Kaufman | Finding groups in data : an introduction to cluster analysis[END_REF] pursued analogous ideas, measuring the average intra cluster distance and the minimal inter cluster separation as contrasted to variances.

A description of the simulation is given on Fig. 5. For each of the 50 data samples and for every model, the estimation of the number of clusters present in the data set according to the various criteria described above is recorded. The obtained results are reported in Tab.1.

Our method was applied with fixed parameters ε = σ e , and k m = 3, where σ e stands for the standard variation of the edge-length distribution over the MST. (Note the corresponding PFA is obtained by applying equation ( 8)). The method proposed in this paper exhibits improved performance as compared to other methods for all the dimensions of the sample data studied here.

Similarity measures

The choice of a particular distance measure between points has a great impact on the resulting cluster identification. Using different metrics leads to different MSTs, and it makes sense to adapt the dissimilarity measure to the physical structure of the sample data. It is for instance well known (see e.g. [START_REF] Chang | An information-theoretic approach to spectral variability, similarity, and discrimination for hyperspectral image analysis[END_REF]) that for feature vectors that are spectral(-like), information divergences dissimilarity measures can lead to better results. Furthermore, these divergences handle the case of missing data, e.g., missing spectral measurements at a given wavelength. There are many information divergence measures, e.g. Csiszàr [START_REF] Csiszár | Information-type measures of difference of probability distributions and indirect observation[END_REF] or Bregman [START_REF] Bregman | The relaxation method of finding the common point of convex sets and its application to the solution of problems in convex programming[END_REF]).

Here we focus on two particular measures: the symmetrized Kullback-Leibler and the Rényi divergences [START_REF] Cover | Elements of Information Theory[END_REF]. The symmetrized Kullback-Leibler divergence is defined as:

d KL (x, y) = L i=1 ( xi -ỹi ) log xi ỹi . (12)
Alternatively, the symmetrized α Rényi divergence (α satisfies 0 < α < 1) can similarly be used as a dissimilarity measure:

d α (x, y) = 1 α -1 log L i=1 xi α ỹi 1-α + log L i=1 ỹi α xi 1-α . (13)
Properties and advantages of the Rényi divergence have been detailed by Hero et al. [START_REF] Hero | Applications of entropic spanning graphs[END_REF].

Note that when α → 1, the α-divergence (13) converges to the Kullback-Leibler divergence [START_REF] Ball | ISODATA, a novel method of data analysis and classification[END_REF]. Divergences do not satisfy the triangular inequality in general ; they are semi-metrics.

In the remote sensing community, alternative similarity measures have been proposed for hyper-spectral data. One such measure is based on the Spectral Angle Mapper (SAM) [START_REF] Keshava | Distance metrics and band selections in hyperspectral processing with applications to material identification and spectral libraries[END_REF]:

θ(x, y) = arccos x, y x y , (14)
where ., . is the dot product, and . is the Euclidean norm. The angle θ belongs to the interval (0, π/2) This measure is used to define the angle existing between spectra. It measures a similarity rather than a dissimilarity.

These measures may be used for evaluating dissimilarities, in the case where 'probabilitylike' or 'spectral-like' data set is analyzed. Such a choice is discussed in Sections 4.1 and 4.3.

In this paper, the focus is on the k-means algorithm. k-means can be adapted to match the metric used in the construction of the tree, similar to the Bregman clustering approach [START_REF] Dhillon | A divisive information-theoretic feature clustering algorithm for text classification[END_REF][START_REF] Banerjee | Clustering with Bregman divergences[END_REF].

The assignment step in Algorithm 1 then uses the Bregman divergence instead of L2 norm 1 .

Clustering experiments with Prim initialization of k-means

First, we present results of comparisons for real data from the UCI Machine Learning Repository data sets [START_REF] Asuncion | UCI machine learning repository[END_REF] . We compare our method with other k-means initialization methods:

Forgy's approach (randomly) [START_REF] Forgy | Cluster analysis of multivariate data: Efficiency vs. interpretability of classifications[END_REF], k-means ++ [START_REF] Arthur | k-means++: the advantages of careful seeding[END_REF], Global k-means [START_REF] Likas | The global k-means clustering algorithm[END_REF]; Katsavounidis et al. (KKZ) [START_REF] Katsavounidis | A new initialization technique for generalized lloyd iteration[END_REF]; Kaufman and Rousseuw's method [START_REF] Kaufman | Finding groups in data : an introduction to cluster analysis[END_REF] and Bradley and Fayyad's algorithm [START_REF] Bradley | Refining initial points for k-means clustering[END_REF]. Since for these UCI data sets the actual number of clusters is known, it is used for all 1 Information theoretic divergences and Euclidean distance belong to the family of Bregman divergences.

methods requiring it as an input parameter. Note that for our methods, this value is estimated (with k m = 3, ε = σ e ). The comparisons presented on Tables 2,3 The reader will notice that our Prim initialization methods performs very well on the Iris data, both in terms of distortion and accuracy measure and competes with other k-means initialization methods. Information divergences appear to be more interesting for the other data.

On the Wine data set, results obtained with the Prim initialization both with Euclidean and informational divergences compete very well with other standard initialization methods. On the image segmentation data set, which contains more instances and more clusters than the previous data sets, our Prim initialization method produces significantly better cluster performance than obtained with classical initialization methods. For all these tests, the threshold of the Prim trajectory has been set to the standard deviation of the edge length, as explained above. The minimal number of vertices required to accept a detection of a mode was set to maintain the PFA to a value below 10%, where this minimal number was determined using our asymptotic Poisson approximation presented in Section 2.

Multi-spectral image of Paris

We next consider a Landsat Thematic Mapper multispectral (7 bands) image of size 512 × 512 with 30 meter resolution (Fig. 6)2 . Each image was recorded from a device operating at a different wavelength. The images cover an area corresponding to the city of Paris, France. The 7 images are perfectly registered. The cluster analysis is performed on the four largest principal components of the 7 bands, so in the following only four dimensional features are considered instead of the initial seven bands. The affinity measures used were information divergences between the 4-point spectra associated with each pixel, considered as feature vectors. The Nearest Neighbor MST algorithm described in Section 3 was applied sub-sampling the images by a factor of 4. The pixels were assumed to be independent, and no information related to the image structure (e.g. adjacency of pixels) was used in this test. Such prior information could be used as additional feature information to improve clustering but this was not investigated here.

Fig. 7 shows the 8 identified clusters, for which 3 are easy to understand. Cluster 2 is characteristic of trees and grassy regions (one can recognize recreation areas and natural parks in Paris surroundings (Boulogne and Vincennes) and trees along the shores of the Seine).

Cluster 3 exhibits the "water areas" in Paris, and the Seine river together with some known ponds are easily extracted. Cluster 4 is clearly associated with roads, asphalt and concrete.

Other clusters cannot be fairly interpreted without cross validation with, for example, pollution imaging or gas detection systems.

As the "groundtruth" cluster representation of the data is unknown, the quality of the different classifications was assessed by evaluating the Davies-Bouldin (DB) index [START_REF] Davies | A cluster separation measure[END_REF] to complement the average intra cluster variance, which is minimized by K-Means. DB index accounts for the separation of the clusters:

DB(C) = 1 k k l=1 R l , where R l = max j =l R jl (15)
where R ij = si+sj δij , and with δ ij = d(µ i , µ j ), and

s i = 1 |Ci| v∈Ci d(v, µ i ).
Thus, a lower DB index indicatesa better clustering result. Table 5 summarizes the results. Our Prim initialization approach combined with informational divergences gives the best results overall. representation of the results consists of binary masks (the pixel is white if it belongs to this cluster). H 2 O ice compounds are on the peripheral of CO 2 ice compounds, which have some sense since at these places the CO 2 sublimates to reveal H 2 O ice. The corresponding mean spectrum can also be seen in Fig. 8. Our results are fully consistent with the clustering results obtained by Moussaoui et al. [START_REF] Moussaoui | On the decomposition of mars hyperspectral data by ICA and Bayesian positive source separation[END_REF].

Conclusions and Future works

In this paper, we have proposed a new approach for clustering multi-dimensional data. The method is based on the estimation of the number of clusters and the centers of the clusters from the Prim construction of a minimum spanning tree, followed by an initialization of the classical k-means clustering algorithm. New criteria were derived for setting the false alarm rate of a test over the Prim's trajectory associated with a MST built over the set of data. The false alarm rates were derived under the null hypothesis that the data points are distributed as a homogeneous Poisson process. We demonstrated the utility of the information divergence as a dissimilarity measure for astrophysical multi-spectral image analysis. In this paper, the threshold value is constant along the Prim's trajectory. This rate could be easily made a function of the position of the connected vertices. This has not been developed here but opens new perspectives.

Dimension reduction and its relationship to spectral clustering methods applied to graphs using information divergence or MST-based distances (for example, dual rooted tree distances) could also be investigated [START_REF] Belkin | Laplacian eigenmaps for dimensionality reduction and data representation[END_REF][START_REF] Bach | Learning spectral clustering, with application to speech separation[END_REF]. Clusters obtained with our method: (d) H2O ice, (e) CO2 ice and (f) dust.

Algorithm 1

 1 Initilization of k-means with Prim's algorithm Input: V : data set of N points ∈ R L d: dissimilarity matrix (N × N ) Output: Index: matrix of labels 1: Construct the Prim trajectory on the dataset V using d.

2 :

 2 Threshold g with ǫ = σ(e i ).

3 :

 3 Evaluate k for a chosen PFA (e.g. 0.05), (k is the minimum number of vertices that should be connected for a detected mode); (Eventually re-adjust ǫ to obtain k in the desired range).

4 : 5 :

 45 Compute of the number k of modes , after discarding all candidate with less than k connected vertices Compute of the centers of mass µ of every detected mode C.

Let x = {x 1

 1 , . . . , x L } and y = {y 1 , . . . , y L } be two non negative feature vectors in R * L . Let xi = xi L j=1 xj and x = {x 1 , . . . , xL } as the normalized quantities, respectively. We define ỹ similarly.

, 4 ,

 4 are based on the error of k-means defined by equation (1) (referenced as distortion in the Table) and the accuracy measure (percentage of correctly labeled points). For these UCI data labels of data sets are only used to assess the error of the different methods. The results obtained are given in Tables 2, 3, 4. The UCI data sets used are described below. The Iris data contains three classes (versicolor, virginica, setosa) of 50 instances each in four dimensions, where one class is linearly separable from the others and the two other classes interleave.The Wine data: contains 178 instances in thirteen dimensions, with three clusters (different cultivars) of different sizes (59, 71, 48). The features are the quantity of constituent found in every type of wine. The Image Segmentation contains 2310 instances in 19 dimensions which describe aspects from 7 classes of images (brick, sky, foliage, cement, window, path, grass). Each class is represented by 330 instances.

Figure 1 :

 1 Figure 1: (a) Toy data set, (b) construction of a MST, (c) extraction of the relevant groups of vertices present below the threshold in the Prim's trajectory, (d) clustering with k-means , (e,f) Prim's trajectories with different initial vertex used in the MST construction.

Figure 2 :

 2 Figure 2: Left: C k denotes the contour of the support where connected vertices are found. The neighborhood which is considered for finding a vertex that could be connected to vi is shown. In the limit of large N , this latter neighborhood is approximated as the half sphere laying on the tangent (hyper)plane to C k (right).

Figure 3 :Figure 4 :

 34 Figure 3: Minimum number of points in a detected cluster as a function of the threshold value ǫ from sets of uniformly distributed vertices over [0, 1] 2 , for a PF A = 0.05 and N = 128, 256, 512, 1024 vertices respectively.

Figure 5 :

 5 Figure 5: (a) Model 1: N = 150; (b) Model 2: N = 175; (c) Model 3: N = 175; (d) Model 4: N = 202

Figure 6 :

 6 Figure 6: Multi-spectral image of Paris

Figure 7 :

 7 Figure 7: Clusters obtained on the multi-spectral image of Paris: (a) cluster 1, (b) cluster 2: trees, parks , (c) cluster 3: water elements, (d) cluster 4: roads, (e) cluster 5, (f) cluster 6, (g) cluster 7 and (h) cluster 8.

Figure 8 :

 8 Figure 8: Mean and ± standard deviation spectra obtained: (a)H2O ice, (b) CO2 ice, (c) dust.

Table 1 :

 1 Estimation of the number of clusters (k * being the true number of clusters)

Table 2 :

 2 Results obtained on the Iris data set

	Method	Distorsion Accuracy
	Prim Initialization (Euclidean)	1.97	0.8933
	Prim Initialization (Kullback-Leibler)	3.83	0.6667
	Prim Initialization (Rényi)	3.11	0.6667
	Forgy (random)	2.99	0.8933
	Bradley Fayyad	1.41	0.8267
	KKZ	1.42	0.8467
	Kaufman Rousseuw	1.40	0.8933
	k-means ++	2.01	0.8867
	Global k-means	1.41	0.8933

Table 3 :

 3 Results obtained on the Wine data set

	Method	Distorsion Accuracy
	Prim Initialization (Euclidean)	2.05 ×10 4	0.6742
	Prim Initialization (Kullback-Leibler) 1.92 × 10 4	0.6910
	Prim Initialization (Rényi)	1.92 × 10 4	0.6910
	Forgy (random)	4.00 ×10 4	0.6966
	Bradley Fayyad	8.15 ×10 4	0.6966
	KKZ	2.92 ×10 4	0.5843
	Kaufman Rousseuw	8.56 ×10 4	0.6573
	k-means ++	4.15 ×10 4	0.7022
	Global k-means	6.86 ×10 4	0.7022

Table 4 :

 4 Results obtained on the Image Segmentation data set

	Method	Distorsion Accuracy
	Prim Initialization (Euclidean)	6.39 ×10 4	0.5173
	Prim Initialization (Kullback-Leibler) 5.20 ×10 4	0.5446
	Prim Initialization (Rényi)	6.39 ×10 4	0.5173
	Forgy (random)	5.52 ×10 4	0.4251
	Bradley Fayyad	4.91 ×10 4	0.5385
	KKZ	1.44 ×10 5	0.3567
	Kaufman Rousseuw	3.79 ×10 4	0.3355
	k-means ++	3.38 ×10 5	0.4771
	Global k-means	2.38 × 10 4	0.4606

Table 5 :

 5 Results obtained on the multi-spectral image

	Clustering Methods	Davies-Bouldin index
	NN MST (Euclidean) + k-means	111.77
	NN MST (SAM) + k-means	157.27
	NN MST (Kullback-Leibler) + k-means	74.87
	NN MST (Rényi) + k-means	88.02
	Expectation Maximization (Mixture of Gaussians)	160.19
	k-means (randomly initialized)	155.13

The image can be downloaded at http://physics.ship.edu/ mrc/astro/NASA Space Science/observe.arc.nasa.gov/nasa/education-/tools/stepby/arcpage/arc20.html.

Acknowledgments

We thank the OMEGA team at IAS/Orsay for his support with sequencing and data download activities. A. Hero's contribution to this paper was partially supported by the US National Science Foundation grant No. CCF 0830490. The authors thank Dr E. Slezak (OCA, Univ. of Nice, France) for usefull discussions on the application of these methods.

Mars hyper-spectral image

Finally we present the results of applying our method to a hyper-spectral image of Mars. This data is very high dimensional and thus serves as a validation that our proposed method is scalable. The other initialization methods could not be compared due to their excessive computational complexity. This Mars image was provided by the Mars Express (European Space Angency) on board imaging spectrometer instrument OMEGA (Observatoire pour la Minéralogie, l'Eau, les Glaces et l'Activité). It covers the south polar cap of Mars in the Martial summer. Two infrared channels are scanned: 128 spectral planes are recorded at wavelength ranging from 0.93 µm to 2.73 µm with a resolution of 0.013 µm and 128 spectral planes from 2.55 µm to 5.11 µm with a resolution of 0.02 µm. Hence, the image is of size 300 × 120 pixels in 256 wavelengths. Some preprocessing of the raw data was applied by "Laboratoire de Planetologie de Grenoble" in order to remove erroneous or unreliable values , e.g., due to faulty pixels of the CCD camera, bad calibration values etc. More information about this dataset can be found in [START_REF] Schmidt | Martian surface classification using omega hyperspectral images. spatiotemporal study of the co2 and h2o seasonal deposits[END_REF].

In recent years, two methods have been introduced to classify chemical species on this hyper-spectral image of planet Mars. The first method developed by Schmidt et al. [START_REF] Schmidt | WAVANGLET: An efficient supervised classifier for hyperspectral images[END_REF] is a supervised approach based on a wavelet representation called wavanglet. This method requires a priori information, specifically the number of classes and the reference spectra (generated from a physical model of chemicals compounds signature). The classification is made by computing the spectral angle in the best subspace of a wavelet filtered space. The results obtained with this method reveals three binary masks (a pixel belongs to one class) of CO 2 ice, H 2 O ice and dust. Note that some pixels have not been classified. As the method requires reference spectra, it requires an a priori report of experts to identify all physical compounds present in the scene, and subsequently generate (or extract from the image) reference spectra.

The second method was developed by Moussaoui et al. [START_REF] Moussaoui | On the decomposition of mars hyperspectral data by ICA and Bayesian positive source separation[END_REF] and is based on blind source separation method referred to as Bayesian positive source separation (BPSS). BPSS requires costly importance sampling, e.g., Monte Carlo Markov Chain (MCMC), for its implementation.

In Fig. 8 clusters obtained with our Prim initialization of k-means are displayed. The dissimilarity measure used was the symetrized KL divergence. The image contains only a few homogeneous groups, so the computation of centroids does not require the whole image.

We resized the image by column and line sub-sampling (one line over 5 and one column over


Dawood Khan
email: khan@laas.fr

Bilel Nefzi
email: bilel.nefzi@loria.fr

Luca Santinelli
email: luca.santinelli@inria.fr

Yeqiong Song
email: song@loria.fr

Probabilistic Bandwidth Assignment in Wireless Sensor Networks

With this paper we offer an insight in designing and analyzing wireless sensor networks in a versatile manner. Our framework applies probabilistic and component-based design principles for the wireless sensor network modeling and consequently analysis; while maintaining flexibility and accuracy.

In particular, we address the problem of allocating and reconfiguring the available bandwidth. The framework has been successfully implemented in IEEE 802.15.4 using an Admission Control Manager (ACM); which is a module of the MAC layer that guarantees that the nodes respect their probabilistic bandwidth assignment as well as the bandwidth assignment policy applied. The proposed framework also aims to accurately analyze the behaviors of communication protocols for energy-consumption and reliability purposes. We evaluate the probabilistic bandwidth assignment methods using CSMA/CA access protocol of IEEE 802.15.4. Furthermore, we analyze the behavior of the ACM and compare the performance of the network with and without using the ACM against the original standard. The simulation results show that the use of ACM increases the overall performance of the network. 4

Introduction

Wireless Sensor Network (WSN) is considered as one of the key technologies for building the future Cyber-Physical Systems (CPS) as it allows today's information systems to monitor and control the physical environment. Often, monitoring and acting through the WSN forms feedback loops within which the control decision should be made in real-time. So the QoS (Quality of Service) of the existing WSNs must be enhanced since most of them mainly focus on the energy efficiency without performance guarantees in terms of bandwidth allocation and end-to-end delay. Designing and analyzing the QoS of WSNs is challenging because of the highly dynamic behavior of WSNs. Moreover, the WSNs need to operate with energy saving policies (e.g., duty-cycled nodes), which further complicates the design and analysis. In general, WSN analysis aims to evaluate the performance limits of a WSN deployment and to guarantee certain QoS. Therefore, it is essential to determine the performance bounds for end-to-end latency, energy, node buffer size, and reliability with respect to the network density, communication protocols, and network topology.

For preserving scalability, most of existing low-power WSNs adopt contention-based MAC. Providing guaranteed resource is even more challenging due to the random nature of channel access method and radio channel behavior; which can be described by a random probability distribution function. In addition, the node reliability/failure is probabilistic due to sudden sensor loss or decay of battery [START_REF] Sastry | Reliable networks with unreliable sensors[END_REF]. Therefore, it is often impossible to provide the deterministic QoS in WSNs.

Furthermore, in WSNs the event/packet (packet generating process) inter-arrivals follow also probabilistic distribution, for example due to the random back-offs in contention based protocols. Jung et al. [START_REF] Jung | Enhanced markov chain model and throughput analysis of the slotted csma/ca for ieee 802.15.4 under unsaturated traffic conditions[END_REF] have shown that the protocols like Carrier Sense Multiple Access with Collision Avoidance (CSMA/CA), modeled using Markov chain, can only provide probabilistic resource guarantees.

Complexity and uncertainty makes WSNs intrinsically probabilistic. Thus, we require analysis methodologies for WSNs which provide the performance guarantees (QoS) with the probabilistic bounds.

In contrast to some recent probabilistic end-to-end delay analysis works [START_REF] Wang | Cross-layer analysis of the end-to-end delay distribution in wireless sensor networks[END_REF], [START_REF] He | Reliability calculus: A theoretical framework to analyze communication reliability[END_REF] which focus on the performance evaluation methods, the probabilistic approach developed in this paper offers a degree of flexibility, while allowing to cope with the dynamics of the environment (channel) and the applications in WSNs. The approach developed can be divided into two parts: i) the analysis, which is a theoretical framework for the analysis of WSNs based on probabilistic network calculus; ii) the Admission Control Manager (ACM), it is a MAC level implementation which ensures that the performance guarantees provided by the analysis part are met. This is achieved in the ACM by actively droping the packets based on the availability of bandwidth (considered probabilistic in this paper) and accepted performance guarantee (QoS) by a node (at the time of composition of the network); in a probabilistic manner.

Related Work

The IEEE 802.15.4 protocol is the popular standard for WSNs and specifies the Medium Access Control (MAC) sub-layer and the physical layer of Low-Rate Wireless Personal Area Networks (LR-WPANs) [START_REF]LAN-MAN Standards Committee of the IEEE Computer Society: Wireless Medium Access Control (MAC) and Physical Layer (PHY) Specifications for Low-Rate Wireless Personal Area Networks (LR-WPANs)[END_REF]. Some works in WSNs [START_REF] Schmitt | Sensor network calculus -a framework for worst case analysis[END_REF][START_REF] She | Deterministic worst-case performance analysis for wireless sensor networks[END_REF], analyze the system with deterministic Network Calculus (NC) [START_REF] Boudec | Network calculus: A Theory of Deterministic Queuing Systems for the Internet[END_REF] to provide worst-case performance bounds for end-to-end latency and per-hop node buffer size. These frameworks rely on deterministic MAC protocols such as Time Division Multiple Access (TDMA) and deterministic routing protocols with worst-case analysis. Most often this contention free communication requires cluster tree topology [START_REF] Jurcík | Dimensioning and worst-case analysis of cluster-tree sensor networks[END_REF] since the mesh topology uses contention based MAC for its reliability and robustness, thus limiting the mesh topologies to provide deterministic worst-case bounds.

WSNs are also required to be adaptive. Most of the approaches over the IEEE 802.15.4, tackle the problem of point-to-point communication within the star with a deterministic model. Relevant time guaranteed communication examples are iGAME [START_REF] Koubâa | An implicit gts allocation mechanism in ieee 802.15.4 for timesensitive wireless sensor networks: theory and practice[END_REF] where the authors propose, using network calculus, a methodology to study the bandwidth allocation problem, and GSA [START_REF] Na | An optimal gts scheduling algorithm for time-sensitive transactions in ieee 802.15.4 networks[END_REF] that tries to minimize the total number of unallocated time-slots by applying a scheduling-based strategy. WSNs have to cope with reactive paradigms, i.e. change the bandwidth allocation according to the events appearing in the environment. Examples of adaptive systems are [START_REF] Huang | An adaptive gts allocation scheme for ieee 802[END_REF] and [START_REF] Nastasi | BACCARAT: a Dynamic Real-Time Bandwidth Allocation Policy for IEEE 802[END_REF], where the authors propose an adaptive solution designed to satisfy realtime constraints and maximize event detection efficiency under variable network load conditions. Focusing on adaptive WSNs, [START_REF] Santinelli | A component-based architecture for adaptive bandwidth allocation in wireless sensor networks[END_REF] advances on the analysis providing a component-based abstraction of WSNs wherein they apply QoS analysis with no strict timing constraints concerns. All these approaches tend to assume or guarantee a worst-case behavior of the WSN application. Instead, recent trends depicts the probabilistic approach, using stochastic network calculus [START_REF] Jiang | A basic stochastic network calculus[END_REF][START_REF] Xie | Stochastic network calculus models under max-plus algebra[END_REF], as flexible enough to cope with the requirements of complex and adaptive WSNs application.

Contributions of the paper:

We propose a probabilistic methodology to dimension a WSN which is modeled with network components. Each network component is identified by a probabilistic arrival stream of data and a probabilistic communication resource guaranteed to it. By hierarchically composing components and probabilistic resource guarantees, it is possible to ease the overall analysis of the system and appropriately configure the network. In particular, we address the problem of allocating and reconfiguring the bandwidth assigned to each node using a MAC level ACM, and we do analysis of MAC-level protocols. The ACM avoid unnecessary transmissions when the bandwidth is not sufficient, which helps to reduce the number of collisions, thus saving energy. Our flexible methodology enables WSN designer to analyze different communication protocols either contention-based like CSMA/CA or contention-free, such as TDMA.

Component based WSN Modeling

Component based design envisions the system as a composition of components. Each component abstracts a process or a physical entity into a black box with interfaces to observe and drive its behavior. We model the WSN nodes as basic computational units which implement functions with temporal or QoS requirements [START_REF] Santinelli | A component-based architecture for adaptive bandwidth allocation in wireless sensor networks[END_REF]. Therefore, WSN nodes can be abstracted as components and the communication (and associated requirements) between nodes representing the link among components. This abstraction makes it easy to analyze and design a WSN; since we can abstract the cluster as one component resulting from the composition of other components. Moreover, component-based approach to the network modeling decomposes the complexity of WSN systems into parts; which can be individually analyzed, thus reducing the overall complexity.

[α u , α l ] [β ′u , β ′l ] [β u , β l ] [α ′u , α ′l ]
The component has an interface associated with it that describes the functional and non-functional behaviors of the component. The interface plays a central role in the component-based design of complex systems such as WSNs; because it defines the notion of composability: two or more components are compatible if they work together properly. Figure 1 depicts a generic WSN component and its component interface abstraction.

System model

We consider a WSN-based CPS for monitoring application where the system consists of a set of sensor nodes collecting data; which is then sent towards collection points (sinks) with bounded transmission delays. We also consider the hierarchical cluster-tree topology for its flexibility and scalability 5 . At the lowest level there are star topologies where a Coordinator (C) manages End-Devices (EDs) to form a leaf cluster. We can obtain a large scale network by extending and interconnecting the star topology clusters in a hierarchical manner creating the so called cluster-tree topology as shown is Figure 2. In the cluster-tree topologies coordinators manage either EDs (nd i ) or other coordinators C k .

Without the loss of generality we assume that coordinators C i do not sense (i.e., their assigned bandwidth is entirely redistributed to the children) and their main function is to maintain the topology and to hierarchically allocate the bandwidth.

Node Component model

Following the reasoning of Network Calculus and Real-Time Calculus [START_REF] Boudec | Network calculus: A Theory of Deterministic Queuing Systems for the Internet[END_REF][START_REF] Thiele | Real-time calculus for scheduling hard real-time systems[END_REF] we define abstraction for WSN components.

Workload and Bandwidth abstraction

The cumulative functions R and S represent the amount of workload and bandwidth respectively; R(t) and S(t) give respectively the amount of workload required and bandwidth available/requested in the time interval [0, t). While R(t) describes a concrete trace of an event/packet stream; the arrival curve α, [START_REF] Boudec | Network calculus: A Theory of Deterministic Queuing Systems for the Internet[END_REF] provides an abstract model which provides an upper-bound on any admissible trace of event/packets in any time interval of length ∆. The same reasoning is applicable to the service S(t) and its service curve β in the interval domain.

Generic aperiodic events are more appropriately modeled with a distribution function (hence a random variable) and its Cumulative Distribution Function (CDF) X, corresponding to the arrival distribution of workload/events (of packet arrival). With probability distributions it is possible to better capture the aperiodicity of aperiodic events. We do not consider any particular CDF in our framework for the packet generating process (simulations most of the time consider a Poisson process) as it is not in the scope of this paper.

Definition 1 (Bounded Workload R + (t)). The "largest "cumulative work function such that the probability of CDF X having event count larger than or equal to R + (t) is lower than a threshold value of Ω.

R + (t) = sup{R(t)|P [X(t) ≥ R(t)] ≤ Ω}, (1)
The extension of the NC to the aperiodic events results in the probabilistic bound where

α(∆) : R + (t) -R + (s) ≤ α(t -s) ∀s < t, (2)
with α(0) = 0 and ∆ = t -s. The probabilistic arrival curve is the couple curve, probability threshold α, Ω , where the probability threshold Ω represents the accuracy of α being the upper bound: the probability that the events are upper bounded is 1 -α.

Definition 2 (Bounded Bandwidth S -(t)). The "smallest "cumulative bandwidth function such that the probability of CDF Y having bandwidth smaller than or equal to S -(t) is lower than a threshold value Λ, where Y is the distribution for bandwidth availability.

S -(t) = inf{S(t)|P [Y (t) < S(t)] ≤ Λ}. (3)
The service curve β lower bounds the available resources in any time interval of length ∆. Therefore, interval based probabilistic resource provisioning curve is represented as:

β(∆) : S -(t) -S -(s) ≥ β(t -s) ∀s < t. (4)
The probabilistic bandwidth provisioning curve is the couple β, Λ where the probability thresholds Λ gives the probability to find bandwidth provisioning below the lower bound β. The model is capable of modeling any possible bandwidth supply in the interval domain, including the bandwidth provisioning by the control coordinators in WSNs. Figure 3 shows the upper and lower bounds applied for building interfaces; in this work we are only interested in lower bounds on bandwidth and upper bounds on workload, which are sufficient to guarantee schedulability among components. In practice, finding the arrival function may be easy as usually the work-arrivals follow the Poisson distribution function for which there is a closed form expression. However, we can also model some other CDFs with the numerical-approximation for count models of the CDF or based on simulations, see [START_REF] Khan | Aperiodic traffic in response time analyses with adjustable safety level[END_REF]; again this is not the goal of the current work. Component Interfaces To model interfaces of WSN components we make use of an approach similar to the real-time calculus [START_REF] Thiele | Real-time calculus for scheduling hard real-time systems[END_REF] and to the assume/guarantee interfaces [START_REF] Wandeler | Interface-based design of real-time systems with hierarchical scheduling[END_REF] tailored towards guarantees on the bandwidth-availability, and requests. An interface of a generic network component has input and output variables related to event streams (the arrivals) and resource availability (the services). We refer our interface model to Probabilistic WSNs Interface (PWI) with the intention of describing component behaviors in terms of curves and their threshold bounds, outlining its probabilistic nature and behavior.

nd i,j C i C β A i,j , Ω i,j β A i , Ω i α i,j , Ω i,j β G i , Λ i β G i,j , Ω i,j
In our framework, the input interface α i,j , Ω i,j (the node's bounded sensing rate) describes the j-th ED component nd i,j of the i-th cluster. The input event/packet stream translates6 into a resource demand β A i,j , Ω i,j ; which defines the probabilistic bandwidth requested from the coordinator, assumed to transmit properly the messages resulting from the event sensing. This is the upper bounded bandwidth requested by a node from its coordinator, based on the nodes constraints such as QoS, deadline, end-to-end delay etc. The ED receives the bandwidth from its parent as β G i,j , Λ i,j , which is the lower bound bandwidth guaranteed by the coordinator. The tuple ( α i,j , Ω i,j , β A i,j , Ω i,j , β G i,j , Λ i,j ) forms the PWI of an ED nd i,j . Similarly, the component link between the coordinator C i and its parent coordinator is given by β A i , Ω i and β G i , Λ i (in terms of bandwidth). These are the assumed (i.e. requested from parent coordinator) and the guaranteed bandwidth (i.e. by the parent coordinator) respectively. Similarly, the coordinator C i 's PWI formed with child nodes is given by β G i , Λ i and β G i,j , Λ i,j . The Figure 3 depicts the network components, their interfaces and their composition. Moreover, ED components may also have outputs representing residual resource ( β ′ i,j , Λ ′ i,j ) and output-arrivals ( α ′ i,j , Ω ′ i,j ), see [START_REF] Santinelli | A component-based framework for modeling and analysing probabilistic real-time systems[END_REF]. We do not discuss the residual resource and output-arrivals in this paper as the analysis here is resource-oriented.

WSN Modeling

With a probabilistic model for WSN components, the bandwidth assignment problem become a probabilistic problem. The bandwidth is allocated in accordance to the requests from the nodes. Each node communicates the tuple β A , Ω to its parent node; that is the probability of demand/request is bound by Ω. In other words it is going to ask for a service greater than β A in 1 -Ω percent of the cases.

Bandwidth Assignment

With the probabilistic model, we can guarantee that a coordinator will provide β G in Λ percent of the cases. For the remaining cases the resource provisioning is less than β G with probability of 1 -Λ. The coordinator then, computes the resource share for the nodes it manages (lower layer nodes) by normalizing the available resource (allocated to the coordinator by its parent coordinator) with respect to the resource demand from its nodes. The i-th coordinator of a cluster asks for the resource from its parent node; which guarantees a resource β G i , Λ i . The resource demand of C i β A i , Ω i (requested by cluster coordinator C i from its parent), comes form the combination of the resource requests of the nodes underneath,

β A i = Ni j=1 β A i,j
with probability Ω i = min i {Ω i,j }. Finally, the resource that C i provides to its ED nodes nd i,j β G i,j , Λ i,j is

β G i,j = β G i β A i β A i,j Λ i,j = Λ i -Ω i,j -Λ i • Ω i,j (5)
where Λ i,j is the probabilistic of the interface.

For the bandwidth allocation among the coordinators; the parent coordinators normalize the available resource using the accumulated resource requests β A i of its children cluster C i s. The normalized resource is then distributed among children cluster nodes with a probabilistic guarantee computed using the probability guarantee of the resource-request, and the probabilistic guarantee of the coordinators resource availability. The resource guarantee for cluster coordinators follows the same principle as in Equation ( 5) where instead of nodes nd i,j is applied to the coordinators C i .

Bandwidth based Component Composability

The probabilistic model allows us to a define flexible relationship among the curves; consequently, this results in the notion of probabilistic composability. Definition 3 ("Greater than or Equal to "( ). We define the greater than or equal to operator 7 ( ) over two probabilistic curves ω, Ω and λ, Λ with ω and λ the curves and Ω and Λ their respective bounding probabilities, as ω, Ω λ, Λ ⇐⇒ ω ≥ λ ∧ Ω ≥ Λ

The following lemma provides the probabilistic guarantees for the composition of two WSN components; wherein one component is assuming a service and the second component is guaranteeing some service. The Lemma 1 is the probabilistic bound on the service that a component will offer to its workload after the composition, i.e. after a cluster coordinator guarantees some service.

Lemma 1 (Resource Reliability). Given two probabilistic curves β A , Ω and β G , Λ respectively an upper and lower bounding probabilistic curves, then min{1, Ω + Λ -Ω • Λ} is the service-reliability probability.

Proof. The output response of the component depends on two inputs, which are β A and β G . We know that probability of β A being larger is Ω and the probability of β G being smaller is Λ. Therefore, probability of the components service reliability depends on, P [β A ∧ β G ]. This probability for the independently distributed random variables can be found as

P [β A ] + P [β G ] -P [β A ] • P [β G ] = Ω + Λ -Ω • Λ,
and since probability can never be larger than one, we have min{1, Ω + Λ -Ω • Λ}.

With this premise, it is possible to define the probabilistic composability for WNSs.

Theorem 1 (Probabilistic Composability). Given two components i and j and two probabilistic curves β G i , Ω i and β A j , Λ j being respectively the probabilistic lower bound to the resource provisioning (the guarantee) of the i-th component and the probabilistic upper bound to the resource demand (the assumption) of the j-th component; then i and j are composable with a probability p iff

β G i , Ω i β A j , Λ j ∧ p ≤ min{1, Ω i + Λ j -Ω i • Λ j } (6)
Proof. The theorem follows as a consequence of Lemma 1

4 Application to IEEE 802.15.4

4.1 CSMA/CA using IEEE 802.15.4 superframe structure

In the beacon-enabled mode of IEEE 802.15.4, a coordinator in a Personal Area Network (PAN) periodically sends beacon frames to synchronize the associated nodes. The standard defines a superframe structure between two successive beacon transmissions. The superframe structure is divided into an active portion and a low power inactive portion. Figure 4 shows an example of the superframe structure. The superframe structure is specified by two values; the superframe duration (SD) which defines the active portion and the beacon interval (BI) which defines the interval between two consecutive beacons. The SD and BI periods are defined using two parameters; the superframe order (SO) and the beacon order (BO), respectively. Equation [START_REF] She | Deterministic worst-case performance analysis for wireless sensor networks[END_REF] gives the definition of SD and BI as a function of SO and BO, respectively. In this equation, aBaseSuperframeDuration is a constant value defined by the standard which is equal to 15.36ms

BI = aBaseSuperf rameDuration • 2 BO SD = aBaseSuperf rameDuration • 2 SO (7) 1 ≤ SO ≤ BO ≤ 14
The standard supports the cluster tree topology wherein the coordinators form a multi-hop tree. The root coordinator is called PAN coordinator. Every coordinator provides a synchronization services (through beacon transmissions) to other devices or other coordinators. Therefore, to avoid beacon collision a beacon scheduling scheme has to be used among coordinators. The IEEE 802.15.4 standard imposes that BO and SO have to be equal for all superframes on a PAN and suggests the scheduling of beacons so that active periods of neighbor coordinators and two-hop neighbor coordinators do not overlap. Figure 4 gives an example of beacon scheduling of three coordinators.

Each child forwards data only during its parents active period. In particular, to transmit data to the PAN coordinator, the child coordinators and the EDs must transmit their data during the active period of the PAN coordinator. Therefore, the maximum bandwidth to be assigned by the PAN coordinator, can be the maximum number of messages (noted MaxNbMsgs hereafter) that are to be transmitted to the PAN coordinator during the SD of every BI. We have done some preliminary simulations to determine suitable values of MaxNbMsgs. Of course, it depends on the duration of the SD period. Table 1 shows the values of MaxNbMsgs as a function of SO.

Admission control manager

The ACM, Figure 5, is a component of the MAC layer in every node of WSN. The ACM ensures that the bandwidth assignment policy is respected. The network has two working stages; a setup stage and a normal working stage. The setup stage works as follows. First, the PAN coordinator gathers service requests from all the nodes of the network. Then, it runs the bandwidth assignment algorithm 8 . Finally, the bandwidth determined by the algorithm for each node is assigned as the guaranteed service to cluster coordinators and EDs. During the normal working stage the ACM becomes active. We note that the PAN coordinator can decide to re-execute the setup stage sometime later. The decision can be based on the reception of new service requests or changes in the network; like node mobility or death.

The ACM constantly monitors all the packets received by the MAC layer (from upper and lower layers) and decides which ones are accepted or dropped, exercising admission control tests. The decision to drop packets is based on the assigned guaranteed service to the node, β G , Λ . However, the ACM of the ED and the ACM of the cluster coordinators work in a different manner. For EDs, it monitors only the packets received from the application layer and ensures that the number of forwarded packets to the parent node does not exceed β G . The Λ G is not taken into consideration in the packet dropping process of the EDs ACM. For cluster coordinators, the ACM monitors the packets received from children EDs and the application layer. The dropping of the packets received from children EDs is based on Λ G of the transmitting ED; while as it is based on the β G of the cluster coordinator for the packets it receives from its application layer. We note that the ACM of the PAN coordinator monitor all received packets but does not drop any of them in case it is the destination. We now can define a metric based on the packets dropped by the ACM, called as packet-loss ration as: Definition 4 (Packet-Loss-Ratio(PLR)). The packet-lost-ratio (PLR) for i th node is defined as l i (t) = Li(t) Ri(t) , where L i (t) gives the packets lost in an interval (0, t]. This is be achieved by keeping two counter in an ACM; that is to count L and R. This gives us an estimate of PLR at each event (i.e. packet dropped or packet arrived). Therefore, intuitively for a given bandwidth assignment higher PLR means lower energy usage by a node; since the unnecessary transmission due to limited bandwidth were avoided. This is because energy-saving due to dropping dominates the total energy consumption; simply put ACM tries to ensure than energy is only consumed for potentially successful transmissions. Although, higher PLR may seem to indicate a noticeable performance issue (QoS); however, it is important to note that the packets are dropped in an ACM based on the agreed level loss indicated by the probabilistic bounds. Moreover, higher PLR does not always indicate a problem; for example when high PLR at the PAN node is acceptable (as it may be within required QoS level). Let us take the example of an ED which has obtained β G = 2 packets/superframe, Λ G = 50% from its parent, named C F . The ACM of ED is responsible for ensuring that no more than two packets are transmitted to its parent C F every superframe. For example, if the application generates four packets during the actual superframe, two of them will be dropped. The ACM of the C F transmits the two received packets to its own parent only 50% of the time; for other 50% of the time the ACM can transmit zero or one packet to its own parent (i.e. 50% of the time, it may drop one or both packets).

Simulation Setup

As a test case we consider the network shown in Figure 2. It is composed of one PAN coordinator (Network Coordinator), three cluster coordinators and 30 EDs, five EDs attached to the first coordinator, 15 EDs attached to the second coordinator and 10 EDs attached to the third coordinator.

All coordinators, including the PAN coordinator, have a non overlapping active periods. The network uses the ZigBee [START_REF]Zigbee Specification Document[END_REF] tree routing protocol. The simulation is performed using OPNET simulator [22], and Table 2 summarizes all the simulation parameters. We performed two sets of simulations. The first one illustrates the behavior of the ACM. The second one compares the results of our bandwidth allocation scheme using the ACM against the original IEEE 802.15.4 standard. The simulations and the results are described next.

The behavior of the ACM

We fix the value of SO to three and varied the EDs service request. The duration of the active period is constant and β G = 120 packets/superframe, see Table 1. We also note that coordinators do not transmit messages. We report the bandwidth assigned by the PAN coordinator to each ED; the total number of packets sent; the total number of packets received by the PAN coordinator (which is the destination); the total number of packets dropped by the ACM; and the total number of packets dropped due to CSMA/CA failure or collision.

Table 3 illustrates the obtained results. We observe that if the total service request is less than the guaranteed bandwidth, cases of α A i = 0.5(pkts/s), Ω A i = 0.8 and α A i = 1(pkts/s), Ω A i = 0.8 , the ACM drops a few packets. In fact, the ACMs of EDs do not drop any packets, since the guaranteed service is higher than the requested service and the packet inter-arrival distribution is constant. However, since the guaranteed probability is not equal to one, some of the packets will be dropped by the ACMs of cluster heads. When the total service request exceeds the guaranteed bandwidth the ACMs of EDs drop the packets to respect the allocated bandwidth; that is α A i = 2(pkts/s), Ω A i = 0.8 . We also notice that the admission of packets by ACM for transmission does not guarantee the packets successful transmission. This is illustrated by the number of dropped packets due to CSMA/CA failure or collisions. In fact, multiple nodes may try to transmit accepted packets at the same time which results in a collision. We recall also that transmissions are not acknowledged and the duty cycle is equal to 3.125%. Nevertheless, the use of ACM minimizes the number of dropped packets due to collision or CSMA/CA failure as we will see in the next subsection.

Comparison with the original standard

In this simulation set, we fixed the service request of EDs and varied the SO. We measure the number of packets received by the PAN coordinator, the number of packets dropped either due to CSMA/CA failure or due to collisions, and the average end-to-end delay of the received packets. We compare the results of our bandwidth allocation scheme using the ACM against the original IEEE 802.15.4 standard (without the use of any bandwidth allocation scheme and control). For comparison we have set the service request probability (Ω A i ) of each ED to 1.

Figure 6 shows the simulation results. Figure 6(a) corresponds to the number of the received packets by the PAN coordinator. We can notice that when using ACM and bandwidth-allocation: the number of received packets increases as SO increases. In fact, increasing SO means increasing the guaranteed bandwidth as pointed out in Table 1; therefore, each ED receives higher allocated bandwidth. We also observe that we obtain similar results for the number of received packets with and without the use of the ACM. However, both versions (with and without the ACM) do not follow the same path to obtain these similar results. Indeed, without using the ACM, dropped packets are caused by collisions or CSMA/CA failure. However, when the ACM is used, packets are dropped mainly by the ACM (though some packets are dropped by collisions or CSMA/CA failures). Figure 6(b) and figure 6(c) illustrates this result. The use of the ACM, therefore, avoids useless transmissions; as the chance of failure/collision increases due to the limited bandwidth. The ACM ensures that the allocated bandwidth is respected by dropping the packets before they attempt their transmission; which can reduce the number of collisions. Consequently, this decreases the end-to-end delay. Figure 6(d) shows the end-to-end delays of the received packets with and without the use of the ACM. Indeed, the results show that the use of the ACM decreases the end-to-end delay of the received packets. Another direct consequence of actively dropping overloaded packets by ACM instead of leaving them dropped because of too collisions is that there will be less retransmissions, so less energy waste.

Conclusions

In this paper we proposed an approach which applies a component based and hierarchical cluster-tree topologies for WSNs by modeling and analyzing such distributed systems in a probabilistic approach. We also proposed a MAC level Admission Control Manager (ACM) which does bandwidth assignment and control for the required performance guarantees (QoS) in a probabilistic manner using the probabilistic bounds. Our framework offers flexibility of a probabilistic approach, therefore, we are able to provide probabilistic guarantees for the system functionality, or degradation; in harsh environments and complex systems.

The probabilistic bandwidth assignment approach has been successfully implemented in IEEE 802.15.4 through the admission control manager. The simulation results showed that the ACM drops packets to respect the bandwidth assigned by the PAN coordinator (to every device in the network). Consequently, we were able to decrease the number of useless transmissions while providing the same performance; that is the number of received packets. Which ensures that the energy is conserved by avoiding useless transmissions. Moreover, we were able to ameliorate the end-to-end delay. Nevertheless, because of the inherent behavior of CSMA/CA, the dropping of packets due to the collisions cannot be totally prevented by the use of ACM. Moreover, we can apply our framework to a contention-free access protocol with more success.

In the future we intend to improve our framework by: i) applying it to a contention-free access protocol; ii) applying the analysis to different network topologies; iii) inferring network characteristics such as the lifetime; and iv) checking its applicability to real WSN platforms/testbeds in the dynamic network conditions and testing the adaptivity of the probabilistic approach.

Fig. 1 .

 1 Fig. 1. A model of a component with its interface abstracting a part of a system.

Fig. 2 .

 2 Fig. 2. WSN hierarchical architecture with end devices, control coordinators and clusters elements. Our simulation applies this architecture with N1 = 5, N2 = 15 and N3 = 10.

Fig. 3 .

 3 Fig. 3. WSN Component composition.

Fig. 4 .

 4 Fig. 4. The superframe structure of the PAN coordinator and an example of beacon scheduling of three coordinators.

Fig. 5 .

 5 Fig. 5. System view and ACM.

Fig. 6 .

 6 Fig. 6. Results of the second simulation set. "Without ACM" corresponds to the standard IEEE 802.15.4.

Table 1 .

 1 Maximum number of messages that can be transmitted to the PAN coordinator

	SO	2 3 4
	MaxNbMsgs 80 120 220

Table 2 .

 2 General simulation parameters

		Application parameters
	Packet length		100 bits
	Service Request α A i (p/s), Ω A i	first set of sim	variable
	Service Request α A i (p/s), Ω A i	second set of sim	2 p/s, 1
	Inter-arrival distribution	constant
	Destination		PAN coordinator
		PG		0.9
		Network parameters
	Cm, Rm and Lm	20, 3 and 2
		MAC parameters
	superframe		BO SO	8 2, 3 and 4
	Beacon Start Time	C1, C2 and C3	0.49s, 1s and 1.7s
	CSMA/CA	Minimum Back-off Exponent Maximum Back-off Number	3 4
	Acknowledgment	disabled
		PHY parameters
		Data rate		250 kb/s
		Ploss		0.1
		Other parameters
	Simulated duration	900s
	Energy model		Micaz

Table 3 .

 3 Results of the first set of simulations. pkts, sf, and Nb designate packets, superframe, and number respectively.

	α A i (pkts/s), Ω A i	0.5, 0.8	1, 0.8	2, 0.8
	β A i (pkts/sf ), Ω A i	1.96, 0.8 3.96, 0.8 7.86, 0.8
	Total service requests	58.98	117.96	235.92
	β G (pkts/sf ), Λ G	4, 0.98	4, 0.98	4, 0.98
	Nb of created pkts	13500	26970	53940
	Nb of received pkts	5462	8252	8337
	Nb of dropped pkts by the ACM	114	243	26694
	Nb of dropped pkts due to
	CSMA/CA failure or collision	7924	18475	18909

This work has been partially supported by the ANR Quasimodo project under grant ANR 2010 INTB 0206 01.

The assumptions are merely considered for the ease of explanation. However, the approach is valid for other application scenarios and network architectures.

That is αi and β A i have same distribution function.