Alice Le Brigant

Stéphane Puechmorel

Quantization and clustering on Riemannian manifolds with an application to air traffic analysis

Keywords: optimal quantization, clustering, Riemannian geometry, air traffic analysis

The goal of quantization is to find the best approximation of a probability distribution by a discrete measure with finite support. When dealing with empirical distributions, this boils down to finding the best summary of the data by a smaller number of points, and automatically yields a K-means-type clustering. In this paper, we introduce Competitive Learning Riemannian Quantization (CLRQ), an online quantization algorithm that applies when the data does not belong to a vector space, but rather a Riemannian manifold. It can be seen as a density approximation procedure as well as a clustering method. Compared to many clustering algorihtms, it requires few distance computations, which is particularly computationally advantageous in the manifold setting. We prove its convergence and show simulated examples on the sphere and the hyperbolic plane. We also provide an application to real data by using CLRQ to create summaries of images of covariance matrices estimated from air traffic images. These summaries are representative of the air traffic complexity and yield clusterings of the airspaces into zones that are homogeneous with respect to that criterion. They can then be compared using discrete optimal transport and be further used as inputs of a machine learning algorithm or as indexes in a traffic database.

Introduction

Most of the statistical tools developed so far are dedicated to data belonging to vector spaces, since it is the most convenient setting for algorithms as well as theoretical derivations. However, when dealing with real world applications, such a framework may not fit with the structure of the data. It is obviously the case for geostatistics over world-sized datasets, but it is also true in many other fields: shapes in computer vision, diffusion tensor images in neuroimaging, signals in radar processing do not belong to a Euclidean space, but rather to a differentiable manifold. Riemannian geometry provides a convenient framework to deal with such objects. It allows a straightforward generalization of basic statistical notions such as means and medians [START_REF] Arnaudon | Medians and means in riemannian geometry: existence, uniqueness and computation[END_REF][START_REF] Fréchet | Les éléments aléatoires de nature quelconque dans un espace distancié[END_REF][START_REF] Karcher | Riemannian center of mass and mollifier smoothing[END_REF], covariance [START_REF] Pennec | Intrinsic statistics on riemannian manifolds: Basic tools for geometric measurements[END_REF], Gaussian distributions [START_REF] Said | Riemannian gaussian distributions on the space of symmetric positive definite matrices[END_REF], and of usual linear operations such as principal component analysis [START_REF] Fletcher | Principal geodesic analysis for the study of nonlinear statistics of shape[END_REF][START_REF] Sommer | Manifold valued statistics, exact principal geodesic analysis and the effect of linear approximations[END_REF]. The use of these statistical tools has met a growing interest in various fields, including shape analysis [START_REF] Kendall | Shape manifolds, procrustean metrics, and complex projective spaces[END_REF], computational anatomy [START_REF] Fletcher | Principal geodesic analysis for the study of nonlinear statistics of shape[END_REF], medical imaging [START_REF] Fletcher | Riemannian geometry for the statistical analysis of diffusion tensor data[END_REF], probability theory [START_REF] Bigot | Geodesic pca in the wasserstein space by convex pca[END_REF], and radar signal processing [START_REF] Arnaudon | Riemannian medians and means with applications to radar signal processing[END_REF][START_REF] Brigant | Computing distances and geodesics between manifold-valued curves in the srv framework[END_REF].

In air traffic management (ATM), a major concern is the ability to infer an estimation of the complexity as perceived by a controller from the knowledge of aircraft trajectories in a given airspace, as depicted in Figure 1. Many interdependent factors are involved in the cognitive process of a human controller, making the problem extremely difficult to solve, if even possible. However, there is a consensus among the experts on the importance of traffic disorder. As detailed in [START_REF] Mykoniatis | A new representation of air traffic data adapted to complexity assessment[END_REF], a good way to estimate traffic disorder is to assume the spatial distribution of the aircraft velocities to be Gaussian and use the covariance function as an indicator for local complexity. This model requires to estimate a covariance matrix at each sample point of the airspace, and therefore yields a mapping from the plane to the space of symmetric, positive definite (SPD) matrices. Such a mapping will be referred to as an (SPD) image in the sequel. Working directly with the SPD images is an extremely computationally expensive task, that is unrealistic in practice. Moreover, the information provided is highly redundant, making it cumbersome to use in statistical analysis. To cope with this problem, we seek to produce summaries of the images. To do so, we model an SPD image as a realization of a random field with values in the space of SPD matrices. The values collected at each sample point in the image define an empirical probability distribution, supported on the space of SPD matrices. We propose to produce a summary of this empirical distribution using optimal quantization.

Optimal quantization is concerned with finding the best approximation, in the sense of the Wasserstein distance, of a probability distribution µ by a discrete measure μn with finite support |supp μn | ≤ n (see [START_REF] Graf | Foundations of quantization for probability distributions[END_REF] or the survey paper [START_REF] Pagès | Introduction to vector quantization and its applications for numerics[END_REF], and references therin). When dealing with empirical distributions, this boils down to finding the best summary of the data by a smaller number of points. In the same setting, quantization naturally yields a clustering of the data points, which coincides with the solution given by the K-means algorithm. In our application, probability measures are supported on the space of SPD matrices, and therefore points actually refer to SPD matrices. The main challenge lies in the fact that SPD matrices do not form a vector space, but rather a Riemannian manifold [START_REF] Pennec | A riemannian framework for tensor computing[END_REF], whereas most work on optimal quantization is suited for vector data [START_REF] Bouchitté | Asymptotic analysis of a class of optimal location problems[END_REF][START_REF] Pagès | Optimal quantization methods and applications to numerical problems in finance[END_REF], or functional data [START_REF] Biau | On the performance of clustering in hilbert spaces[END_REF][START_REF] Pages | Quadratic optimal functional quantization of stochastic processes and numerical applications[END_REF]. Nonetheless, the case of probability distributions on Riemannian manifolds has recently received attention [START_REF] Iacobelli | Asymptotic quantization for probability measures on riemannian manifolds[END_REF][START_REF] Kloeckner | Approximation by finitely supported measures[END_REF]. In particular, the asymptotic behavior of the quantization error, i.e., the evolution of the error made by approximating µ by μn as n → ∞, was studied for the manifold case in [START_REF] Iacobelli | Asymptotic quantization for probability measures on riemannian manifolds[END_REF]. However, to the best of our knowledge, no numerical schemes have yet been introduced in this setting.

In this work, we introduce Competitive Learning Riemannian Quantization (CLRQ), a Riemannian counterpart of Competitive Learning Vector Quantization (see for example [START_REF] Pagès | Introduction to vector quantization and its applications for numerics[END_REF]). It is a gradient descent algorithm that computes the best approximation μn of a probability measure µ over a Riemannian manifold using observations sampled from µ. For empirical distributions, this allows to summarize a manifold-valued dataset of size N by a smaller number n N of points, which additionally correspond to the centers of a K-means-type clustering. Recall that the regular Riemannian K-means algorithm computes at each step a partition of the data space (of the size of the desired number of clusters) and computes the geometric (Karcher) mean of each cell of the partition [START_REF] Turaga | Statistical computations on grassmann and stiefel manifolds for image and video-based recognition[END_REF]. These two steps are usually very computationally expensive in the Riemannian setting, where the geodesic distance doesn't always have a closed form and the mean can only be obtained by an iterative procedure. Other clustering procedures relying on the estimation of distance-based kernels or graphs, such as mean-shift [START_REF] Subbarao | Nonlinear mean shift over riemannian manifolds[END_REF], or persistence-based clustering [START_REF] Chazal | Persistence-based clustering in riemannian manifolds[END_REF], require at least the computation of the matrix of pairwise-distances between the data points, if not the distances from the current (new) point to all the data points at each step. It is also the case of kernel-based approaches which first transport the data into a Reproducing Kernel Hilbert Space (RKHS) before applying a clustering procedure, such as K-means. In [START_REF] Zhao | Efficient clustering on riemannian manifolds: A kernelised random projection approach[END_REF], the authors propose to accelerate the procedure through random projections, which require to compute the distance matrix of a small subset of points instead of the whole dataset. In contrast, CLRQ is an on-line algorithm (i.e., processing one data point at a time) and requires at each step as many distance computations as the number of desired clusters, which is usually very small compared to the number of data. This leads us to argue that CLRQ is well-suited for large datasets.

After proving the convergence of CLRQ under the relevant assumptions, we use it to perform air traffic complexity analysis as presented above. Applying the algorithm to the empirical distribution of an SPD image yields two desirable results: (1) a summary of the image through μn , and (2) a clustering of the image, and thereby of the corresponding airspace, into different zones homogeneous with respect to complexity. The latter means that to each point of the image is assigned a level of complexity according to the class it belongs to, and the former, that different traffic images can easily be compared through the comparison of their summaries, using the Wasserstein distance. This is an interesting prospect, since it allows for the indexation of air traffic situations in view of creating automatized decision-making tools to help air traffic controllers. Naturally, CLRQ can conceivably be applied to many other application-driven problems, requiring either the construction of summaries or the clustering of geometric data. A natural example of such data in positive curvature is given by GPS positions on the Earth (spherical data); in negative curvature, one can think of Gaussian distributions, which can be parameterized by the hyperbolic half-plane in the univariate case, and by higher dimensional Riemannian symmetric spaces in the multivariate case [START_REF] Calvo | An explicit solution of information geodesic equations for the multivariate normal model[END_REF][START_REF] Lovrić | Multivariate normal distributions parametrized as a riemannian symmetric space[END_REF].

The paper is organized as follows. In Section 2, we introduce the Riemannian setting and give an example of a statistical object on a manifold (the Riemannian center of mass), before introducing the context of optimal quantization. In Section 3, we present the CLRQ algorithm and show its convergence. After showing some simple illustrations on the circle, the sphere and the hyperbolic half-plane in Section 4, we present our main application in air traffic management in Section 5.

Mathematical setup

Notions of Riemannian geometry

Let us begin by introducing some notations and reminding some basic notions of Riemannian geometry. We consider a differentiable manifold M of dimension d equipped with a Riemannian metric, i.e., a smoothly varying inner product •, • x defined on each tangent space T x M at x ∈ M. Recall that T x M is a linear approximation of M at point x, and contains all tangent vectors to M at x. The subscript x in the metric will often be omitted and the norm associated to the Riemannian metric •, • will be denoted by • . Vector fields on M are mappings X : M → T M that associate to each point x a tangent vector X x ∈ T x M. It is possible to take the derivative of a vector field with respect to another using an affine connection, that is, a functional (X, Y) → ∇ X Y that acts on pairs of vector fields according to the following rules

∇ f X+Y Z = f ∇ X Z + ∇ Y Z, ∇ X ( f Y) = X( f )Y + f ∇ X Y, ∇ X (Y + Z) = ∇ X Y + ∇ X Z.
If a vector field is defined along a curve γ(t), i.e., V(t) := X γ(t) , then its covariant derivative along γ is denoted by DV dt (t) := ∇ γ (t) V. There is only one affine connection that is symmetric, meaning XY -Y X = ∇ X Y -∇ Y X, and compatible with the Riemannian metric, that is

d dt U, V = DU dt , V + U, DV
dt , for any vector fields U, V along a curve γ. It is called the Levi-Civita connection and we will always take covariant derivatives with respect to that connection. The geodesics of M are the curves γ satisfying the relation Dγ dt = 0, which implies that their speed has constant norm γ (t) = cst. They are also the local minimizers of the arc length functional l:

l : γ → 1 0 γ (t) dt
where in the previous expression curves are assumed, without loss of generality, to be defined over the interval [0, 1]. The exponential map at point x is the mapping, denoted by exp x , that maps a tangent vector v of an open ball

B 0 ⊂ T x M centered in 0 ∈ T x M to the endpoint γ(1) of the geodesic γ : [0, 1] → M verifying γ(0) = x, γ (0) = v, exp x (v) = γ(1).
Intuitively, the exponential map moves the point x along the geodesic starting from x at speed v and stops after covering the length v . The image exp x (B 0 ) is called a normal neighborhood of x. Conversely, the inverse of the exponential map associates to any y in this normal neighborhood the vector exp -1 x (y) that sends x on y through the exponential map, and will be denoted by the more intuitive notation -→ xy. We assume that M is complete, i.e., that the exponential map at x is defined on the whole tangent space T x M. By the Hopf-Rinow theorem, we know that M is also geodesically complete, that is, any two points x, y ∈ M can be joined by a geodesic of shortest length. This minimal length defines the geodesic distance between x and y, denoted in the sequel by d(x, y). For further details, we refer the reader to standard textbooks such as [START_REF] Carmo | Riemannian geometry[END_REF] or [START_REF] Jost | Riemannian geometry and geometric analysis[END_REF].

The Riemannian center of mass

Before introducing optimal quantization, let us briefly give an example of generalization of a fundamental statistical notion to the manifold setting: the Riemannian center of mass of a probability distribution, also called the Fréchet mean [START_REF] Fréchet | Les éléments aléatoires de nature quelconque dans un espace distancié[END_REF]. Consider a probability measure µ on M with density with respect the the Riemannian volume form and a compact support K = suppµ. As a compact set, K is contained in a geodesic ball B(a, R) = {x ∈ M, d(a, x) < R}, and we assume that this ball is convex, i.e., that for all points x, y ∈ B(a, R), there exists a unique minimizing geodesic from x to y (sufficiently small geodesic balls are convex [START_REF] Cheeger | Comparison theorems in Riemannian geometry[END_REF]Theorem 5.14]). Let X be a random variable of law µ. The Fréchet mean x of X generalizes the minimum-square-error property of the Euclidean mean

x = E µ (X) = arg min a∈M M d(x, a) 2 µ(dx). (1)
It can be characterized as the point where the gradient of the above functional vanishes, which gives

M -→ xxµ(dx) = 0. ( 2
)
This is due to the fact that the gradient of the functional f : a → d(x, a) 2 is given by grad a f = -2 -→ ax (see Lemma 1 in the Appendix), and to the assumptions we made on M and supp µ. The same assumptions guarantee that x exists and is unique [START_REF] Pennec | Intrinsic statistics on riemannian manifolds: Basic tools for geometric measurements[END_REF]. When µ is an empirical probability measure equally distributed on N data points x 1 , . . . , x N , the mean can be computed using a gradient descent algorithm, where the update is given by

x ← exp x        1 N N i=1 -→ xx i        .
This algorithm is often denoted by the Karcher flow.

Optimal Riemannian quantization

Optimal quantization addresses the problem of approximating a random variable X with distribution µ by a simplified version q(X) where q : M → Γ ⊂ M is a measurable function with an image Γ of cardinal at most n. The function q is called an n-quantizer and is chosen to minimize the L p criteria, called the distortion function, E µ {d(X, q(X)) p } .

(

) 3
Since q takes only a finite number of values, the distribution of q(X) will be a finite sum of point measures. It is well known that, since any n-quantizer q of image Γ ⊂ M verifies for all x ∈ M, d(x, q(x)) ≥ inf a∈Γ d(x, a), with equality if and only if q(x) = arg min a∈Γ d(x, a), the optimal quantizer is the projection to the nearest neighbor of Γ. Moreover, if |Γ| < n and |supp µ| ≥ n, one easily checks that q can always be improved, in the sense of criteria (3), by adding an element to its image. This means that an optimal n-quantizer has an image of exactly |Γ| = n points. Therefore, the optimal n-quantizer for criteria (3) is of the form q Γ : M → Γ = {a 1 , . . . , a n }, where the a i 's are pairwise distinct, and

q Γ (•) = n i=1 a i 1 C i (Γ) (•), where C i (Γ) = {x ∈ M, d(x, a i ) ≤ d(x, a j ) ∀ j i}.
The ties in C i (Γ) are arbitrarily broken. The set C i (Γ) is the i th Voronoi cell associated to Γ and the union of all these cells form the Voronoi diagram. The quantization problem is therefore equivalent to the approximation of supp µ by an n-tuple (a 1 , . . . , a n ) ∈ M n minimizing the cost function F n,p :

M n → R + , F n,p (a 1 , . . . , a n ) = E µ min 1≤i≤n d(X, a i ) p = M min 1≤i≤n d(x, a i ) p µ(dx). (4)
This cost function is obtained by evaluating the distortion function (3) at q = q Γ . Notice that if we seek to approximate µ by a single point a ∈ M (i.e., n = 1) with respect to an L 2 criterion (p = 2), we retrieve the definition (1) of the Riemannian center of mass.

Finally, there is a third way, in addition to (3) and (4), of expressing the quantization problem: it is also equivalent to the approximation of the measure µ by the closest discrete measure μn supported by n points, with respect to the Wasserstein distance of order p W p (µ, μn ) = inf P d(u, v) p dP(u, v).

(

) 5
Here the infimum is taken over all measures P on M × M with marginals µ and μn . One can construct an optimal discrete approximation μn (i.e., a minimizer of ( 5)) from an optimal n-tuple (a 1 , . . . , a n ) (i.e., a minimizer of ( 4)), and vice-versa, using

μn = n i=1 µ {C i (Γ)} δ a i , (6)
and then we have

F n,p (a 1 , . . . , a n ) = W p (µ, μn ).
This is well known for the vector case [START_REF] Graf | Foundations of quantization for probability distributions[END_REF] and applies verbatim to measures on manifolds. In the sequel, we will focus on the second formulation (4) of the quantization problem. The first question that arises is the existence of a minimizer of ( 4). Since we have assumed that µ has compact support, this existence is easily obtained. Proposition 1. Let M be a complete Riemannian manifold and µ a probability distribution on M with density and a compact support. Then the distortion function F n,p is continuous and admits a minimizer.

Proof. Just as in the vector case [31, Proposition 1.1], for any x ∈ M, the function

M n → R + , α = (a 1 , . . . , a n ) → min 1≤i≤n d(x, a i ) is 1-Lipschitz for the distance d (α, β) := max 1≤i≤n d(a i , b i ), where β = (b 1 , . . . , b n ).
Therefore it is continuous, and so is its p th power. Since K = supp µ is compact, for all α ∈ M n and all β in a neighborhood

B(a 1 , r 1 ) × . . . × B(a n , r n ) of α, we have ∀x ∈ K, min 1≤i≤n d(x, b i ) p ≤ min 1≤i≤n sup y∈K d(y, a i ) + r i p < ∞.
So by dominated convergence, F n,p is continuous. Recall that as a compact set, K is contained in a geodesic ball B(a, R). If α = (a 1 , . . . , a n ) ∈ M n is such that d(a, a i ) > 2R for at least one a i , then for all x ∈ K, d(x, a i ) ≥ d(a, a i ) -R > R, and so the same n-tuple where a replaces a i is a better candidate to minimize F n,p . We can therefore limit our search to B(a, 2R), which is a closed and bounded subset of the complete manifold M, and thus compact. The continuous function F n,p reaches a minimum on this compact, which is an absolute minimum.

The elements of a minimizer α = (a 1 , . . . , a n ) are called optimal n-centers of µ. The minimizer α is in general not unique, first of all because any permutation of α is still a minimizer, and secondly because any symmetry of µ, if it exists, will transform α into another minimizer of F n,p . For example, any rotation of the optimal n-centers of the uniform distribution on the sphere preserves optimality.

The second question that comes naturally is: how does the error one makes by approximating µ by μn (as given by ( 6)) evolve when the number n of points grows ? The n th quantization error is defined by

V n,p (µ) = inf (a 1 ,...,a n )∈M n F n,p (a 1 , . . . , a n ) = inf (a 1 ,...,a n )∈M n M min 1≤i≤n d(x, a i ) p µ(dx).
In the vector case, Zador's theorem [START_REF] Graf | Foundations of quantization for probability distributions[END_REF]Theorem 6.2] tells us that it decreases to zero as n -p/d , and that the limit of n p/d V n,p (µ) can be expressed in terms of the p th quantization coefficient, i.e., the limit (which is also an infimum) when µ is the uniform distribution on the unit square of

R d Q p [0, 1] d = lim n≥1 n p/d V n,p U [0, 1] d .
More precisely, if µ has an absolutely continuous part (with respect to the Lebesgue measure) with density h, then the limit of the n th quantization error is equal to

Q p ([0, 1] d ) h d/(d+p)
. Moreover, when µ has no singular part, the asymptotic empirical distribution of the optimal n-centers is proportional to h d/(d+p) .

In the case of a Riemannian manifold M, the moment condition of the flat case generalizes to a condition involving the curvature of M. The following term measures the maximal variation of the exponential map at x ∈ M when restricted to a (d -1)-dimensional sphere S ρ ⊂ T x M of radius ρ

A x (ρ) = sup v∈S ρ ,w∈T v S ρ , w =ρ d v exp x (w) .
The following generalization of Zador's theorem to Riemannian quantization was proposed by Iacobelli [START_REF] Iacobelli | Asymptotic quantization for probability measures on riemannian manifolds[END_REF].

Theorem 1 ([18, Theorem 1.4 and Corollary 1.5]). Let M be a complete Riemannian manifold without boundary, and let µ = h dvol + µ s be a probability measure on M, where dvol denotes the Riemannian volume form and µ s the singular part of µ. Assume there exist x 0 ∈ M and δ > 0 such that

M d(x, x 0 ) p+δ dµ(x) + M A x 0 {d(x, x 0 ) p }dµ(x) < ∞.
Then

lim n→∞ n p/d V n,p (µ) = Q p [0, 1] d h d/(d+p) ,
where • r denotes the L r -norm. In addition, if µ s = 0 and (a 1 , . . . , a n ) are optimal n-centers, then

1 n n i=1 δ a i D -→ λh d/(d+p) dx as n → ∞,
where D → denotes convergence in distribution and λ is the appropriate normalizing constant.

In this work, we are interested in finding numerical schemes to compute the optimal n-centers α = (a 1 , . . . , a n ) in practice from potentially large sets of data. To do so, we will search for the critical points of the distortion function.

Competitive Learning Riemannian Quantization

Differentiability of the distortion function

We assume that the only knowledge that we have of the probability measure µ that we want to approximate is through an online sequence of i.i.d. observations X 1 , X 2 , . . . sampled from µ. A classical algorithm used for quadratic (p = 2) vector quantization is the Competitive Learning Vector Quantization algorithm, a stochastic gradient descent method based on the differentiability of the distortion function F n,2 . We propose here a natural extension of this method to our setting, i.e., a compactly-supported probability measure on a complete Riemannian manifold. It relies on the differentiability of the distortion function.

Proposition 2. Let α = (a 1 , . . . , a n ) ∈ M n be an n-tuple of pairwise distinct components and p > 1. Then F n,p is differentiable and its gradient in α is

grad α F n,p =        -p Ci (α) -→ a i x p-1 -→ a i x -→ a i x µ(dx)        1≤i≤n ∈ T α M n ,
where Ci (α) is the interior of the i th Voronoi cell of α and

-→ xy := exp -1
x (y) denotes the vector that sends x on y through the exponential map. In particular, the gradient of the quadratic distortion function is given by

grad α F n,2 = -2 Ci -→ a i x µ(dx) 1≤i≤n = -2 E µ 1 {X∈ Ci } -→ a i X 1≤i≤n . ( 7
)
Remark 1. The first observation we can make is that optimal n-centers are Riemannian centers of mass of their Voronoi cells, as characterized by [START_REF] Arnaudon | Riemannian medians and means with applications to radar signal processing[END_REF]. Hence, the term n-centers is as appropriate in the Riemannian setting as it is in the vector case. More generally, for any value of p, each a i , i = 1, . . . , n, is the p-mean of its Voronoi cell, i.e., the minimizer of a →

Ci (α) d(x, a) p µ(dx).
Therefore, the optimal n-centers are always contained in the compact support of µ.

Remark 2. The second observation is that the opposite direction of the gradient is, on average, given by the vectors inside the expectation. Competitive learning quantization consists in following this direction at each step k, that is, updating only the center a i corresponding to the Voronoi cell of the new observation X k , in the direction of that new observation. In the Riemannian setting, instead of moving along straight lines, we simply follow geodesics using the exponential map.

Proof. Let α = (a 1 , . . . , a n ) ∈ M n be an n-tuple of pairwise distinct components, w = (w 1 , . . . , w n ) ∈ T α M n a tangent vector to α, and let (-, ) t → α(t) = (a 1 (t), . . . , a n (t)) be a variation of α in the direction given by w, i.e., such that a i (0) = a i et a i (0) = w i for all i = 1, . . . , n. The functional t → min 1≤i≤n d(x, a i (t)) p is differentiable for all x ∪ 1≤i≤n ∂C i (α(t)), that is µ-almost everywhere since µ(∂C i (α(t))) = 0 for all i (this is shown in the Appendix, Lemma 2). Its derivative in t = 0 is given by

d dt t=0 min 1≤i≤n d(x, a i (t)) p = n i=1 1 {x∈ Ci } d dt t=0 d(x, a i (t)) p = - n i=1 1 {x∈ Ci } p 2 d(x, a i (0)) 2 p/2-1 2 ----→ a i (0)x, a i (0) = - n i=1 1 {x∈ Ci } p --→ a i x p-2 --→ a i x, w i .
To go from the second to the third line, we have used the well known property that for any given x, the gradient of the function f : a → d(x, a) 2 is given by grad a f = -2 exp -1 a x = -2 -→ ax (see Lemma 1 in the Appendix). We obtain by Cauchy-Schwarz, since x ∈ B(a, R) and a i ∈ B(a, 2R) (recall that any n-tuple containing a coordinate outside of B(a, 2R) is a worse candidate than the same n-tuple where a replaces a i ),

d dt t=0 min 1≤i≤n d(x, a i (t)) p ≤ p n i=1 --→ a i x p-1 w i ≤ p(3R) p-1 n i=1 w i .
Therefore, by dominated convergence, t → ψ n,p (α(t)) is differentiable and its differential in α at w is given by

T α F n,p (w) = n i=1 -2 1 {x∈ Ci } -→ a i x p-1 -→ a i x -→ a i x µ(dx), w i = grad α F n,p , w ,
where •, • denotes the L 2 metric on M n , which gives the desired result for the gradient.

The algorithm

Assume that we have access to independent and identically distributed (i.i.d.) observations X 1 , X 2 , . . . sampled from µ. We choose a sequence of positive steps (γ k ) k≥1 ⊂ (0, 1) satisfying the usual conditions

k≥1 γ k = +∞, k≥1 γ 2 k < +∞. (8)
We propose the following algorithm.

Algorithm 1 (Competitive Learning Riemannian Quantization). Input: i.i.d. observations X 1 , X 2 , . . . sampled from µ.

1. Initialization: set k = 0 and select α(0) = (a 1 (0), . . . , a n (0)) at random.

a i (k + 1) = exp a i (k) γ k+1 -------→ a i (k)X k+1 , a j (k + 1) = a j (k) ∀ j i.
3. Compute the Vornoi cells C 1 , . . . , C n associated to the optimal solution α = (a 1 , . . . , a n ).

Output: optimal approximation μ = n i=1 µ(C i )δ a i .

The steps are chosen to be in [0, 1] so that at each iteration, the center that is updated stays in the same Voronoi cell, guaranteeing that the centers stay pairwise distinct (if initially pairwise distinct). Now we show the convergence of Algorithm 1, using a theorem from Bonnabel [START_REF] Bonnabel | Stochastic gradient descent on riemannian manifolds[END_REF].

Proposition 3. Assume that the injectivity radius of M is uniformly bounded from below by some I > 0, and let (α(k)) k≥0 be computed using Algorithm 1 and samples from a compactly supported distribution µ. Then F n,2 (α(k)) converges a.s. and grad α(k) F n,2 → 0 as k → ∞ a.s.

Proof. The proof relies on [5, Theorem 1], and therefore we adopt the same notations. The cost function F n,2 and its gradient can be respectively expressed as expectations of the functions

Q(x, α) = min 1≤i≤n d(x, a i ) 2 , H(x, α) = grad α Q(x, α) = -2 1 Ci (x) --→ a i x 1≤i≤n .
This is once again due to Lemma 1, given in the Appendix. As stated in Remark 1, the n-centers α(k) = (a 1 (k), . . . , a n (k)) are each the barycenter of their Voronoi cell and therefore always stay in the same compact ball B(a, R) as the data

X 1 , X 2 , . . .. Since 1 Ci (x) - → yx ≤ 2R for all y, x ∈ K ⊂ B(a, R) and i = 1, . . . , n, the gradient H is bounded on K × K. That
is, all the assumptions of [5, Theorem 1] are verified and Algorithm 1 converges.

Link to clustering

Notice that if the measure µ one wants to approximate is the empirical measure of a dataset x 1 , . . . , x N , N n,

µ = 1 N N k=1 δ x k ,
then the weight of the atom a i in the optimal approximation is simply given by the proportion of data points that fall inside its Voronoi cell, that is,

μn = n i=1 |C i (Γ)| N δ a i .
The atoms a 1 , . . . , a n minimize the distortion function, i.e., in this case, the sum of the squared geodesic distances within each Voronoi cell

F n,2 (a 1 , . . . , a n ) = n i=1 x k ∈C i d 2 (x k , a i ). ( 9
)
Notice that this is exactly the cost function of K-means clustering, to which algorithm 1 (where the observations X k are simply chosen randomly among x 1 , . . . , x N ) therefore provides an efficient alternative. Just as for classical K-means, the clusters are given by the Voronoi cells.

Examples

Now let us show some toy examples on manifolds of constant sectional curvature: the circle, the 2-sphere and the hyperbolic plane. Figure 2 shows results of optimal quantization of size n = 6 of the uniform (top row) and the von Mises (bottom row) distributions on the circle. The considered von Mises distribution is centered in 0 with concentration parameter K = 5. From left to right, Figure 2 shows the initialization of Competitive Learning Riemannian Quantization (CLRQ), the optimal centers obtained after convergence, and the clustering (or Voronoi diagram) of the data taken as input of the algorithm. The centers are initialized uniformly on (-π, π] and [-π/2, π/2] respectively, S 1 being identified with (-π, π]. In order to reduce dependency on the initialization, each step k is repeated a certain number m of times. In other words, the same step size is used several times. In the uniform case, m = 10 is sufficient to obtain a visually satisfying result. For the von Mises distribution, we choose m = 50.

On the right-hand side of Figure 2, we plot the evolution of the Wasserstein distance between the initial distribution and its quantized version. (Recall that the quantization cost function (5) involves the L 2 -Wasserstein distance.) As shown in [START_REF] Rabin | Transportation distances on the circle[END_REF], the computation of the L p -Wasserstein distance between two measures on the circle can be reduced to the same operation in the unit interval [0, 1] by "cutting" the circle at a certain point s ∈ S 1 , i.e., by identifying it with a unit length fundamental domain for R/Z. However, when the two measures are not evenly distributed, the optimal cut is easier to find in the L 1 case, therefore we choose to merely compute the weaker L 1 -Wasserstein (or Kantorovich-Rubinstein) distance using the algorithm introduced in [START_REF] Cabrelli | The kantorovich metric for probability measures on the circle[END_REF]. We plot the distance between the measure µ and its approximation at each step k of Algorithm 1

ν(k) = n i=1 µ{C i (k)}δ a i (k) ,
where (a 1 (k), . . . , a n (k)) are the n-centers at step k, C 1 (k), . . . , C n (k) are the corresponding Voronoi cells, and δ x is the Dirac distribution at x ∈ M. Assuming that N n, we can approximate µ by the empirical measure of the observations x 1 , . . . ,

x N μ = 1 N N i=1 δ x k .
In order to compare two discrete measures with the same number of points, we then identify μ and ν(k) with the measures obtained on the reunion of their supports by completing with zero masses. For both the uniform and the von Mises examples, the Wasserstein distance decreases as expected.

Next, we show examples on the sphere. Figure 3 displays, in the same order as Figure 2, the initialization and results of the CLRQ algorithm applied to the von Mises distribution with concentration parameter K = 5, where each step contains m = 20 iterations at the same step size. Finally, to show an example in the negative curvature setting, we use the work of Said et al. [START_REF] Said | Riemannian gaussian distributions on the space of symmetric positive definite matrices[END_REF] regarding the definition and generation of Gaussian distributions on the space of SPD matrices to generate Gaussian samples on the hyperbolic half-plane. Recall that the hyperbolic half-plane is one of the models of 2-dimensional hyperbolic geometry, and can be identified with the complex upper-half plane H 2 = {z = x + iy, (x, y) ∈ R × R * + } where the length of an infinitesimal element dx + idy at point x + iy is measured by

ds 2 = dx 2 + dy 2 y 2 .
The special linear group SL 2 acts on H 2 from the left through the Moebius transformation:

SO 2 × H 2 → H 2 , defined by a b c d • z = az + b cz + d .
This action is transitive since for all (x, y) ∈ R × R * + ,

y 1/2 xy -1/2 0 y -1/2 • i = x + iy.
Noticing that SO 2 is the stabilizer of i, we can identify H 2 with SL 2 /SO 2 , which is also homeomorphic to the space of 2 × 2 SPD matrices of determinant 1 [START_REF] Jost | Riemannian geometry and geometric analysis[END_REF]. The space of SPD matrices of determinant 1 is therefore homeomorphic to H 2 , and the homeomorphism is given by Φ : P → L • i, where L is the upper-triangular matrix of the Cholesky decomposition of the SPD matrix P = L L. To generate a Gaussian sample in H 2 , we generate a Gaussian sample of SPD matrices using [START_REF] Said | Riemannian gaussian distributions on the space of symmetric positive definite matrices[END_REF]Proposition 6], renormalize and transport them to the hyperbolic half-plane using Φ. The first row of Figure 4 shows the initialization and results of the CLRQ algorithm applied to observations sampled from the Gaussian distribution centered in i and with standard deviation σ = 0.5. Each step contains m = 20 iterations.

Application to air traffic complexity analysis

Context

This work was motivated by complexity analysis in air traffic management (ATM). ATM deals with the means of organizing the aircraft trajectories in a given airspace so as to ensure both safety and efficiency. One of the most important part of ATM is the air traffic control (ATC) that acts on flying or taxiing aircraft in such a way that separation norms are satisfied at all time. Nowadays, most of the ATC is surveillance based, relying primarily on information coming from the Radars to give instructions to planes. Even in country-sized airspaces the amount of traffic to be controlled is far beyond the limits of a single operator and the area under the responsibility of an air traffic controller (ATCO) has to be kept tractable. As a consequence, the whole airspace must be divided into elementary cells, known as control sectors, that are alloted to a team of ATCOs. A major concern arising in the design of such sectors is to ensure that the ATCOs workload is equally balanced over the airspace. Highly complex areas, involving many flows crossings and altitude changes, like those encountered close to the largest airports, must be kept small, while low complexity ones, with either a small aircraft density or a simple airspace structure may be made large. Finding a complexity indicator that can be computed using only airspace and traffic information and that closely mimics the perceived workload is a difficult, still unsolved problem [START_REF] Cook | Applying complexity science to air traffic management[END_REF][START_REF] Prandini | Toward air traffic complexity assessment in new generation air traffic management systems[END_REF].

In France, the current operational indicator is the number of aircraft present in a control section at a given time. It is used by the European network manager to compute departure slots so that the value given is unlikely to be exceeded. This crude indicator is known to be too coarse as in some situations, ATCOs are willing to accept more aircraft than the threshold, while in others, they experience control overload with fewer aircraft. This is mainly due to differences in the traffic structure, that may be highly organized, thus easy to control, or on the contrary close to random and requiring a lot of attention. It is worth noticing that an incomplete knowledge of the complexity perceived by the ATCOs induces unnecessary delays, as some take-off slots would have been improved if the real acceptable number of aircraft in control sectors had been known.

In research, one of the most widely used indicators is the dynamic density [START_REF] Laudeman | Dynamic density: An air traffic management metric[END_REF]. It combines influential factors, like the number of maneuvering aircraft, the number of level changes and so on, to output a single positive real value representing the complexity level. Although quite pertinent from an operational point of view, the dynamic density is a tool requiring a lot of tuning, involving experiments with a panel of ATCOs. These experiments would be costly to implement in practice and cannot be adapted to different airspaces without having to re-tune from scratch. For the same reason, it is quite difficult to use dynamic density for assessing the performance of new concepts, since in such a case there is no reference situation or existing operational context that may be used to perform the tuning phase.

On the other hand, purely geometrical indicators have been introduced [START_REF] Delahaye | Air traffic complexity based on dynamical systems[END_REF][START_REF] Lee | Air traffic complexity : An input-output approach[END_REF], that are free of the dynamic density limitations. While perfectly suited to infer an intrinsic notion of complexity, they do not model all the aspects of the workload, as perceived by a human operator. The approach taken in the present work may be viewed as a mix between the two previous ones: it relies on an intrinsic description of traffic, but does not produce a complexity value. Instead, a summary of the situation is issued, that serves as an input to a classification or machine learning algorithm. Even if this last phase looks very similar to a kind of dynamic density evaluation, it is far less complex:

• The traffic summary itself requires no tuning, while the influential factors taken into account in the dynamic density have weights that must be adjusted during the learning phase.

• Complexity classes are computed instead of complexity levels: a clustering algorithm will first segment the traffic dataset into homogeneous subsets, then a workload value will be associated by experts to each of them. This process is lightweight, since only the representative in each class has to be expertized.

• Adaptation to different airspaces is an easy task for the same reason: experts will evaluate only the representative situation in each class.

The first step is to model the spatial distribution of the aircraft velocities as a Gaussian law. Then, the covariance function is used as an indicator of traffic complexity. This assumption makes sense from an operational point of view as it represents the local organization, that is the most influential factor on ATCCOs workload.

Estimating the covariance matrices

Although we will in practice consider time-windows, we start by considering a given airspace at a fixed time t containing N aircraft either flying or taxiing. We respectively denote by z i and v i the position and speed of the aircraft i, 1 ≤ i ≤ N, at time t. Since the altitude of an aircraft plays a special role and does not appear on controllers displays, we choose to disregard it and adopt a planar representation through a stereographic projection. An underlying Gaussian field model is assumed for the relation between the velocity and the position, whose variance function will be interpreted as a pointwise measure of the traffic complexity.

A non parametric approach of type Nadaraya-Watson [START_REF] Nadaraya | On estimating regression[END_REF][START_REF] Watson | Smooth regression analysis[END_REF] was taken to estimated the mean and variance functions at point z:

m(z) = N i=1 V i K h (z -Z i ) N j=1 K h (z -Z j ) , Σ(z) = N i=1 {V i -m(z)}{V i -m(z)} K h (z -Z i ) N j=1 K h (z -Z j ) . ( 10
)
The weights are given by a kernel function K, i.e., a positive, symmetric function of unit area, scaled by a factor h > 0:

K h (x) = h -1 K(x/h).
Since most kernels have compact support, the estimations are based in practice on a number of observations that is very inferior to the size N of the sample. The estimator Σ has been studied in [START_REF] Yin | Nonparametric covariance model[END_REF] where it is shown to be asymptotically normal. Evaluating it at positions z 1 , . . . , z N yields a series of symmetric, positive definite matrices Σ(z 1 ), . . . , Σ(z N ) with empirical distribution

μ = 1 N N i=1 δ Σ(z i ) ,
where δ Σ denotes the Dirac mass at Σ. In order to obtain a summary of the traffic complexity, we propose to quantize μ using the CLRQ algorithm on the space of SPD matrices.

The geometry of SPD matrices

For the sake of completeness, let us briefly recall the most commonly used Riemannian structure [START_REF] Pennec | A riemannian framework for tensor computing[END_REF] on the space P d of d × d symmetric, positive definite matrices. Note that in this application, we are simply interested in the case d = 2. The Euclidean dot product on the space M d of square matrices of size d is given by the Frobenius inner product

Σ 1 • Σ 2 = tr(Σ 1 Σ 2 )
, where tr denotes the trace. As an open subset of the vector space M d , P d is a manifold where the tangent vectors are symmetric matrices. It can be equipped with a Riemannian metric invariant with respect to the action of the general linear group GL d × P d → P d , (A, Σ) → A ΣA. At the identity, this metric is given by the usual Euclidean scalar product W 1 , W 2 Id = W 1 • W 2 = tr(W 1 W 2 ), and at Σ, we ask that the value of the scalar product does not change when the tangent vectors are transported back to the identity via the action of A = Σ -1/2 , i.e.,

W 1 , W 2 Σ = Σ -1/2 W 1 Σ -1/2 , Σ -1/2 W 2 Σ -1/2 Id = tr(Σ -1/2 W 1 Σ -1 W 2 Σ -1/2 ). (11)
The associated geodesic distance is given by

d(Σ 1 , Σ 2 ) = d i=1 log 2 λ i (Σ -1/2 1 Σ 2 Σ -1/2 1
) ,

where we use the notation λ i (Σ), i = 1, . . . , d, to denote the eigenvalues of Σ. Recall that in order to update the centers of the discrete approximation in the CLRQ algorithm, we need the exponential map, i.e., a mapping that associates to each point Σ and tangent vector W at Σ the end point of the geodesic starting from Σ at speed W. In the case of metric [START_REF] Cheeger | Comparison theorems in Riemannian geometry[END_REF], it is given by exp

Σ (W) = Σ 1/2 exp Σ -1/2 WΣ -1/2 Σ 1/2 ,
where the exp on the right-hand side denotes the matrix exponential. Finally, we also need the inverse mapping, i.e., the logarithm map log

Σ 1 (Σ 2 ) = ---→ Σ 1 Σ 2 = Σ 1/2 1 log Σ -1/2 1 Σ 2 Σ -1/2 1 Σ 1/2 1 , 12
where the log on the right-hand side denotes the matrix logarithm. Note that the matrix logarithm is well defined for any symmetric matrix Σ due to the factoring out, in the logarithm series, of the rotation matrices of the spectral decomposition Σ = UDU , yielding log(UDU ) = U log(D)U . Endowed with this Riemannian structure, the space of SPD matrices is a geodesically complete Riemannian manifold with infinite injectivity radius at every point (i.e., no cut locus), and so all the exposed theoretical results are valid in this case.

Real data analysis 5.4.1. Segmenting and constructing summaries

We now have all the tools to construct summaries of the traffic complexity in a given airspace during a certain time period. As input, we consider an image such as the ones displayed in the first row of Figure 5, showing the traffic over Paris, Toulouse and Lyon during a one-hour period of time. The color is related to the norm of the velocity, increasing from yellow to red. To simplify, we center and reduce the velocities v i . The samples (z i , v i ) are seen as observations arriving in a random order, and the covariance matrix at z i is estimated according to (10) using a truncated Gaussian kernel K(x) = 1/ √ 2π e -x 2 /2 1 |x|<r . The truncation of size r avoids useless computations.

Algorithm 2 (CLRQ for air traffic complexity analysis). Input: observations of position z 1 , . . . , z N and speed v 1 , . . . , v N .

1. Choose z (0) 1 , . . . , z (0) n randomly among the z i 's and compute Σ(z (0) 1 ), . . . , Σ(z (0) n ) using (10). 2. Set k = 0 and A 1 (0) = Σ(z (0) 1 ), . . . , A n (0) = Σ(z (0) n ).

Repeat until convergence:

(a) Choose z (k+1) at random among the z i 's and compute Σ k+1 = Σ(z (k+1) ) using [START_REF] Chazal | Persistence-based clustering in riemannian manifolds[END_REF].

(b) Find j min = argmin j d Σ k+1 , A j (k) . (c) Set A(k + 1) = (A 1 (k + 1), . . . , A n (k + 1))
where

A j min (k + 1) = exp A j min (k) γ k+1 ----------→ A j min (k)Σ k+1 , A j (k + 1) = A j (k) ∀ j j min .
(d) Compute the Voronoi cells C 1 , . . . , C n associated to the optimal solution A = (A 1 , . . . , A n ).

Output: density approximation μ = n i=1

|C i ]
N δ A i and clustering given by the Voronoi diagram {C 1 , . . . , C n }. In practice, we usually look for n = 3 centers, i.e., the best approximation by three points. Indeed, we have found that generically, the centers can be ordered for the Loewner order when n = 3 but not always for n > 3. (Recall that the Loewner order is a partial order on the space of SPD matrices, according to which A ≥ B if A -B is positive semi-definite.) This can be explained by the fact that the Riemannian metric [START_REF] Cheeger | Comparison theorems in Riemannian geometry[END_REF] sorts by rank, and therefore the covariance matrices are segmented into those with both eigenvalues close to zero (since the velocities are centered with respect to the neighborhood average, parallel velocities yield an estimated covariance matrix close to zero), those with only one eigenvalue close to zero, and the others. In the second row of Figure 5, these clusters are respectively shown in green, blue and red. As could be expected, the first cluster corresponds to zones with either an isolated trajectory or parallel trajectories, the second to simple crossings or variations of speed in a given trajectory, and the third to zones with high density and crossings involving many trajectories. Naturally, the choice of the size r of the kernel's support has a great influence on the clustering, and it should be adjusted according to the minimum distance authorized between two aircraft in a zone considered as non conflictual.

Further results are shown in Figure 6, where we consider the whole French airspace over different one-hour periods of time. The clusterings obtained using CLRQ is shown in the middle column. To illustrate the importance of the Riemannian setting with respect to the Euclidean one, we show results of Competitive Learning Vector Quantization on the same datasets, i.e., the same algorithm where the centers are updated using straight lines (linear interpolations between the matrix coefficients) and the distances are computed using the Frobenius norm. These results are shown in the right column of Figure 6, and are less convincing. Many crossings and convergence points of different trajectories, such as the area around Paris, are put in the same class as isolated trajectories by the algorithm. This is contrary to what we want, since the former present a high level of complexity for controllers in contrast to the latter. In Figure 7, we show that the initialization has little influence on the segmentation of the airspace, which is satisfactory.

Comparing summaries

Let us now propose a way to compare our different summaries. A natural way to do so is through discrete optimal transport, which allows one to compute the distance between two discrete measures

µ = m i=1 µ i δ A i and ν = n j=1 ν j δ B j .
In our case, the A i 's and B j 's are SPD matrices. Optimal transport seeks to transport the mass from µ to ν in a way that minimizes a certain cost. Formally, a transport plan between µ and ν is a matrix π = (π i j ) i, j with non-negative entries that verifies the two following properties for all i = 1, . . . , m and j = 1, . . . , n, n j=1 π i j = µ i and m i=1 π i j = ν j .

The set of transport plans between µ and ν is denoted by Γ(µ, ν). Intuitively, the value π i j represents the part of the mass µ i transported from A i to B j to reconstruct ν j . Here, we measure the cost of transporting an infinitesimal unit of mass from A i to B j by the corresponding geodesic distance d(A i , B j ). The optimal transport plan is chosen to minimize the global cost, i.e., to be a solution of the following minimization problem

W 2 (µ, ν) = min          m i=1 n j=1 π i j d(A i , B j ) 2 ; π ∈ Γ(µ, ν)          . ( 12
)
The value of this minimal cost is the discrete L 2 -Wasserstein distance between µ and ν. In general, the linear programming problem ( 12) is difficult and many different algorithms have been developed to solve it in various special cases [START_REF] Mérigot | Discrete optimal transport: complexity, geometry and applications[END_REF]. However, in our setting, the discrete measures involved are supported by a very small number of points (3 in the examples shown) and problem [START_REF] Cook | Applying complexity science to air traffic management[END_REF] presents no difficulties. The distances between the summaries of the different traffic situations shown in the middle column of Figure 6 are given in Table 5.4.2. As expected, the first two situations are deemed similar but very different from the third one which has a much more complex traffic. The last situation is intermediary. In comparison, the different summaries of Figure 7 obtained for the same traffic situation but different initializations are at small distances from one another, as shown in Table 5.4.2.

Evaluation of the quantization results

The French territory is divided in five en-route control centers (North, East, West, South-East and South-West), which are responsible for controlling aircraft in the corresponding area of the French airspace. Each of them is partitioned into airspace blocks, known as elementary control sectors, and designed in such a way that the heaviest expected traffic is still manageable. For that purpose, the geometry of the elementary sectors is carefully chosen, so that the number of aircraft is small enough, but at the same time crossing points are far apart the boundaries. Since elementary control sectors are adapted to worst-case situations, they are no longer optimal for light to medium loaded traffic. Grouping of sectors is then triggered, to obtain larger airspace blocks thus reducing the number of ATCOs needed. On the opposite, when the traffic starts to densify, degrouping occurs to cope with the increased load. This mechanism has a major effect on controllers workload, and is easily observed when conducting experiments with control simulators. It is also a good a posteriori indicator of complexity.

We evaluate the performance of the CLRQ algorithm on a dataset corresponding to the air traffic over the Eastern control center during a time period of a few hours. First, we compare the clustering given by quantization to other clustering methods. Figure 8 shows, from left to right, the clusterings obtained for a fifteen-minute time slot by CLRQ, K-means, and Kernelized Gaussian Random Projection (KGRP) [START_REF] Zhao | Efficient clustering on riemannian manifolds: A kernelised random projection approach[END_REF]. KGRP projects the data points in a RKHS (Reproducing Kernel Hilbert Space) generated by a randomly selected subset of the data points, before performing K-means. As suggested in [START_REF] Zhao | Efficient clustering on riemannian manifolds: A kernelised random projection approach[END_REF], we use a Gaussian kernel, we set the dimension of the RKHS to 300, and we use a subset of 100 points to generate it. Without optimizing any of the corresponding codes, CLRQ improves computation time by a factor of at least 10 with respect to K-means and KGRP. Both of these could be made faster with a better optimized computation of the Voronoi diagrams computed at each step, but CLRQ could also be optimized through parallel computations. Comparison of optimized versions of these codes is postponed to future work. The clustering results of CLRQ and K-means are very similar, and both algorithms attain approximately the same value of intra-cluster variance (9) (6.489 and 6.481 respectively). In contrast, the KGRP method yields a much larger value (approximately 12.3), since the K-means is performed in the projection space instead of the space of SPD matrices.

Finally, we confront our results to the perception of an ATCO. For this dataset, we have access to an ATCO's perception of the complexity of the traffic, evaluated every five minutes, in real time, in the form of an index between 0 and 5. In Figure 9, the evolution of this perceived complexity index is shown by the blue curve. Above, we display the results of optimal quantization performed on fifteen time slots of fifteen minutes each. Each bar diagram describes a cluster, or equivalently, an atom of the discrete distribution that best approximates the data, with the same color code as previously defined. In each bar diagram, the heights of the bars indicate the proportions of the data in the cluster, or equivalently, the weights associated to the atom, while their widths correspond to the distance from the atom (or cluster center) to the smallest atom among all those displayed. Finally, the bottom curve displays the evolution of the number of aircraft present in the control zone for 5-minute time slots.

Since the ATCO's perception of the traffic complexity can be influenced by the recent past as well as by the anticipation of upcoming situations, it is not easy to infer the true real-time complexity of the traffic from it. In other words, a non-zero time interval between reality and human reaction seems inevitable. However, it seems that the two drops of complexity perceived by the ATCO around 7:21 and just after 8:51 are visible from the quantization summaries, in particular from the weights of the first cluster. The first of these two time periods corresponds to a decision to regroup, probably motivated by the simplification of the traffic. Peaks of perceived complexity around 8:00, just before 8:51 and around 9:51 also seem to coincide with the quantization results. The first of these times corresponds to a decision to degroup, probably in reaction to the high complexity.

Also, it is important to notice that the complexity perceived by the ATCO is not directly indexed on the number of aircraft in the considered airspace, and neither is the quantization summary. For example, using the Wasserstein distance we can notice that the summary of the 10 th time slot (from 8:36 to 8:51) is slightly more similar to (at distance 1.49 of) the summary of the last time slot (from 9:51 to 10:06), despite an increase of 42% in the number of aircraft, than to the summary of the following 11 th time slot (8:51 to 9:06, at distance 1.57), which represents an increase of only 8%. This is already a significant improvement on the operational complexity index based on the number of aircraft.

Conclusion and future work

We have proposed a gradient descent type algorithm to find the best finite discrete approximation of a probability measure on a Riemannian manifold. This algorithm is adapted to large sets of data as it is online, and yields a clustering on top of a finite summary of the non-Euclidean data of interest. It is convergent when the manifold is complete and when its injectivity radius is uniformly bounded from below by a positive constant. We have used it to compute summaries of air traffic images in the form of a small number of covariance matrices representing different levels of local complexity, with associated weights corresponding to the occurrence of these levels of complexity in the images. These summaries can then be compared through discrete optimal transport, yielding a simple way to compare and index air traffic situations.

A still open question is to find an optimal number of centers for the quantization procedure. A possible approach is to use the associated clustering and to apply standard quality indicators like the silhouette [START_REF] Peter | Silhouettes: A graphical aid to the interpretation and validation of cluster analysis[END_REF]. In such a case, one starts with a minimal guess on the number of centers and increases it until the clustering is good enough.

In future work, we will also consider best finite constrained quantization, i.e., restriction to finite approximations with equal weights. Within this frame, one wants to find an optimal set of centers (a 1 , . . . a n ) such that the distribution:

1 n n i=1
δ a i best approximates the true sample density. In the air traffic complexity application, this may be used for finding an optimal segmentation of the airspace in areas of equal complexity, thus making the process of control sector design fully automatic. It may also be used in a dynamic way, since the quantization process presented above is an on-line procedure and will self-adapt to changing situations. Changes in the shape of control sectors must however be limited to allow ATCOs to seamless switch to new airspace partitions. This last point is currently under investigation in the context of manual or semi-automatic dynamic airspace design, so that our procedure may benefit from the results of this study.

Finally, only the covariance information was used in the application, but the mean function that yields the expected direction, and the local density bear some important information. Working on the product manifold R + ×R d ×P d (with d = 2 or d = 3 for the air traffic framework) allows taking them into account. The gradient algorithm is readily adapted to this case, since geodesics on products are easily seen to be curves whose projections on factors are geodesics. A fully fledged complexity indicator may then be issued. Proof. Let a ∈ M, u ∈ T a M, and construct a curve (-, ) t → a(t) passing through a(0) = a at speed a (0) = u, so that the differential of f at a in u is written, in terms of the norm associated to the Riemannian metric, starting at γ(t, 0) = x and ending at γ(t, 1) = a(t). The vector field J(s) := ∂γ/∂t(0, s) then defines a Jacobi field along γ(s) := γ(0, s), measuring the way the geodesics spread out from x, and its initial velocity is given by ( Now for any i j, the subset ∂H(a i , a j ) = {x ∈ K, d(x, a i ) = d(x, a j )} defined as the kernel of the submersion x → d(a i , x)-d(a j , x) is a submanifold of M of codimension 1, yielding µ{∂H(a i , a j )} = 0 and therefore µ{∂C i (α)} = 0.

Figure 1 :

 1 Figure 1: Traffic over the French airspace during one-hour periods of time. The speed of the aircraft increases from yellow to red.

2 .

 2 For each new observation X k+1 , repeat until convergence: (a) find i = argmin j d(X k+1 , a j (k)), (b) set α(k + 1) = (a 1 (k + 1), . . . , a n (k + 1)) where

Figure 2 :

 2 Figure 2: Quantization of the uniform (top) and von Mises (bottom) distributions. From left to right: initial positions of the n-centers, final positions of the n-centers, associated Voronoi diagram, and evolution of the L 1 -Wasserstein distance between the initial distribution and its quantized version.

Figure 3 :

 3 Figure 3: Quantization of the von Mises distribution on the 2-sphere. From left to right: initialization of the centers, optimal centers, and the associated Voronoi diagram.

Figure 4 :

 4 Figure 4: Quantization of the Gaussian distribution on the hyperbolic half-plane. From left to right: initialization of the centers, optimal centers, and the associated Voronoi diagram.

Figure 5 :

 5 Figure 5: Clustering of the airspaces over Paris (left), Toulouse (middle) and Lyon (right).

Figure 7 :

 7 Figure 7: Clustering of the French airspace with 3 different initializations.

Figure 8 :

 8 Figure 8: Clustering of the French Eastern control center from 6:36 to 6:51 AM, obtained by 3 different clustering methods (from left to right): CLRQ, K-means, and Kernalized Gaussian Random Projection.

Figure 9 :

 9 Figure 9: Results of quantization of the French Eastern control center for 15 time slots of 15 minutes each.

2 ,

 2 where v(t) := ---→ xa(t) = log x a(t). Consider the sheave of geodesics γ(t, s) := exp x {sv(t)}, (t, s) ∈ (-, ) × [0, 1],

Lemma 2 .

 2 x ) sv(0) {sv (0)} = D ds s d(exp x ) sv(0) {v (0)} .Identifying the tangent spaces T x M T v(0) (T x M), we getDJ ds (0) = d(exp x ) 0 {v (0)} + 0 × D ds s=0 d(exp x ) sv(0) {v (0)} = v (0),and so,T a f (u) = 2 v (0), v(0) = 2 DJ ds (0), γ(0).Since J verifies the classical Jacobi equation in terms of the curvature tensor R of M, D 2 J/ds 2 = -R(J, γ )γ , we get D 2 J/ds 2 , γ = 0, leading to the scalar product DJ/ds, γ being constant (since Dγ /ds = 0). This gives, sinceJ(0) = 0, J(s), γ (s) = DJ ds (0), γ (0) s + J(0), γ (0) = DJ ds (0), γ (0) s,yielding DJ/ds(0), γ (0) = J(1), γ (1) . To conclude, notice that J(1) = ∂γ/∂t(0, 1) = u and, setting γ(s) = γ(1s), that γ (1) = -γ (0) = --→ ax and so finally T a f (u) = -2 u, -→ ax , which completes the proof. Let µ be a probability measure with density with respect to the Riemannian volume form on M, α = (a 1 , . . . , a n ) ∈ M n and C i (α) denote the Voronoi cell associated to a i for all i = 1, . . . , n. No mass is assigned by µ to the boundaries of the Voronoi diagram µ{∂C i (α)} = 0.Proof. For any n-tuple α = (a 1 , . . . , a n ), the i th Voronoi cell can be writtenC i (α) = ∩ j i H(a i , a j ), where H(a, b) = {x ∈ K, d(x, a) ≤ d(x, b)},which gives ∂C i (α) = C i (α) ∩ Ci (α) c = C i (α) ∩ {∩ j i H(a i , a j )} c = ∪ j i ∂H(a i , a j ) ∩ C i (α).

Table 1 :

 1 Distances between the summaries of Figure6.

	0.00 1.92 6.74 4.55
	1.92 0.00 8.31 6.07
	6.74 8.31 0.00 1.22
	4.55 6.07 1.22 0.00
	0.000 0.033 0.031
	0.033 0.000 0.016
	0.031 0.016 0.000

Table 2 :

 2 Distances between the summaries of Figure7.

Figure 6: Clustering of the French airspace using Riemannian quantization (middle) versus vector quantization (right).

Acknowledgements

We would like to thank Salem Said and Lionel Bombrun for sharing with us their code to generate Gaussian samples in the hyperbolic half-plane.

Pattern Recognition, 51:333-345, 2016.


G D Mcbain

S G Mallinson

Impulsively Generated Incompressible Two-Phase Flow and the Asai Thermal Ink-jet Model

A widespread model of thermal ink-jet actuation prescribes the pressure of the expanding bubble as a rapidly decreasing function of time, this being taken as a finite version of a Dirac impulse with only the integral with respect to time mattering. Here it is shown that the impulse is actually easier to treat numerically than the bounded approximation and an implementation in OpenFOAM is described; however, solutions so obtained reveal that the flow does depend on the timescale over which the pressure is applied, i.e. on first and higher moments of the pressure function. An asymptotic expansion of a lumped linear model is derived which predicts the dependence of the induced flow on these moments for an arbitrarily shaped pressure pulse.

Introduction

A thermal ink-jet is a device for ejecting microscopic droplets, e.g. for printing [START_REF] Asai | Three-dimensional calculation of bubble growth and drop ejection in a bubble jet printer[END_REF][START_REF] Zhao | Pressure and power generation during explosive vaporization on a thin-film microheater[END_REF]. As sketched in figure 1, it is a vessel filled with ink, open at one end through a nozzle to the air and at the other to a reservoir. It contains a heater which suddenly boils ink near it so that the expanding bubble pushes ink forwards out through the nozzle. Following the nucleation of the bubble, the flow consists of three distinct phases: the vapour bubble, the liquid ink, and the ambient air. The ink and air are essentially incompressible but the overall dynamics is complicated by strong interaction between the thermodynamics of the bubble and its heat, mass, and momentum transfer with the ink.

To alleviate the complexity of the bubble, a much simpler decoupled model was introduced [START_REF] Asai | Three-dimensional calculation of bubble growth and drop ejection in a bubble jet printer[END_REF] in which it is replaced with a pressure which is a prescribed function of time; we also mention a number of precursors [START_REF] Beasley | Model for fluid ejection and refill in an impulse drive jet[END_REF][START_REF] Kyser | Design of an impulse ink jet[END_REF][START_REF] Fromm | A numerical study of drop-on-demand ink jets[END_REF][START_REF] Deshpande | Fluid mechanics of bubble growth and collapse in a thermal ink jet printhead[END_REF]. Although more recent computational models have been developed which include a coupled bubble [START_REF] Cornell | Nucleation quality and bubble momentum and their effect on droplet velocity and stability[END_REF][START_REF] Zhao | Pressure and power generation during explosive vaporization on a thin-film microheater[END_REF][START_REF] Sen | Droplet ejection performance of a monolithic thermal inkjet print head[END_REF]24,[START_REF] Zhou | Simulation model and droplet ejection performance of a thermal-bubble microejector[END_REF][START_REF] Tan | An adaptive mesh refinement based flow simulation for free-surfaces in thermal inkjet technology[END_REF], the uncoupled model continues to usefully simulate aspects of ink-jet actuation [START_REF] Deshpande | Significance of inertia and resistance in fluidics of thermal ink-jet transducers[END_REF][START_REF] Chen | Bubble growth and ink ejection process of a thermal ink jet printhead[END_REF][START_REF] Wu | Development of a three-dimensional simulation system for micro-inkjet and its experimental verification[END_REF][START_REF] Lindemann | Droplet Generation From the Nanoliter to the Femtoliter Range[END_REF][START_REF] Glatzel | Computational fluid dynamics (CFD) software tools for microfluidic applications-a case study[END_REF][START_REF] He | The roles of wettability and surface tension in droplet formation during inkjet printing[END_REF][START_REF] Jiang | One dimensional model for droplet ejection process in inkjet devices[END_REF].

Asai has written much on ink-jet technology [START_REF] Asai | One-dimensional model of bubble growth and liquid flow in bubble jet printers[END_REF][START_REF] Asai | Bubble generation mechanism in the bubble jet recording process[END_REF][START_REF] Asai | Application of the nucleation theory to the design of bubble jet printers[END_REF][START_REF] Asai | Bubble dynamics in boiling under high heat flux pulse heating[END_REF]; here 'the Asai model' is merely used as a concise name for the 'Method of Calculation' [START_REF] Asai | Three-dimensional calculation of bubble growth and drop ejection in a bubble jet printer[END_REF], and more specifically for modelling the bubble by a prescribed history of pressure. Originally it was suggested that 'the bubble pressure can be modeled by an impulsive function' [START_REF] Asai | Three-dimensional calculation of bubble growth and drop ejection in a bubble jet printer[END_REF], and this had been done in an earlier lumped model [START_REF] Asai | Bubble dynamics in boiling under high heat flux pulse heating[END_REF], but then recourse to generalized functions was avoided by arguing that 'numerically, such impulsive pressure change can be described by a rapidly decreasing function of time' [START_REF] Asai | Three-dimensional calculation of bubble growth and drop ejection in a bubble jet printer[END_REF] with the same pressure impulse (integral of pressure over time). While this is asymptotically true, the first purpose of the present paper is to show that actually the direct use of an impulsive pressure is by no means difficult.

The second purpose is to see whether the actuating pressure can indeed be characterized purely in terms of its impulse [START_REF] Deshpande | Fluid mechanics of bubble growth and collapse in a thermal ink jet printhead[END_REF] or its timing is significant [START_REF] Kyser | Design of an impulse ink jet[END_REF]; i.e. how slow the pressure has to be before its higher moments with respect to time become relevant.

The Asai Ink-Jet Model

The method presented here differs from the original 'Method of Calculation' [START_REF] Asai | Three-dimensional calculation of bubble growth and drop ejection in a bubble jet printer[END_REF] in which the equations of motion were only solved in the ink with a constant pressure imposed beyond the meniscus outside the nozzle and a varying pressure in the bubble. These days, it is more usual to fill out the domain with a second immiscible incompressible air-like fluid outside the nozzle and inside the bubble using the 'one-fluid formulation' [START_REF] Bierbrauer | A solenoidal initial condition for the numerical solution of the Navier-Stokes equations for two-phase incompressible flow[END_REF][START_REF] Tryggvason | Direct Numerical Simulation of Gas-Liquid Multiphase Flows[END_REF] with density ρ = ρ(α) and viscosity µ = µ(α) functions of the volume-of-fluid α interpolating between the values for air and water at α = 0 and 1, respectively. The surface tension, previously part of the boundary condition on the interface, is represented by a continuous volumetric force [START_REF] Brackbill | A continuum method for modeling surface tension[END_REF] . The essential feature of the model remains the prescription of the pressure in the bubble; here this is shifted from the interface to the surface of the actuator inside the bubble.

Governing Equations

The system for the evolution of the velocity u in time t is thus

ρ Du Dt = ρF -∇p + ∇ • µ ∇u + [∇u] T (1)
∇ • u = 0 (2) Dα Dt = 0, (3)
where F includes capillarity.

In the model, the pressure p is prescribed far off in the air as atmospheric, at the upstream ink reservoir as some other constant, and at the surface of the heater as a rapidly decreasing function of time such as [START_REF] Asai | Bubble dynamics in boiling under high heat flux pulse heating[END_REF][START_REF] Asai | Three-dimensional calculation of bubble growth and drop ejection in a bubble jet printer[END_REF]

p(t) = p 0 e -(t/t 0 ) λ (4)
plus a slower varying more bounded part which is of less interest here. Where the velocity crosses the boundary into the domain, a value for α is prescribed; α = 0 for the heater. On the other walls, u = 0 and a contact angle is set. This is a standard set of equations and canned solvers are available; here interFoam from the free OpenFOAM suite was used in figure 2 to repeat Asai's [START_REF] Asai | Three-dimensional calculation of bubble growth and drop ejection in a bubble jet printer[END_REF] case λ = 1 2 , t 0 = 0.17 µs.

Impulsively Generated Flow

The form of the pressure equation ( 4) is a bounded approximation to a Dirac δ function, approaching Πδ(t) as t 0 → 0 with λ fixed and p 0 varying inversely with t 0 so that the pressure impulse Π ≡ m 0 is constant, where m k is the k-th moment:

m k ≡ ∞ 0 t k p(t) dt. (5)
In the Navier-Stokes equation ( 1), if the energy and therefore velocity is to remain finite, a Dirac pressure spike can only be balanced by the temporal derivative of the velocity with the latter having a Heaviside step [START_REF] Lamb | Hydrodynamics[END_REF]; thus asymptotically,

ρ ∂ {uH(t)} ∂t = -∇ {Πδ(t)} (6)
or u = -ρ -1 ∇Π.

An equation for the pressure impulse field Π is obtained by substituting equation ( 7) for the impulsive velocity into the continuity equation ( 2):

-∇ • ρ -1 ∇Π = 0. (8)

Impulsively Generated Flow in Inhomogeneous Fluid

For a single homogeneous fluid, ρ can be taken outside the divergence in equation ( 8) and the pressure impulse satisfies Laplace's equation [START_REF] Lamb | Hydrodynamics[END_REF]; for inhomogeneous fluid, it must be retained as a spatially varying coefficient.

Equation ( 8) was previously derived for the impulsive pressure in a stratified fluid [START_REF] Saffman | Kinetic energy generated by the incompressible Richtmyer-Meshkov instability in a continuously stratified fluid[END_REF]. It also arose in another discussion of computational multiphase flow in trying to filter out the nonsolenoidal component of an initial condition [START_REF] Bierbrauer | A solenoidal initial condition for the numerical solution of the Navier-Stokes equations for two-phase incompressible flow[END_REF].

Laplace's equation would suffice for a post-impulse nonzero initial condition for the velocity in Asai's [START_REF] Asai | Three-dimensional calculation of bubble growth and drop ejection in a bubble jet printer[END_REF] liquid-only simulation, but if the air is to be included as well, as it must be with one-fluid solvers like interFoam, the density must be retained as a spatially varying coefficient. Although the physical effect of the air on the ink is indeed essentially limited to its pressure, difficulties will arise if the initial condition on velocity does not satisfy the continuity equation (2) [START_REF] Bierbrauer | A solenoidal initial condition for the numerical solution of the Navier-Stokes equations for two-phase incompressible flow[END_REF].

Solving for the Impulsively Generated Flow in OpenFOAM

While none of the standard OpenFOAM solvers can solve equation [START_REF] Brackbill | A continuum method for modeling surface tension[END_REF], little had to be done to generalize potentialFoam; indeed, only one line in the main C++ source file was modified, to include the varying coefficient of equation ( 8):

fvm::laplacian(1.0/rho, Phi)

although other subsidiary files had to be included from interFoam in order to set up the spatially varying rho.

The Two Approaches to Initial Impulses

The concept of having two equivalent sets of initial conditions, before and after an impulse, is common in electrical circuit analysis [19, pp. 470-480]. Mechanically, consider golf: one might contrast, on the one hand, simulating a club hitting a ball, the contact problem and elastic deformation, with, on the other, the comparatively simple ballistics of the rigid ball after it has received its momentum. Similarly in the bubble jet, it is easier to compute the impulsively generated velocity field by solving the potential-like equation ( 8) than to resolve in time the rapidly varying pressure of equation ( 4) and the correspondingly rapid acceleration of ink and air.

The Complex Lamellar Impulsively Generated Flow

Whereas flow impulsively generated in a uniform fluid is potential, equation ( 7), being the product of one scalar and the gradient of another is complex lamellar [START_REF] Thomson | A mathematical theory of magnetism-continuation of part I[END_REF]. Thus it is not necessarily irrotational. The vorticity is ∇ × u = -∇(1/ρ) × ∇Π; the impulse might set up a vortex sheet at an interface across which the density jumps. A couple of one-dimensional examples might be instructive. First consider two fluids in series: the impulsive pressure drops in proportion to the density and they set off with the same velocity, the whole being irrotational since the impulsive pressure gradient is normal to the interface. Second, two in parallel: the impulsive pressure drops are the same but the imparted velocities are inversely proportional to the density, with a vortex sheet along the interface.

Application to Modeling Ink-Jets

With a solver for equation [START_REF] Brackbill | A continuum method for modeling surface tension[END_REF] in hand, simulating Asai's prototype bubble jet becomes simpler: the preliminary solution provides the initial condition for velocity and the part of the pressure imposed on the heater corresponding to the impulse is omitted. Qualitatively the overall results are similar, but a closer quantitative examination in figure 3 using the flow-rate through the heater (as done by Deshpande [12, figures 3, 5, 7]) reveals that the flow-rate for Asai's pressure equation ( 4) with t 0 = 0.17 µs is not the same as that generated impulsively; however, the theory is asymptotically validated, being approached as the pressure is delivered more rapidly (t 0 → 0).

The pure impulse (t 0 = 0) provides a nonzero velocity at t = 0 and therefore a jump in flow-rate. The ratio of Π to this jump is as the inertance L. If the flow-rate in figure 3 is fitted to the log-linear form

q(t) ∼ Π L e -Rt/L , (9)
the quantity R has the dimensions of hydraulic resistance and L/R ≈ 6 µs is a hydraulic timescale.

This might seem long compared to the 0.17 µs proposed for t 0 [START_REF] Asai | Three-dimensional calculation of bubble growth and drop ejection in a bubble jet printer[END_REF] but unless λ = 1 in equation 4, t 0 is a misleading indicator of the timescale of the pressure input. A more physical value might be based on the ratio of the first two moments from equation ( 5); for equation ( 4) with = 1 2 , this is m 1 /m 0 = 6t 0 ≈ 1 µs. [START_REF] Deshpande | Fluid mechanics of bubble growth and collapse in a thermal ink jet printhead[END_REF] in the linear lumped model with pressure from (4) and ratios of timescales m 1 R/m 0 L chosen to roughly match figure 3 Impulsive and Phased Pressure A Lumped Model Equation ( 9) solves the lumped hydraulic equation [START_REF] Deshpande | Significance of inertia and resistance in fluidics of thermal ink-jet transducers[END_REF] L q + Rq = p(t) [START_REF] Cornell | Nucleation quality and bubble momentum and their effect on droplet velocity and stability[END_REF] with impulsive p(t) = Πδ(t) and q(0) = 0. For arbitrary p(t):

q(t) = L -1 t 0 e R(t ′ -t)/L p(t ′ ) dt ′ . ( 11
)
This is plotted in figure 4 with pressure equation ( 4) for values of the ratio of timescales (L/R)/(m 1 /m 0 ) roughly corresponding to those of figure 3.

As the lumped equation [START_REF] Deshpande | Fluid mechanics of bubble growth and collapse in a thermal ink jet printhead[END_REF] qualitatively reproduces the full multiphase three-dimensional dynamics of figure 3, it is used to investigate the failure of the hypothesis that only the zeroth moment of p(t) matters at times long after the pressure input has subsided-as seen there, the flow-rate after the passing of the pressure is a quarter higher than that due to impulse alone.

Asymptotic Expansion

Rewriting equation [START_REF] Deshpande | Fluid mechanics of bubble growth and collapse in a thermal ink jet printhead[END_REF] with p(t) = p 0 f (t/t 0 ) where t 0 is a characteristic timescale for the pressure as in equation ( 4),

q(t) = p 0 R e -Rt/L Rt/L 0 e τ f L Rt 0 τ dτ, (12)
an asymptotic expansion for t 0 ≪ L/R can be obtained by repeated integration by parts:

q(t) ∼ p 0 t 0 L ∞ ∑ k=0 Rt 0 L k e -Rt/L i k+1 f (0) -i k+1 f t t 0 ( 13
) where i 0 f ≡ f and for k = 1, 2, . . . the iterated integrals are

i k f (θ) ≡ ∞ θ i k-1 f (θ) dθ. ( 14
)
The final decay of the flow-rate is

q(t) ∼ p 0 t 0 R e -Rt/L ∞ ∑ k=0 Rt 0 L k i k+1 f (0). (15)
Cauchy's formula for iterated integrals [21, p. 5]:

k!p 0 i k+1 f (0) = p 0 ∞ 0 θ k f (θ) dθ ≡ m k t k+1 0 ( 16
) so that q(t) ∼ Π L e -Rt/L ∞ ∑ k=0 M k k! Rm 1 Lm 0 k , (17)
where the dimensionless small parameter Rm 1 /Lm 0 of the asymptotic expansion is the ratio of the pressure and hydraulic timescales and the moments are normalized as

M k ≡ m k /m 0 (m 1 /m 0 ) k . (18)
Equation [START_REF] Kyser | Design of an impulse ink jet[END_REF] shows that the higher temporal moments of the pressure input increase the later flow-rate relative to the impulsively driven equation ( 9), exactly as seen in figures 3 and 4.

For Asai's λ pressure equation ( 4)

M k = Γ k+1 λ Γ 1 λ Γ 1 λ Γ 2 λ k . (19)
With this, equation ( 17) is compared with the interFoam simulations of figure 3 at 10 µs in figure 5.

Conclusions

The impulsively generated incompressible multiphase flow can be computed by solving a Laplace equation with variable coefficient. This can be used to drive Asai's thermal ink-jet model without having to approximate the Dirac-δ pressure impulse with a bounded drawn-out function; however, it turns out that these finite pressure functions introduced for supposed numerical convenience are not separated in timescale from the hydraulics of typical devices. The flow does approach the impulsively generated one as the pressure-time tends to zero but for typical finite durations, there is an enduring effect. This can be explained in terms of a linear first-order inertance-resistance response which shows how the higher moments of the pressure pulse with respect to time affect the later flow.

Figure 1 :

 1 Figure 1: Sketch of Asai's [3] 'prototype bubble jet printer'. It is filled with ink from to the left of the inlet through to the nozzle, to the right of which is air. The top of the heater is 30 µm × 150 µm.

Figure 2 :

 2 Figure 2: Elevation in the plane of left-right symmetry showing ink (grey) and airy 'bubble' (white) after 30 µs, simulated with interFoam; cf. [3, figure 3]

Figure 3 :Figure 4 :

 34 Figure 3: Bubble expansion, with the same pressure impulse applied (solid) instantaneously, solving equations (2) and (7) for the initial velocity and (dashed) over time by equation[START_REF] Asai | One-dimensional model of bubble growth and liquid flow in bubble jet printers[END_REF]

Figure 5 :

 5 Figure5: Excess of flow-rate at 10 µs driven by finite pressure pulse relative to impulsively driven value; comparison with the first three terms of the asymptotic equation[START_REF] Kyser | Design of an impulse ink jet[END_REF]

Acknowledgements

We thank Dr Paul Croaker of ESI Pacific for much assistance with a previous effort at this problem and Dr Darrin Stephens of Applied CCM for expert advice in programming OpenFOAM.


Christian Frisson

Colin Gallacher

Marcelo M Wanderley

Haptic

Haptic techniques for browsing sound maps organized by similarity

come

M C

Our research lies in a subset of the field of haptics for information visualization [ ] focusing on multimedia data, in our case sound material. Sound designers either record foley sounds or rely on pre-recorded material to compose soundtracks. In both cases, they need to search and browse their collection of sounds during their composition. Stateof-the-art solutions for sound collection browsing combine text-based search using queries on tags and metadata, with content-based exploration of sounds displayed in maps using audio features extracted from the sounds, so as for instance to display files that sound similar in the same visual neighborhoods. With this demo, we expore the design space on how haptic feedback can support sound collection browsing and sequencing.

I D

We propose two haptic techniques for browsing sound maps organized by similarity. ) We enhance the exploratory navigation of sound maps by adding physical effects, such as: magnetizing the pointer to its closest sound item, rendering viscosity when hovering sounds. ) We assist loop-based musical creation by having user-defined sound paths followed by the force-feedback pointing device.

I C

Our system reuses off-the-shelf opensource and openhardware technologies. For organizing sound samples maps by content-based similarity, we forked the Freesound Explorer [ ] that we can run off-the-grid with local servers. For providing force feedback interaction, we used the * christian@frisson.re -supported by Mitacs Elevate grant IT † colin@haply.co ‡ marcelo.wanderley@mcgill.ca

  Haply development kit [ ] set up with the -DOF Pantograph mechanism [ ]. For rendering force feedback effects, we used the Chai d library [ ].

Figure : Freesound

 : Figure : Freesound Tracker with haptic controls tab and sound path