Report on Calculation of Jacobian Matrix of Poincaré Return Map for

Self-synchronized 3D Walking Gaits

Qiuyue Luo 1 , Victor De-León-Gómez 1 , Anne Kalouguine 1,2 , Christine Chevallereau 1 , and Yannick Aoustin 1

Fig. 1: A simplified model of a 3D biped robot. Each leg is massless and has variable length. At the contact point, the stance leg rotates passively around axes x and y, the rotation around axis z is not considered since this rotation is usually inhibited by friction in normal biped locomotion. The swing leg has a fully actuated spherical joint with respect to the concentrated mass of the hip.

I. INTRODUCTION

A. Modeling of walking via LIP model

For simplification, the 3D linear inverted (3D LIP) model is used in this paper to analyze the dynamics of a humanoid robot. In this paper, the robot is approximated as a point mass with point feet and a CoM trajectory that is constrained to a horizontal plane for 3D walking.

In Fig. 1, a simplified model of a 3D biped robot is illustrated. The lines connecting the two feet to the CoM are regarded as the two legs of the simplified model. The stance leg is free to rotate about axes x and y at the ground contact and the length of each leg can be modified through actuation, allowing a desired vertical motion of the pendulum to be obtained.

The configuration of the robot is defined via the position of the CoM (x m ; y m ; z m ) with respect to the reference frame attached to the stance foot and the position of the swing foot denoted by (x s ; y s ; z s ). In order to explore simultaneously the synchronization and stability of periodic orbits for many step length and width, a dimensionless dynamic model of the pendulum will be used. The normalized scaling factors applied along axes x and y depend on the desired step length S and desired step width D. Thus, a new set of viariables is defined:(X, Y, z, X s , Y s , z s ) = ( xm S , ym D , z m , xs S , ys D , z s ). For a pendulum with a constant height of CoM z = z m , the motions in the sagittal and frontal planes are decoupled. Thus, the equations of motion of the 3D LIP with respect to the reference frame attached to the stance foot are [START_REF] Kajita | The 3d linear inverted pendulum mode: A simple modeling for a biped walking pattern generation[END_REF]:

Ẍ = ω 2 X Ÿ = ω 2 Y ( 1
)
where ω = g zm characterizes the LIP and varies with the height of the CoM. As the legs of the robot are assumed to be massless, the swing leg motions X s , Y s , z s do not affect the equation of dynamic of the 3D LIP. The solution to this system is:

X(t) = X + cosh(ωt) + Ẋ+ ω sinh(ωt) Y (t) = Y + cosh(ωt) + Ẏ + ω sinh(ωt) Ẋ(t) = - ωX + 2 sinh(ωt) + Ẋ+ cosh(ωt) Ẏ (t) = ωY + 2 sinh(ωt) + Ẏ + cosh(ωt) (2)
where X + and Y + denote the initial position of the CoM in x direction and y direction respectively during a step, while Ẋ+ and Ẏ + denote the initial velocity if it.

The orbital energies [START_REF] Kajita | Dynamic walking control of a biped robot along a potential energy conserving orbit[END_REF]:

E x = Ẋ2 -ω 2 X 2 E y = Ẏ 2 -ω 2 Y 2 (3)
and the synchronization measure

L = Ẋ Ẏ -ω 2 XY (4)
are conserved during a single support phase [START_REF] Razavi | Restricted discrete invariance and self-synchronization for stable walking of bipedal robots[END_REF]. We can say that the solution in one step is synchronized if and only if the synchronization measure is zero. In fact, this condition L(X, Y, Ẋ, Ẏ ) = 0 defines a one-dimensional submanifold. Any solution starting from this submanifold is synchronized and leads to periodic motion.

B. Transition between support.

Due to the hypothesis that the contact between the swing foot and the ground does not affect the velocity of the CoM, the velocity of CoM will be conserved at each transition of stance leg. Since the reference frame is always attached to the stance foot and the y axis is directed toward the CoM, the sign of the velocity along y axis will be changed from positive to negative [START_REF]Symmetry in legged locomotion: a new method for designing stable periodic gaits[END_REF], i.e.

Ẋ+ k+1 = Ẋ- k , Ẏ + k+1 = -Ẏ - k .
(5)

The state before the transition, i.e. at the end of a step, is expressed by superscript -and that after the transition, i.e. at the beginning of a step, is expressed by + . The variables corresponding to the step k, are denoted with index k, while those of the next step are denoted with k + 1.

After transition, the swing foot placement becomes the new stance foot placement. Thus the CoM position after transition along x axis equals to the CoM position before transition minus the swing foot position. Similar result can be obtained for the CoM position along y axis:

X + k+1 = X - k -X - s,k Y + k+1 = -Y - k + Y - s,k (6)
Knowing the final state of the single support phase, the transition model ( 5) and ( 6) determines the initial state of the ensuing single support phase.

C. Cyclic motion.

For a normalized system, periodic symmetric motions varies from -1/2 to 1/2 along x axis and from 1/2 to 1/2 along y axis. That is:

X * + = -0.5 X * -= 0.5 Y * + = 0.5 Y * -= 0.5 (7)
The superscript * denotes the cyclic motion. Since the orbital energy is constant during one step, the relationship between the final velocity and the initial velocity for a cyclic motion is: Considering the relationship between the final velocity and the initial velocity for cyclic motion described in equation(8), the initial velocity of CoM for a periodic motion of period T can be pointed out from equation(2):

Ẋ * -= Ẋ * + Ẏ * -= -Ẏ * + (8) -0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.5 X 0 0.2 0.4 0.6 0.8 1 Y T = 0.1 T = 1.0
Ẋ * + = ω 1 + cosh(ωT ) 2 sinh(ωT ) Ẏ * + = ω 1 -cosh(ωT ) 2 sinh(ωT ) (9)
In normalized variables, the cyclic motion for different values of time duration T is presented in Fig. 2.

We characterize the orientation of the velocity at the end of the periodic single support phase by α = Ẏ * -Ẋ * -. For cyclic motion in normalized coordinates: α = -1-cosh(ωT ) 1+cosh(ωT ) , 0 < α < 1. During a step, the motion should start with an initial velocity that Ẋ+ > 0 and Ẏ + < 0.

For the 3D LIP model, the eigenvalues of the Poincaré map at the fixed points are {λ λ λ, 1}, where λ λ λ is the set of eigenvalues except of the one respect to kinetic energy. And if ∀|λ i | ∈ λ λ λ < 1, the symmetric periodic orbits are selfsynchronized. This means that the period of oscillations in the x direction eventually matches that in the y direction, and the 3D LIP biped follows a periodic orbit. The other eigenvalue, which is 1, corresponds to neutral stability in kinetic energy. That is, if a small perturbation is applied to the 3D LIP, it will still become synchronized but will eventually follow a periodic orbit with a different level of kinetic energy. For stability, all the eigen values must be strictly less than one in norm.

D. The swing foot motion

A normalized variable Φ varying from zero to one, named phasing variable needs to be defined to describe the trajectory of the swing foot.

The position where to land the swing foot must be chosen cautiously. In this paper, in order to analyze all the possible cases in between, the pose of the swing foot at the end of a step is expressed in a generalized form:

X - s,k = (1 -k S )(X - k -0.5) + 1 Y - s,k = (1 -k D )(Y - k -0.5) + 1 (10)
where 0 ≤ k S ≤ 1 and 0 ≤ k D ≤ 1. How the parameters k S and k D affect the placement of swing foot is illustrated in Fig. 3. With this choice, we can deduce the expression of the error in position of the CoM (δX = X -X * ) through the change of support:

δX + k+1 = k S δX - k δY + k+1 = -k D δY - k (11)
The case k S = k D = 0 allows to nullify the error in position through the change of support, while the case k S = k D = 1 corresponds to a desired fixed step length and width independently of the position of the CoM at the end of the step.

E. Switching condition

In this paper, we chose that the robot accepts to change of support when the CoM crosses the switching manifold S defined to be a line as function of C, as shown in Fig. 4.

S = {(X, Y )|(X -0.5) + C(Y -0.5) = 0} (12)
Many other sets of position can be considered but since stability studied in this paper is a local condition, a line is a convenient choice. In presence of perturbation, the final CoM position error satisfies:

δX - k + CδY - k = 0 (13)
For a chosen evolution of the swing leg, the initial position of CoM can be deduced.

II. CALCULATION OF JACOBIAN MATRIX FOR GENERAL

CASES

To study the stability of the walking gait, the Poincaré return map is used. Since X - k and Y - k are coupled via the switching manifold (12), and the dynamic is autonomous (the time does not appear), the chosen independent state variables are (X - k , L - k , K - k ). The Jacobian is defined for a given periodic motion, and is characterized by the direction of the velocity of the CoM α.

Due to the existence of disturbance, the initial state of the robot after step k is written as:

X + k = X * + + δX + k Y + k = Y * + + δY + k Ẋ+ k = Ẋ * + + δ Ẋ+ k Ẏ + k = Ẏ * + + δ Ẏ + k (14)
where

X 0 = -1 2 , Y 0 = 1 2 .
At the end of the step, the state of the robot is denoted as:

X - k = -X * -+ δX - k Y - k = Y * -+ δY - k Ẋ- k = Ẋ * -+ δ Ẋ- k Ẏ - k = -Ẏ * -+ δ Ẏ - k ( 15
)
Using the fact that Ẋ0 and Ẏ0 define a synchronized motion and neglecting the second order terms, we obtain:

L k = δ Ẋk Ẏ0 + Ẋ0 δ Ẏk - ω 2 2 (δX k -δY k ) (16)
We will now express the final error in velocity as function of the initial error for the step k. As the orbital energies, E x and E y and synchronizaiton measure L are conserved quantities, we have:

( Ẋ- k ) 2 -ω 2 (X - k ) 2 = ( Ẋ+ k ) 2 -ω 2 (X + k ) 2 (17) ( Ẏ - k ) 2 -ω 2 (Y - k ) 2 = ( Ẏ + k ) 2 -ω 2 (Y + k ) 2 (18) Ẋ- k Ẏ - k -ω 2 X - k Y - k = L - k ( 19
)
Submit equation ( 14) and (15) into equation ( 17), ( 18) and ( 19) and neglect the second order terms, we can obtain:

(δ Ẋ- k ) 2 + 2δ Ẋ- k Ẋ0 + Ẋ2 0 - ω 2 4 -ω 2 δX - k -ω 2 (δX - k ) 2 = (δ Ẋ+ k ) 2 + 2δ Ẋ+ k Ẋ0 + Ẋ2 0 - ω 2 4 + ω 2 δX + k -ω 2 (δX + k ) 2 (20) (δ Ẏ - k ) 2 -2δ Ẏ - k Ẏ0 + Ẏ 2 0 - ω 2 4 -ω 2 δY - k -ω 2 (δY - k ) 2 = (δ Ẏ + k ) 2 + 2δ Ẏ + k Ẏ0 + Ẏ 2 0 - ω 2 4 -ω 2 δY + k -ω 2 (δY + k ) 2 (21)
L - k = Ẋ0 δ Ẏ - k -δ Ẋ- k Ẏ0 - ω 2 2 (δX - k + δY - k ) (22)
Using these equations and neglecting the second order terms, we obtain:

δ Ẋ- k = ω 2 2 Ẋ0 (δX - k + δX + k ) + δ Ẋ+ k (23) δ Ẏ - k = ω 2 2 Ẏ0 (-δY - k + δY + k ) -δ Ẏ + k (24)
Submit equation ( 23) and (24) into equation ( 22) and rearrange it, we obtain:

(1 - 1 α )δY - k + (1 -α)δX - k = αδX + k - δY + k α - 4L k ω 2 -δX + k + δY + k (25)
Because of the fact that:

δX - k = -CδY - k (26) δX + k = k S C k D δY + k (27)
Thus replace δY with δX in equation ( 25) and rearrange it, we obtain:

(-1 + α)(1 + αC) α δY - k = (-1 + α)(αCk S + k D ) αk D δY + k - 4Li ω 2 (28)
From the geometric relationship, we know that:

δX + k+1 = k S δX - k (29) δY + k+1 = -k D δY - k (30)
Thus, from equation (28), we can obtain:

-[1 - 1 α -C(1 -α)]δX + k+1 = ( k S C k D α - 1 α - k D C k D + 1)δX + k - 4Lik S C ω 2 (31)
and then:

δX + k+1 = - k D + αCk S 1 + αC δX + k + 4αCk S ω 2 (-1 + α)(1 + αC) L k
(32) From equation (28), we know that:

δX - k = -C( -4αL k ω 2 (-1 + α)(1 + αC) + k D + αCk S k D + αCk D δY + k (33) δY - k = -4αL k ω 2 (-1 + α)(1 + αC) + k D + αCk S k D + αCk D δY + k (34)
Submit equation ( 33) and (34) into equation ( 23) and (24), we obtain:

δ Ẋ- k = δ Ẋ+ k + 2CαL k Ẋ0 (-1 + α)(1 + αC) + ω 2 δX + k 2 Ẋ0 - Cω 2 (k D + αCk S )δY + k 2(k D + αCk D ) Ẋ0 (35) δ Ẏ - k = -δ Ẏ + k + 2αL k Ẏ0 (-1 + α)(1 + αC) + αω 2 C(k D -k S ) 2(1 + αC)k D Ẏ0 δY + k ( 36
)
For the step k + 1, the synchronization measure is:

L - k+1 = Ẋ0 δ Ẏ - k+1 -δ Ẋ- k+1 Ẏ0 - ω 2 2 (δX - k+1 + δY - k+1 ) (37)
Submit equation ( 35) and ( 36) into (37) and rearrange it, we obtain the relationship between L k+1 , L k and δX + K :

L k+1 = (1 + α)(αC -1) -2α(k D -Ck S ) (1 -α)(αC + 1) L k + [(-2 + C)k D + Cl](k D + αCk S )ω 2 2(1 + αC)Ck S δX + k (38)
In conclusion, the Jacobian matrix is:

J =    -k D +αCk S (k S +αC) 4αCk S ω 2 (1+αC)(α-1) 0 J 21 2α(k D -Ck S )+(1-αC)(α+1) (1+αC)(α-1) 0 * * 1    (39) with J 21 = (k D +αCk S )(-1+C+Ck S -k D )ω 2 2(1+αC)Ck S . III. STABILITY STUDY. A. Particular case : k D = k S = 0
In order to nullify the error in pose for the CoM at the beginning of step, k D and k S are supposed to be zero. Thus the Jacobian becomes

J =    0 0 0 α(-1+C)ω 2 2(1+αC) (-1+αC)(α+1) (1+αC)(1-α) 0 * * 1    (40)
It has to be noted that with the proposed control strategy (the initial error in position of CoM is null), the decrease of synchronization measure, at the first order, is independent of the direction of the velocity error.

The eigen values are λ 1 = 0, λ 2 = 1 and λ 3 = (-1+αC)(α+1) (1+αC)(1-α) . The term (1+α) (1-α) is always greater than 1 for 0 < α < 1, its diverging effect increases when α increases. A condition on C to have synchronization can be easily deduced to meet the condition:

-1 < (-1 + αC)(α + 1) (1 + αC)(1 -α) < 1 (41)
which gives:

1 < C < 1 α 2 (42)
The condition λ 3 < 1 ensures convergence toward a synchronized motion, thus is a condition of self-synchronization.

B. Particular case : k

D = k S = 1
When k D = k S = 1, the step width and length are not modified. The Jacobian becomes

J =    -1 4αC ω 2 (1+αC)(α-1) 0 (-1+C)ω 2 C 2α(1-C)+(1-αC)(α+1) (1+αC)(α-1) 0 * * 1    (43)
The eigen values of the Jacobian are λ 1 = 1 and λ 2 , λ 3 =

(1+α)(1-αC)±2 √ α(-1+C)(-1+α 2 C) (1+αC)(α-1)
. When the eigen-value are complex, thus when α(-1 + C)(-1 + α 2 C)(α + 1) < 0 , it can be shown that the norms of the eigenvalues are strictly equal to 1, in the other cases at least one eigen value has a norm greater than 1. No synchronization can be achieved.

Proof: When α(-1 + C)(-1 + α 2 C)(α + 1) < 0, the square of the norm of the second eigenvalue equals to (λ 2 ) 2 = ((1+α)(1-αC)) 2 -4α(-1+C)(-1+α 2 C) In this case, since the step width is not modified, the robot walks along the axis x without deviation.

The Jacobian becomes In this case, the error in x direction at the beginning of a step is nullified, there is no connection between X + k and X + k+1 . Thus instead of X, Y is considered in the Jacobian matrix and it needs to be calculated specifically. (46) For the case where k S = k D = 0, it is very simple to analyze the synchronization because the Jacobian of Poincaré return map is a triangular matrix. However, the feet will not converge to their initial position. On the other hand, if we chose a constant position for landing the swing foot, the synchronization or stability is hard to achieve. For all these cases, the eigenvalues of the Jacobian matrix can be calculated symbolically.

J =    -1+αCk S

Fig. 2 :

 2 Fig. 2: Periodic motions in normalized variables for several values of T , orientation of the final velocity is defined by α.

Fig. 3 :

 3 Fig. 3: Influence of k S and k D on the foot locations. a) Step length and width are fixed; b) The initial CoM position error is nullified; c) The general case. The black and the red dots represent respectively the stance feet during the current and the next steps. The curved line represents the CoM trajectory, and the cross the CoM position at the end of the current step.

Fig. 4 :

 4 Fig. 4: The step is stopped when the CoM crosses the switching manifold (X -0.5) + C(Y -0.5) = 0, which is shown by the red line in the figure. The dashed line is the cyclic motion.

  case : k D = 1, k S < 1

  whereJ 21 = (-2+C+Ck S )(1+αCk S )ω 2 2C(1+αC)k S . D. Particular case : k D ≤ 1, k S = 0

  k D (ω 2 +αω 2 ) (-1+α)(1+αC)ω 2 4αk D (-1+α)(1+αC)ω 2 0 (-1+C-k D )ω 2 2+2αC 1+α-α(1+α)C+2αk D (1+αC)(α-1)


Christophe Abraham

Rémi Servien
email: remi.servien@inra.fr

Nicolas Molinari

A Clustering Bayesian Approach for Multivariate Non-Ordered Circular Data

Keywords: Circular data, Dirichlet process, Non-ordered multivariate data, Projected Normal Distribution, Radiotherapy machine data, Unsupervised clustering

This paper presents a Bayesian model for the clustering of non-ordered multivariate directional or circular data. The particular trait of our data is that each single observation is made up of k ≥ 2 non-ordered points on the circle. We introduce a hierarchical model that combines a symmetrization technique, Projected Normal distributions and a Dirichlet Process. One parameter is introduced to model the non-ordered trait and another one to control the variability of the angles on the circle. An informative prior on the relative locations of the k angles is also provided. The gain of the symmetrization is highlighted by a theoretical study. The parameters of the model are then inferred using a Metropolis-Hastings within Gibbs algorithm. Simulated datasets are analyzed to study the sensitivity to hyperparameters. Then, the benefits of our approach are illustrated by clustering real data made up of the positions of five separate radiotherapy x-ray beams on a circle.

Introduction

Circular and directional data arise in a number of different fields such as oceanography (wave direction), meteorology (wind direction), biology (animal movement direction). The present paper is motivated by circular data in medicine. Nowadays, intensity-modulated radiation therapy (IMRT) has demonstrated its effectiveness for cancer treatment. The latest generation of radiotherapy machines projects multiple rays. Multiplying beams allows to concentrate radiation on the tumor while avoiding the massive irradiation of healthy areas. However, the selection of the incident angles of the treatment beams may be a crucial component of IMRT planning. Due to variations in tumor locations, size and patient anatomy, repositioning for the multiple beams takes a long time and is based on the planner's experience to find an optimal set of beams. So, establishing a small set of standardized beam bouquets for planning could be of valuable help. The set of beam bouquets could be determined by learning the beam configuration features from previous IMRT datasets. The multiple beams are fixed on a circle in the transverse plane around the patient. Consequently, an observation is composed of the -1.0 0.0 1.0 -1.0 0.5

Patient 1 q q q q q -1.0 0.0 1.0 -1.0 0.5

Patient 2 q q q q q -1.0 0.0 1.0 -1.0 0.5

Patient 3 q q q q q -1.0 0.0 1.0 -1.0 0.5

Patient 4 q q q q q -1.0 0.0 1.0 -1.0 0.5

Patient 5 q q q q q -1.0 0.0 1.0 -1.0 0.5 Patient 6 q q q q q -1.0 0.0 1.0 -1.0 0.5 Patient 7 q q q q q -1.0 0.0 1.0 -1.0 0.5 Patient 8 q q q q q -1.0 0.0 1.0 -1.0 0.5 Patient 9 q q q q q -1.0 0.0 1.0 -1.0 0.5

Patient 10 q q q q q -1.0 0.0 1.0 -1.0 0.5

Patient 11 q q q q q -1.0 0.0 1.0 -1.0 0.5

Patient 12 q q q q q -1.0 0.0 1.0 -1.0 0.5

Patient 13 q q q q q -1.0 0.0 1.0 -1.0 0.5 Patient 14 q q q q q Figure 1: Real data set of 14 patients with k = 5 angles. A point on the circle represents the location of a treatment beam.

k beams of a patient, that is k circular measurements. A real data set from post-operative treatment of liver cancer at the Institute of Sainte Catherine in Avignon, France, is represented in Figure 1. One actual observation consists of a (non-ordered) set of k angles rather than of a vector (ordered) of length k but to cope with the technical difficulty of dealing with sets, it is convenient to store the angles of each patient in a vector in increasing order (or in any other given order). Of course, the derived vectors may be very different even for similar sets of angles. This is easily seen by considering a simple case of two patients with angles {1 • , 60 • , 100 • , 150 • , 180 • } and {60 • , 100 • , 150 • , 180 • , 359 • }: the two patients should share the same cluster as the sets of angles are very similar (modulo 360) although the derived vectors are very different and, if any classical clustering method was applied, are not likely to share the same cluster. [START_REF] Abraham | Unsupervised clustering of multivariate circular data[END_REF] proposed a first tool to assist the selection of beam orientations to enhance the therapist's experience. A suitable distance on the circle was defined and, for a fixed number of clusters, an algorithm based on simulated annealing was proposed. [START_REF] Yuan | Standardized beam bouquets for lung IMRT planning[END_REF] generalized the precedent approach using k-medoids to cluster beam configuration features with different numbers of beams. These methods suffer from some major flaws. First, the number of clusters has to be supplied by the user. An additional procedure of model selection (AIC, BIC, RIC, silhouette index, ...) can be used to select the number of clusters but an appropriate methodology that automatically finds this number would be very useful. Second, the final result is only one unique clustering whereas there are probably other clusterings that could be acceptable. A final result with all possible clusterings and a probability of appearance for each could be of great help for the practitioner. These problems can naturally be solved with a Bayesian clustering method based on Dirichlet Process as it does not require a preselected number of clusters and provides different clusterings (possibly with different numbers of clusters) with their posterior probabilities. Also note that the Bayesian framework is well adapted to our application as the sample size is low and can be compensated to some extent by prior information.

Circular data have first been studied using classical non-Bayesian approaches. Three main models for circular data can be found in the litterature: the von-Mises distributions, the wrapped distributions and the projected normal distributions. The von-Mises distributions, first introduced by Von Mises [1918] and extended by [START_REF] Singh | Probabilistic model for two dependant circular variables[END_REF] and [START_REF] Mardia | A multivariate von Mises distribution with applications to bioinformatics[END_REF], are the natural analogues of the normal distribution on the sphere. The wrapped distributions [START_REF] Mardia | Directional Statistics[END_REF] are based on a simple fact that a probability distribution on a circle can be obtained by wrapping a probability distribution defined on the real line. Projected normal distributions are obtained by projecting multivariate normal random variables radially onto the sphere [START_REF] Presnell | Projected multivariate linear models for directional data[END_REF]]. These latter distributions allow for asymmetric and possible bimodal models. We refer the reader to [START_REF] Mardia | Directional Statistics[END_REF] for a complete review on probability distributions of circular data.

Bayesian litterature on circular data is more recent. Von Mises distributions are used in the univariate case in [START_REF] Damien | A full bayesian analysis of circular data using the von Mises distribution[END_REF] and are applied to a change-point problem in SenGupta and Laha [2008]. Wrapped distributions appear in [START_REF] Ravidran | Bayesian analysis of circular data using wrapped distributions[END_REF], with a data augmentation algorithm to overcome some computational difficulties, and in [START_REF] Jona-Lasinio | Spatial analysis of wave direction data using wrapped gaussian processes[END_REF], to handle structured dependences between spatial measurements. Nuñez-Antonio and Gutiérrez-Peña [2005], [START_REF] Wang | Directional data analysis under the general projected normal distribution[END_REF] adapted the projected normal distributions in a Bayesian framework. A more sophisticated model was considered in [START_REF] Wang | Modeling space and space-time directional data using projected gaussian processes[END_REF] to capture structured spatial dependence for modeling directional data at different spatial locations. This model was then upgraded to capture joint structured spatial and temporal dependence [START_REF] Wang | Joint spatio-temporal analysis of a linear and a directional variable: space-time modeling of wave heights and wave directions in the adriatic sea[END_REF]. Then, it was extented to the important spherical case and to any dimension [START_REF] Hernandez-Stumpfhauser | The general projected normal distribution of arbitrary dimension: Modeling and bayesian inference[END_REF]. Also, it was adapted to a multidimensional time series forecasting framework coupled with a Dirichlet process by [START_REF] Mastrantonio | Distributions-oriented wind forecast verification by a hidden markov model for multivariate circular-linear data[END_REF].

Note that, for all the models cited above, each observation is simply a point on a circle or on a sphere while in our case, a single observation is made up of k (k ≥ 2) non-ordered points on the circle. For this reason these models cannot straightforwardly be adapted to our dataset. We propose an extension of the projected normal distribution to our data. This extension does not reduce to a simple projection of a multivariate normal distribution but enables us to model the multivariate and the non-ordered features of our data. We also provide an informative prior distribution on the relative locations of the k angles on the circle. This prior distribution expresses that the k angles are a priori regularly spaced on the circle. A new parameter is also introduced to control the variability of the angles on the circle. Inference on the variabiliy of the angles is of particular interest for a clustering purpose as an inadequate value of this parameter can alter the final results. The projected normal distribution is then associated with a Dirichlet process to perform clustering. Therefore, the proposed method includes an automated selection for the number of clusters.

In the present paper, the Bayesian model is described in the next section. Section 3 is devoted to the inference of the parameters of the model. Section 4 provides a theoretical study to highlight the adaptability of our model to the multivariate and the non-ordered features of the data. Section 5 and 6 provide empirical results first on simulated data and then on the real data set that motivated the present work. A short conclusion is given in Section 7.

Model

A simple way of generating distributions on the p-dimensional unit sphere S p is to radially project probability distributions originally defined on the p-dimensional space R p [START_REF] Presnell | Projected multivariate linear models for directional data[END_REF]]. Let x be a random p-dimensional vector, then x/||x|| is a random point on S p . If x has a p-variate Normal distribution N p (µ, Σ) then x/||x|| is said to have a projected normal distribution, denoted by P N p (µ, Σ). The literature has been first confined to the special case where p = 2 and Σ = I [START_REF] Presnell | Projected multivariate linear models for directional data[END_REF][START_REF] Nuñez-Antonio | A bayesian analysis of directional data using the projected normal distribution[END_REF][START_REF] Nuñez-Antonio | A bayesian regression model for circular data based on the projected normal distribution[END_REF]. Then, [START_REF] Wang | Directional data analysis under the general projected normal distribution[END_REF] studied the projected normal family with a general covariance matrix Σ and refer to this richer class P N p (µ, Σ) as the general projected normal distribution. This general version allows asymmetry and bimodality [see Figure 2. in [START_REF] Wang | Modeling space and space-time directional data using projected gaussian processes[END_REF]. The general projected normal distribution is not identifiable because x/||x|| is invariant to scale transformation. To overcome this problem [START_REF] Wang | Directional data analysis under the general projected normal distribution[END_REF] fixed some variance parameters in Σ to provide identifiability.

In a first step of simplification, we assume that the ith of the n observations is given by a vector of k angles θ i = (θ i1 , . . . , θ ik ) ∈ [0, 2π[ k instead of a non-ordered set {θ i1 , . . . , θ ik }. Using a projected normal distribution, we denote by x i = (x i1 , . . . , x ik ) ∈ (R 2 ) k a random vector with distribution N 2k (µ i , I 2k ) where θ ij is defined as the radial projection of x ij on the unit circle of R 2 . In other words, we have x ij = (x ij1 , x ij2 ) = (r ij cos θ ij , r ij sin θ ij ) for all i ∈ {1, . . . , n} and all j ∈ {1, . . . , k} where r ij denotes the Euclidean norm of x ij . Note that θ i is observed while r i = (r i1 , . . . , r ik ) is not and is treated as an unknown parameter. We denote by P N 2k (µ i , I 2k ) the joint distribution of (θ i , r i ). Clustering analysis will be based on a Dirichlet process mixture (DPM) model described as follows:

θ i , r i |µ ∼ P N 2k (µ i , I 2k ) µ i |P ∼ P P ∼ DP (n 0 P 0 ), (1)
where µ = (µ 1 , . . . , µ n ) and where DP (n 0 P 0 ) denotes the Dirichlet process (DP) introduced by [START_REF] Ferguson | A bayesian analysis of some nonparametric problems[END_REF] with center P 0 = N 2k (0, Σ 0 ) and precision parameter n 0 . The clustering properties of the DP are well known and date back to [START_REF] Blackwell | Ferguson distributions via Polya urn schemes[END_REF]. It is shown that the parameter µ = (µ 1 , . . . , µ n ) follows the Pólya urn scheme:

µ 1 ∼ P 0 µ i+1 |µ 1 , . . . , µ i ∼ 1 n 0 +i i j=1 δ µ i + n 0 n 0 +i P 0 , for i ≥ 2.
(

) 2
with δ µ i indicating the point measure on µ i . So, µ i+1 may be equal to one of the previous µ i 's or may be drawn from P 0 . This results in a positive probability of sharing the parameter value with previous observations; hence the clusters. In the sequel, we will denote by Pólya(n 0 P 0 ) the distribution of µ given by (2). Although the DPM is very popular for Bayesian clustering, other model-based cluster methods exist. For a review of these methods, we refer the reader to [START_REF] Quintana | A predictive view of bayesian clustering[END_REF], [START_REF] Lau | Bayesian model-based clustering procedures[END_REF], [START_REF] Fritsch | Improved criteria for clustering based on the posterior similarity matrix[END_REF] and references therein.

Note that the DPM model does not require choosing the number of clusters. On the other hand, it is well known that the number of clusters can be controlled by n 0 . Learning about n 0 from the data may be addressed by assuming a Gamma prior distribution n 0 ∼ G(a n 0 , b n 0 ) [START_REF] Escobar | Bayesian density estimation and inference using mixtures[END_REF].

Now recall that the actual ith observation consists of a (non ordered) set of the form {θ i1 , . . . , θ ik } rather than of a vector (ordered) θ i = (θ i1 , . . . , θ ik ) . The impact of this simplification is quite easy to understand. Using model (1), two observations i 1 and i 2 with the same angles but in different orders would have a very low posterior probability of sharing the same cluster, that is µ i1 = µ i2 . We treat the observations as vectors for convenience but we have to introduce a permutation parameter τ i to compensate this simplification. More precisely, for all µ i = (µ i1 , . . . , µ ik ) and all permutation τ i of {1, . . . , k}, we set µ τ i i = (µ iτ i (1) , . . . , µ iτ i (k) ) ; µ τ i i can be viewed as a random permutation of the coordinates of µ i . Therefore, the clustering model becomes:

θ i , r i |µ, τ ∼ P N 2k (µ τ i i , I 2k ) µ i |P ∼ P P ∼ DP (n 0 P 0 ), (3)
where τ = (τ 1 , . . . , τ n ) and µ = (µ 1 , . . . , µ n ). The permutations τ i are assumed to be a priori independent with a uniform distribution U P on the set P of permutations of {1, . . . , k}. The posterior probability that two observations i 1 and i 2 with the same angles but in different orders would share the same cluster is increased with model (3) as there exist some values of τ i1 and τ i2 such that µ τ i1 i1 = µ τ i2 i2 . A theoretical study of the impact of the symmetry introduced by τ i is given in Section 4.

Prior information

It is natural to assume that the k angles θ i1 , . . . , θ ik are a priori distributed so that the radial projections x i1 / x i1 , . . . , x ik / x ik are roughly equally spaced on the unit circle. As x ij / x ij can be seen as a random change of µ ij / µ ij , this prior information can be incorporated into the covariance matrix Σ 0 of P 0 as follows. From (3), it is well known that the marginal distribution of µ i is P 0 = N 2k (0, Σ 0 ). Denote by R the 2 × 2-matrix of the rotation in R 2 with angle 2π/k and center 0 and by ε ij , j ∈ {2, . . . , k}, k -1 random variables with distribution ε ij ∼ N 2 (0, I 2 ). Assume that µ i1 , ε i2, , . . . , ε ik are independent and set µ i1 ∼ N 2 (0, ρI 2 ) where ρ is a positive number and µ ij = R j-1 µ i1 + ε ij for j ∈ {2, . . . , k}. It is important to note the influence of ρ on the distribution of µ i . First, it is shown in the Appendix that the expected value of µ ij is ρπ/2 for j = 1 and (ρ + 1)π/2 for j ≥ 2 and that the components µ ij become more correlated and equally spaced as ρ increases. Furthermore, if we denote by

v ij = µ ij / µ ij the radial projection of µ ij onto the unit circle, it is shown that (v i1 , . . . , v ik ) tends to (v i1 , Rv i1 , . . . , R k-1 v i1 ) almost surely as ρ → ∞ and tends to (ε i1 / ε i1 , . . . , ε ik / ε ik )
almost surely as ρ → 0 where ε i1 , . . . , ε ik are independent and identically distributed random variables with distribution N 2 (0, I 2 ). In other words, the radial projections of the components of µ i are highly correlated and equally spaced for large values of ρ but approximately independent and uniformly distributed for small ones. The influence of ρ is also studied through some simulations in the sequel (see Subsection 5.1 and Figure 2).

Inference on ρ can be performed using an inverse gamma prior ρ ∼ IG(a ρ , b ρ ) for which the full posterior conditional distribution will be calculated in the following section.

It is worth pointing out that the weak dissymmetry introduced in the definition of P 0 between µ i1 and the other components of µ i leads to a convenient closed-form expression of Σ 0 from which closed-form expressions of the inverse Σ -1 0 and the determinant |Σ 0 | can be obtained as well. Such closed-form expressions and the full posterior conditional distribution of ρ could not have been obtained so easily with a perfect symmetric distribution P 0 . The calculations of Σ 0 , Σ -1 0 and |Σ 0 | are given in the Appendix. To highlight the dependence on ρ, Σ 0 will be also denoted by Σ 0 (ρ) in the sequel.

Finally, the complete Bayesian model can be expressed as follows:

θ i , r i |µ, τ ∼ P N 2k (µ τ i i , I 2k ) µ|n 0 , ρ ∼ Pólya(n 0 P 0 (ρ)) τ i ∼ U P ρ ∼ IG(a ρ , b ρ ) n 0 ∼ G(a n 0 , b n 0 ). ( 4
)
where P 0 (ρ) = N 2k (0, Σ 0 (ρ)). By convention, it is assumed that the random variables at a stage of the hierarchy are independent.

Inference

We set θ = (θ 1 , . . . , θ n ), r = (r 1 , . . . , r n ), µ = (µ 1 , . . . , µ n ), τ = (τ 1 , . . . , τ n ) and ξ = (r, µ, τ, ρ, n 0 ). Thus, the parameter is ξ and the observation is θ. We sample from the posterior distribution of ξ with a Metropolis-Hastings-Within-Gibbs algorithm. In what follows, p stands for a generic notation for a density distribution.

Simulations of µ

We can restrict our attention to model (3) instead of the full model ( 4) for the simulations of µ as every component of ξ except µ remains fixed. An alternative parameter setting of µ, θ and ρ will prove useful. Denote x = (x 1 , . . . , x n ) where x i = (x i1 , . . . , x ik ) . Firstly, note that the full conditional distribution of µ reduces to the conditional distribution of µ given (x, n 0 , ρ, τ ) as there is a natural bijection between x ij and (θ ij , r ij ). Secondly, if we denote by N 2k (x i ; µ i , I 2k ) the value of the density of N 2k (µ i , I 2k ) at x i , it is easy to check that:

N 2k (x i ; µ τ i i , I 2k ) = N 2k (x τ -1 i i ; µ i , I 2k ) (5)
where τ -1 i is the permutation such that x

τ i oτ -1 i i = x i . Consequently, if we set y i = x τ -1 i i
, sampling from the posterior distribution of µ in the DPM model (3) reduces to sampling from the posterior distribution of µ in the following conjugate DPM model:

y i |µ ∼ N 2k (µ i , I 2k ) µ i |P ∼ P P ∼ DP (n 0 P 0 ). (6)
There are several samplers for conjugate DPM models; for a review, we refer the reader to [START_REF] Maceachern | Computational methods for mixture of Dirichlet process models[END_REF], [START_REF] Neal | Markov chain sampling method for Dirichlet process mixture models[END_REF], [START_REF] Griffin | Computational issues arising in bayesian nonparametric hierarchical models[END_REF]. Following the notations of [START_REF] Dahl | An improved merge-split sampler for conjugate Dirichlet process mixture models[END_REF], we use a parameter setting of µ in terms of:

• a set partition η = {S 1 , . . . , S q } for {1, . . . , n} where each S j represents a cluster, i.e., µ i = µ j if there exists j 1 ∈ {1, . . . , q} such that i, j ∈ S j 1 and µ i = µ j if there exist j 1 , i 1 ∈ {1, . . . , q}, i 1 = j 1 such that i ∈ S i 1 , j ∈ S j 1 ,

• a vector φ = (φ 1 , . . . , φ q ) composed of the distinct values of µ, i.e., φ j = µ i for all i ∈ S j .

Then, the conjugate DPM model ( 6) may be expressed as:

y i |η, φ ∼ N 2k ( q j=1 φ j 1 {i∈S j } , I 2k ) φ j |η ∼ P 0 η ∼ p(η) ∝ q i=1 n 0 Γ(|S j |), (7)
where |S j | is the cardinal of S j , 1 A is the indicator function for the event A, Γ denotes the gamma function and p stands for the generic notation for any density. We can integrate over the cluster location parameter φ analytically in (7) as P 0 is conjugate to the normal distribution of y i given η and φ. Then, we run the SAMS sampler of [START_REF] Dahl | An improved merge-split sampler for conjugate Dirichlet process mixture models[END_REF] for simulating η. This sampler may improve the merge-split sampler initially proposed by [START_REF] Jain | A split-merge Markov chain Monte-Carlo procedure for the Dirichlet process mixture model[END_REF]. Once a simulation of η is obtained, it is easy to simulate the cluster location parameter φ from its full conditional which reduces to sample independently each φ j from a N 2k (Σ j i∈S

j y i /|S j |, Σ j ) distribution with Σ -1 j = |S j | -1 I 2k + Σ -1 0 (ρ).
As recommended by the previous authors, we combine three runs of the Metropolis-Hastings update of the SAMS sampler with a full scan of Gibbs sampling for µ [see [START_REF] Maceachern | Estimating normal means with a conjugate style Dirichlet process prior[END_REF], for a presentation of this particular Gibbs sampler]. Some details of the SAMS and the Gibbs samplers used in this article are given in the Appendix.

Simulations of r It is shown in the Appendix that the r ij are independent given (θ, τ, µ, ρ, n 0 ) with density:

p(r ij |θ, τ, µ, ρ, n 0 ) ∝ r ij e -1 2 (rij-u ij µ iτ i (j)) 2 , (8)
with u ij = (cos θ ij , sin θ ij ). If we denote by N + 1 (m, v) the univariate normal distribution truncated to [0, ∞), we remark that (8) is close to the value of the density of N + 1 (u ij µ iτ i (j) , 1) at r ij . It is then natural to simulate from (8) by a Metropolis-Hastings step with a N + 1 (u ij µ iτ i (j) , 1) as the proposal distribution. Clearly, the probability of acceptance reduces to the ratio min{r new ij /r old ij , 1} where r old ij and r new ij are, respectively, the current and the proposed values of r ij in the algorithm. Alternative methods of simulations could have been used at this step such as, for example, the slice sampler proposed by [START_REF] Hernandez-Stumpfhauser | The general projected normal distribution of arbitrary dimension: Modeling and bayesian inference[END_REF].

Simulations of τ As the prior distribution of τ is uniform, we have:

p(τ |θ, r, µ, ρ, n 0 ) = p(τ |x, µ) ∝ p(x|τ, µ) ∝ n i=1 N 2k (x i ; µ τ i i , I 2k ).
Thus, given (θ, r, µ, ρ, n 0 ), the τ i are independent with density (with respect to the counting measure on the set T of permutations of {1, . . . , k}):

p(τ i |x, µ) = N 2k (x i ; µ τ i i , I 2k ) t∈T N 2k (x i ; µ t i , I 2k ) . ( 9
)
Simulations of ρ From (4), it is clear that the full conditional distribution of ρ reduces to the conditional distribution of ρ given µ. Then, using the parametrization of µ in terms of (η, φ) and ( 7), we note that η and ρ are independent and that:

p(ρ|θ, r, µ, τ, n 0 ) = p(ρ|η, φ) ∝ p(φ|η, ρ)p(ρ|η)

∝ q j=1 p(φ j |ρ) p(ρ). ( 10
)
We show in the Appendix that |Σ -1 0 (ρ)| = ρ -2 and that the components of the matrix Σ -1 0 (ρ) are independent (constant) of ρ except the components of the first 2 by 2 diagonal submatrix (lines and columns 1 and 2). As this submatrix is equal to (ρ -1 + (k -1))I 2 , it is easily seen that

φ i Σ -1 0 (ρ)φ i = (ρ -1 + (k -1))φ i1 φ i1 + constant = ρ -1 φ i1 φ i1 + constant.
where constant stands for a generic notation for an expression independent of ρ. Since φ j |ρ ∼ P 0 (ρ) = N 2k (0, Σ 0 (ρ)) and ρ ∼ IG(a ρ , bρ), we have:

q j=1 p(φ j |ρ) ∝ ρ -q e -1 2 ρ -1 q j=1 φ i1 φ i1 ,
and it is easy to conclude from (10) that the full conditional of ρ is

IG a ρ + q, b ρ + 1 2 q i=1 φ i1 φ i1 . ( 11
)
Simulations of n 0 Using the arguments of [START_REF] Escobar | Bayesian density estimation and inference using mixtures[END_REF], under the G(a n 0 , b n 0 ) prior, n 0 is updated at each Gibbs iteration by sampling first an additional variable ζ from a Beta distribution and then a new value of n 0 from a mixture of Gamma distributions as follows:

ζ|n 0 ∼ B (n 0 + 1, n) n 0 |ζ, q ∼ π n G(a n 0 + q, b n 0 -log ζ) + (1 -π n )G(a n 0 + q -1, b n 0 -log ζ), (12)
with weights π n defined by

π n /(1 -π n ) = (a n 0 + q -1)/[n(b n 0 -log ζ)].
The whole procedure is summarized in Algorithm 1.

Theoretical study of the symmetrized model

To investigate the impact of the symetrization induced by the variables τ i , we consider a simple model of the following form:

x i |η, φ ∼ N 2k ( q j=1 φ j 1 {i∈S j } , I 2k ) φ j |η ∼ P 0 η ∼ G (I)
Algorithm 1 Require: Data set θ = (θ 1 , . . . , θ n ).

Require: Hyperparamaters a ρ , b ρ , a n 0 , b n 0 . Repeat :

1. Simulate η.

(a) Run the SAMS sampler three times.

(b) Run the Gibbs sampler.

2. Simulate φ j ∼ N 2k (Σ j i∈S j y i /|Σ j |, Σ j ) for each cluster j.

Propose r

new ij ∼ N + 1 (u ij µ iτ i (j)
, 1), accept with probability min(r new ij /r old ij , 1).

4. Simulate new τ i from 9.

5. Simulate new ρ from 11.

6. Simulate n 0 from 12.

and its symmetrized version :

x i |η, φ ∼ N 2k ( q j=1 φ τ i j 1 {i∈S j } , I 2k ) φ j |η ∼ P 0 η ∼ G τ i ∼ U P , (II)
where φ τ i j = (φ jτ i (1) , . . . , φ jτ i (k) ) is obtained by random permutation of the coordinates of φ j = (φ j1 , . . . , φ jk ) ∈ (R 2 ) k . In both models, P 0 = N 2k (0, Σ 0 ) and G is any distribution of the partition η = {S 1 , . . . , S q } of {1, . . . , n}. Such distributions include the distribution derived from the Dirichlet process given by (7). Model (II) can be viewed as a simplified and reparametrized version of (4). Now consider an idealized sample x 1 , . . . , x n for which every observation x i is simply a random permutation of one unique observation x 0 = (x 01 , . . . , x 0k ) ∈ (R 2 ) k ; in other words, for every i, there exists a permutation α i such that x i = (x 0α i (1) , . . . , x 0α i (k) ) . As the coordinates x ij of all the x i are the same but in a different order, it is expected that all the observations are put together in one unique cluster. The aim of this section is to study whether model (II) is more appropriate than model (I) for this purpose.

Let p 0 and p I (x|η) denote respectively the density of P 0 and the conditional density of x = (x 1 , . . . , x n ) given η for model (I). We have:

p I (x|η) = q j=1 i∈S j N 2k (x i ; φ j , I 2k )p 0 (φ j )dφ j = q j=1 m(x S j ),
where x S j = (x i , i ∈ S j ) and m(x S j ) = i∈S j N 2k (x i ; φ j , I 2k )p 0 (φ j )dφ j .

Denote by p II (x|η) the conditional density of x given η for model (II). By ( 5) and noting that {τ -1 i , τ i ∈ P} = P, we have:

p II (x|η) = τ 1 (k!) n q j=1 i∈S j N 2k (x i ; φ τ i j , I 2k )p 0 (φ j )dφ j = τ 1 (k!) n q j=1 i∈S j N 2k (x τ i i ; φ j , I 2k )p 0 (φ j )dφ j = 1 (k!) n τ q j=1 m(x τ S j ),
where the sum above is taken for all the values of τ = (τ 1 , . . . , τ n ) in P n , x τ S j = (x τ i i , i ∈ S j ) and x τ i i = (x iτ i (1) , . . . , x iτ i (k) ) . Therefore, models (I) and (II) reduce to

x|η ∼ q j=1 m(x S j ) η ∼ G, (I') and x|η ∼ 1 (k!) n τ q j=1 m(x τ S j ). η ∼ G. (II')
For all partition η = {S 1 , . . . , S q } and all observation x, we set

f (x, η) = 1 (k!) n τ ∈P n exp 1 2 q j=1   i∈S j x τ i i 2 S j - i∈S j x i 2 S j   ( 13
)
where Σ S = Σ -1 0 + |S|I 2k -1 for all subset S ⊂ {1, . . . , n} and t 2 S = t Σ S t for all t ∈ (R 2 ) k . Proposition 1. a) For all partition η = {S 1 , . . . , S q } and all observation x = (x 1 , . . . , x n ), we have:

p II (x|η) p I (x|η) = f (x, η).
b) For all distribution G, there exists a positive number B G such that:

p II (η|x) p I (η|x) = B G f (x, η),
for all partition η and all observation x.

c) For all distribution G, all partition η and all observation x, we have:

p II (η|x) p I (η|x) ≥ f (x, η) 1 max η f (x, η) (14)
where the maximum is taken over all partitions of {1, . . . , n}.

From a) of Proposition 1, we see that f (x, η) is the likelihood ratio of models (II') and (I'). From b), we know that the posterior odds ratio is large when f (x, η) is large. It would be of interest to know whether this ratio is greater than one. Unfortunately, this is not an easy task except for a few particular cases given below. Indeed, although the factor B G is actually known (see the proof of Proposition 1 in the Appendix), it is rather intractable. From c), we deduce that the posterior odds is actually greater or equal to one at least for the partition η x that maximizes f (x, η). This partition does exist for any observation x and is independent of G. In other words, for any x, there exists a partition η x such that p II (η x |x) ≥ p I (η x |x) for all prior G. Finally, we can remark from the proof of the theorem that the equality in ( 14) is obtained when G is a Dirac distribution; a meaningless prior.

Consider the partition η with a single cluster: q = 1 and S 1 = {1, . . . , n}. From ( 13), the posterior odds ratio when η = η is likely to be large when n i=1 x i ≈ 0 and small when all the x i ≈ x 0 for all i ∈ {1, . . . , n}. Assume from now that n i=1 x i = 0 and that Σ 0 = I 2k . Remenber that Σ 0 models the prior information about the mutual positions of the angles on the circle. Therefore Σ 0 = I 2k can be viewed as a non informative prior. In this case,

t 2 S j = (1 + |S j |) -1 t t = (1 + |S j |) -1 t for all t ∈ (R 2
) k and we have:

f (x, η) = 1 (k!) n τ ∈P n exp 1 2(n + 1) n i=1 x τ i i 2 . ( 15
)
Example 1 below provides a typical sample x = (x 1 , . . . , x n ) for which the posterior probabilty of a unique cluster is greater with model (II) than with model (I) independently of the prior distribution G.

Example 1. First, by noting that:

n l=0 e ilθ = sin θ(n+1) 2 sin θ 2 e i θ 2 n ,
for all θ ∈ R, we deduce that:

k-1 l=0 cos 2πl k = 0 and k-1 l=0 sin 2πl k = 0. ( 16
)
Assume n = k and set x ij = (cos(i + j -2)2π/k, sin(i + j -2)2π/k) for i ∈ {1, . . . , k} and j ∈ {1, . . . , k}. In other words, x 1 = (x 11 , . . . , x 1k ) ∈ (R 2 ) k is made up of k consecutive points on the unit circle separated from an angle of 2π/k, x 2 is obtained by a rotation with angle 2π/k of each point of x 1 and so on. Therefore, it is easy to see from (16) that n i=1 x i = 0. Our conjecture is that max η f (x, η) = f (x, η) for all integer k which implies, from c) of Proposition 1, that the probability of a unique cluster is greater for model (II) than for model (I) for any distribution G. For n = k = 2 the conjecture reduces to f (x, η) ≤ f (x, η) for a single partition η = {{x 1 }, {x 2 }}. As x i S j = x τ i i S j for all i and τ i , it is easily seen from ( 13) that f (x, η) = 1. On the other hand, as x 1 2 = k and x 1 = -x 2 , we see from (15) that

f (x, η) = 1 4 2 exp 1 6 x 1 + x 2 2 + 2 exp 1 6 2x 1 2 = 1 2 1 + 2 exp 4 3 ,
hence the proof of the conjecture for n = k = 2. We also proved the conjecture for n = k = 3 with a rather large amount of calculations (not given here) to take into account all the partitions η and all the permutation τ = (τ 1 , τ 2 , τ 3 ). We are not in a position to provide general proof of the conjecture for n = k ≥ 4.

Simulations

With small data sets, a misspecification of the prior could have a strong negative impact on the final results. Therefore, special attention has to be paid to the prior and the hyperparameter specifications. Consequently, we test our algorithm on two simulation studies to evaluate the influence of some hyperparameters. The performances of our method are investigated using the Adjusted Rand Index (ARI), proposed by [START_REF] Hubert | Comparing partitions[END_REF], to compare our obtained partition to the actual one. The Rand Index [START_REF] Rand | Objective criteria for the evaluation of clustering methods[END_REF] is a well known measure of the similarity between two partitions. If we denote by N 00 the numbers of pairs that are in the same cluster in both partitions and by N 11 the number of pairs that are in different clusters in both partitions, then the Rand Index is defined by the ratio (N 00 + N 11 )/( n 2 ). The ARI is a corrected-for-chance version of the Rand index. Its expected value (under the generalized hypergeometric model) is equal to 0 and its maximum is 1 while the expected value of the Rand Index depends on the number of clusters. For a presentation of the different criteria for clustering comparison and for a study investigating the usefulness of the adjusted measures, we refer the reader to [START_REF] Fritsch | Improved criteria for clustering based on the posterior similarity matrix[END_REF] and [START_REF] Nguyen | Information theoretic measures for clustering comparison: Is a correction for chance necessary ?[END_REF].

Influence of the Precision Parameter ρ

First we choose to simulate data using a procedure which is close to our model in order to investigate the influence of the precision parameter ρ. We set q = 3 clusters of 10 data. We simulate the coordinates µ ij of each center µ i approximately on a circle with a fixed radius. Since it is shown in Section 2 that E µ i1 = ρπ/2 ≈ 1.25 √ ρ, the first coordinate µ i1 is simulated according to a uniform distribution on the circle with radius 1.25 √ ρ. The other coordinates µ ij , j = 2, . . . , 5 (k = 5) are generated according to a noisy rotation with angle 2πj/5 of µ i1 . For each cluster i, we generate 10 data according to P N 10 (µ i , I 10 ). A comparison of the generated data is provided in Figure 2 with different values for ρ; for the clarity of the picture we choose to represent only q = 2 clusters of 5 observations. It is clear from Figure 2 that large values of ρ provide small variability for the projected observations. This observation can be confirmed by a simulation study. We simulate some datasets according to the above procedure with different values for ρ. For each value of ρ, a hundred datasets are simulated. Then, our Bayesian methodology is applied using a n 0 = 10 and b n 0 = 1. The mean values for the ARI are given in Table 1. As expected, this parameter has an important influence on the obtained results. We choose a noninformative prior for ρ by setting a ρ = b ρ = 0.01. This prior is close to the Jeffreys prior for the model N 2 (0, ρI 2 ) whose density reduces to ρ -1 .

Robustness to the Hyperparameters a n 0 and b n 0

It is well-known that the number of clusters does depend on n 0 whose prior distribution is fixed by the hyperparameters a n 0 and b n 0 . In this subsection we investigate the sensitivity of the ARI with respect to these hyperparameters. We apply the same simulation strategy as in the previous subsection with a fixed value of ρ such that 1.25 √ ρ = 20. Note that the parameters a n 0 and b n 0 are not at all involved in the simulation of the dataset. The mean values for the ARI over 100 simulated data sets are given in Table 2.

Table 2: Adjusted Rand Index (Proportion of clusterings with the actual number of clusters) according to a n 0 and b n 0 .

b n 0 = 0.1 b n 0 = 1 b n 0 = 10 b n 0 = 100 b n 0 = 1000 a n 0 = 0.

Influence of the number of clusters

Using the previously fixed hyperparameters a simulation study is provided with different numbers of clusters. For each number of clusters, one hundred datasets are simulated according to the above procedure. Each cluster has a number of points randomly chosen between 10 and 20. The results are gathered in the following Table 3. The results are rather good even with a high number of clusters.

Real Data

We then apply the methodology to a real data set from post-operative treatment of liver cancer at the Institute of Sainte Catherine in Avignon, France (see Figure 1 and Table 4). Let us recall that no other competing methods exist for these kind of multivariate circular data except the method described in [START_REF] Abraham | Unsupervised clustering of multivariate circular data[END_REF] with a fixed number of clusters. Consequently, our results are compared to those of [START_REF] Abraham | Unsupervised clustering of multivariate circular data[END_REF] in which the number of clusters was preselected to q = 2. Let us remind you that the a priori distribution of n 0 is a gamma distribution with parameter a n 0 and b n 0 with an expected value equal to a n 0 /b n 0 (if a n 0 > 1) and a variance equal to a n 0 /b 2 n 0 (if a n 0 > 2). Remember that the expected number of clusters given n 0 is approximately equal to n 0 log(1+n/n 0 ) [START_REF] Teh | Dirichlet processes[END_REF]. According to the results of Section 5, the results are robust with respect to the choice of the hyperparameters a n 0 and b n 0 with 1 ≤ a n 0 /b n 0 ≤ 10. We choose a rather non-informative prior by setting a n 0 = 3 and b n 0 = 0.3 which leads to a distribution of n 0 centered around 3 with a large variance. Other values for a n 0 and b n 0 have been tested and give nearly the same results. As in Section 5, we choose a non-informative prior by setting a ρ = b ρ = 0.01.

MCMC convergence diagnostics was investigated with the clustering entropy

- q i=1 |S i | n log |S i | n .
Traceplots for this quantity and for other parameters of the model suggest a good mixing and the convergence of our chain. On a classical personal laptop and using a nonoptimized code, the time of mixing was approximately twelve minutes whereas the whole procedure lasts ninety minutes.

The majority clustering (mode of the posterior distribution of the clusterings) is the same as in [START_REF] Abraham | Unsupervised clustering of multivariate circular data[END_REF] (two clusters: one containing data 1,2,6,9 and 12, the second containing data 3,4,5,7,8,10,11,13 and 14) with a posterior probability equal to 30.5%. This result was Figure 3: Barplot of the proportion of 32 more probable clusterings. For a better apprehension, the clusterings that have a posterior probability lower than 0.02% have been omitted. awaited and is coherent with the choice of 2 clusters in the previous method. But the real gain from our Bayesian approach is to look beyond this majority clustering. Here there are 3 more clusterings that are significant and that could give some information on this real dataset. The second majority clustering is nearly the same as the previous one : the clusters are the same but data 6 is alone in a third cluster. Indeed, this data is very atypical because it is the only one that contains an angle near 1.69π. The posterior probability for this clustering is 14.9%. The third majority clustering gives nearly the same information with a posterior probability of 13.5%. There are two clusters: one with data 6 and a second with all the others. Finally, another clustering with a posterior probability of 12.0% is made up of only one cluster. Even with other choices for the hyperparameters a n 0 and b n 0 , the posterior probability of this clustering remains high. It highlights the fact that all the data share some common traits and the main difference in the two clusters of the majority clustering only concerns one angle. All the clusterings are included in Figure 3 sorted by their posterior probabilities. It can be noted that a credible region with a posterior probability of 71% is composed of the 4 previous clusterings.

We give in Figure 4 the posterior distribution of the number of clusters. The posterior probabilities of 1, 2 or 3 clusters are respectively 12%, 65% and 21%. Consequently, the number of clusters is certainly (with probability 98%) less than or equal to 3. As expected, these results are in line with the clusterings obtained in [START_REF] Abraham | Unsupervised clustering of multivariate circular data[END_REF]. The final choice of 2 clusters (that is not made here a priori) could provide, using the centers of the clusters, preset positions for praticians. As explained above, these two centers share only one main difference on one unique angle. This is highlighted by the important posterior probability of the clustering with only one cluster. Thus, using these preset positions should be fairly easy for praticians, with four fixed values and only two choices for the last one. Furthermore, the results suggest another preset position that should be added and tested if the two previous one do not fit: the beam angles of data 6. As explained in [START_REF] Yuan | Standardized beam bouquets for lung IMRT planning[END_REF], the definition of such presets will help to save time (at least 30 minutes for each patient) and will allow more people to be treated with this technology; it is also shown that the beams generated with our methodology show dosimetric qualities comparable to their manually generated clinical counterpart, even if no adjustments were allowed around the fixed presettings.

Conclusion

We present a full Bayesian framework for the clustering of multivariate circular and non-ordered data. It is based on a hierarchical model that combines Projected Normal distributions and the Dirichlet Process. Two original parameters are also introduced in this model: the parameter ρ to infer the variance of the angles and the symmetrization parameter τ to model the nonordered feature of the data. The parameters of the model are then inferred using a Metropolis-Hastings within Gibbs algorithm and a theoretical study of the impact of the symmetrization parameter is provided. The simulation study and the real data example show the benefits of this approach. Indeed, the number of clusters is chosen automatically by the method and the final result is much more complete than the majority clustering which is usually provided by classical clustering algorithms. However some improvements could be considered, such as, incorporating covariates (shape or size of the tumor, stage of the cancer, sex, age, ...) to preselect the beam positions.

Figure 2 :

 2 Figure 2: Two data sets are generated with two different values for the parameter ρ to highlight the influence of this parameter. Two clusters of five data are represented on each plot. Each data is composed of k = 5 angles on the circle. One cluster is represented by the black cross, the other by the red square.

Figure 4 :

 4 Figure 4: Posterior distribution of the number of clusters.

Table 1 :

 1 Adjusted Rand Index according to ρ.

	ρ 1.25 √ ρ	0.0064 0.1024 0.64 5.76 256 0.1 0.4 1 3 20
	ARI	0.35	0.39 0.45 0.62 0.78

  Table 2 suggests that a choice of a n 0 /b n 0 approximately between 1 and 10 provides good and similar results.

	1	0.73 (0.80) 0.71 (0.79) 0.62 (0.72) 0.63 (0.75) 0.59 (0.67)
	a n 0 = 1	0.76 (0.91) 0.72 (0.84) 0.65 (0.79) 0.67 (0.76) 0.64 (0.71)
	a n 0 = 10	0.72 (0.76) 0.78 (0.96) 0.69 (0.84) 0.67 (0.80) 0.65 (0.74)
	a n 0 = 100	0.70 (0.70) 0.68 (0.79) 0.79 (0.92) 0.72 (0.82) 0.62 (0.75)
	a n 0 = 1000 0.66 (0.69) 0.62 (0.72) 0.68 (0.79) 0.75 (0.88) 0.65 (0.76)

Table 3 :

 3 Adjusted Rand Index according to the number of clusters.

	Number of clusters	1	3	5	8	10
	ARI	0.97 0.74 0.71 0.66 0.60

Table 4 :

 4 Real data set (radians).Patient 1 st angle 2 nd angle 3 rd angle 4 th angle 5

	th angle

Supplementary Material

Supplementary Material containing the Appendix is attached to this submission.


Anh Khoa Nguyen
email: anh-khoa.nguyen@univ-orleans.fr

Eric Blond
email: eric.blond@univ-orleans.fr

Thomas Sayet
email: thomas.sayet@univ-orleans.fr

Athanasios Batakis
email: athanasios.batakis@univ-orleans.fr

Emmanuel De Bilbao
email: emmanuel.debilbao@univ-orleans.fr

Minh Duc Duong

NUMERICAL SIMULATION OF IMPREGNATION IN POROUS MEDIA BY SELF-ORGANIZED GRADIENT PERCOLATION METHOD

Keywords: impregnation, porous media, capillary pressure profile, gradient percolation

d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

INTRODUCTION

The numerical modelling of the impregnation process requires a multiphysics model taking into account the material properties. Yet, it often demands large computing facilities.

The goal of this paper is to propose a new approach, which doesn't use classical modelling by partial differential equations and the associated numerical methods, to predict the capillary pressure profiles without spurious oscillations [START_REF] Verruijt | An accuracy condition for consolidation by finite elements[END_REF] and with reduced computational cost. To develop and present the basics of this method, the simplest case of non-reactive impregnation, for quasi one-dimensional problem, is developed.

A new numerical algorithm based on the gradient percolation theory is proposed. The initialization of the algorithm is driven by an analytic solution of the homogeneous diffusion equation, which is a convolution between a Probability Density Function (PDF) and a smoothing function [START_REF] Jähne | Handbook of Computer Vision and Applications[END_REF]. The evolution of the capillary pressure profiles with time is reproduced by the self-evolution of the standard deviation of the PDF. This model is therefore named Selforganized Gradient Percolation model (SGP). In order to test this model, its solutions (i.e., the capillary pressure profiles and the mass gain curve) are compared with those obtained by F.E.M and by experiments, respectively.

SELF-ORGANIZED GRADIENT PERCOLATION (SGP) MODEL

Self-organized Gradient Percolation (SGP) model

The Self-organized Gradient Percolation (SGP) model defines that the porous medium is considered as a random network model [START_REF] Grimmett | Percolation[END_REF][START_REF] Nolin | Critical exponents of planar gradient percolation[END_REF] where each site has a local state ( ) which should be transformed into the local average saturation and has Gaussian distribution [START_REF] Nguyen | Simulations of the growth of cities[END_REF] with mean ( ) and variance as:

( )~ ( ( ), ) (1)
The Eq. [START_REF] Verruijt | An accuracy condition for consolidation by finite elements[END_REF] shows that the local average saturation is driven by the PDF of Gaussian distribution that is used to construct the self-organization of the local profile such that:

(i) its initialization is proved as an analytic solution of the homogeneous diffusion equation [START_REF] Jähne | Handbook of Computer Vision and Applications[END_REF]; (ii) for further time steps, the evolution of the capillary pressure profile should be reproduced by the self-evolution of the PDF.

Initialization of the SGP model

The SGP model aims at predicting the capillary pressure profile at any time that is inferred from a capillary pressure profile at initial time-step (i.e., initialization of the SGP model). As mathematically proven [START_REF] Jähne | Handbook of Computer Vision and Applications[END_REF], the initialization is the analytic solution of the homogeneous diffusion equation (i.e. Richard's equation [START_REF] Celia | A general mass-conservative numerical solution for the unsaturated flows equation[END_REF]) as the following formulation:

, = , * (2)
where 〈 * 〉 is convolution operator; stands for a smoothing point-spread function; , is the function of the local average saturation at initial time-step , and is described as the PDF of a distribution with standard deviation , maximum saturation and residual saturation as:

, = + ( -) - 2 (3)
where characterizes the type of a distribution: = 1 or = 2 indicate that , is the PDF for a Normal distribution or for a Laplace distribution by means of statistics, respectively. The Eq. ( 3) might be presumed to be more simple and convenient in application by the Figure 1. To determine the capillary pressure profiles, it is needed to identify the evolution of the local average saturation. It is thus proposed to reproduce the evolution of the local average saturation by that of the standard deviation of the PDF.

Time evolution

Let us remark that at each time-step , a standard deviation function is denoted by to be different from constant standard deviation ; the evolution of the abscissa with time-step is denoted by , . The evolution of standard deviation is thus defined by:

= + , (4)
The evolution of the standard deviation function at the second time-step, i.e. = 1, is expressed by the following form:

= + - = + ( -) (5)
where is a function of the capillary pressure [START_REF] Christiansen | Evaluation of methods for measuring relative permeability of Anhydrite from the Salado Formation: Sensitivity Analysis and Data Reduction[END_REF]. The general standard deviation function, i.e. = , is therefore deduced as follows:

= + - = + + , - (6)
At the local scale, the pore space can be assumed as a capillary (vertical for 1D model proposed herein). Thus, the evolution of the abscissa with time-step , i.e. , , is driven by Poiseuille's equation [START_REF] Extrand | Forces, pressures and energies associated with liquid rising in nonuniform capillary tubes[END_REF] as follows:

, = = • = • + = • ∆ (7)
where depends on capillary diameter and liquid properties, and are the total and the hydrostatic pressure at time-step , respectively; ∆ is the difference between capillary and hydrostatic pressures. In the case of the vertical capillary tube, the Eq. ( 7) points out that the transient simulation stops at the steady-state when there is no longer difference between capillary pressure and hydrostatic pressure (i.e., ∆ = 0).

Continuity and Boundary Conditions

As shown in the Eq. ( 2), a convolution procedure is first employed to introduce the continuity of the SGP model. Indeed, the procedure is related to the local average saturation and the smoothing point-spread function as the following formulation:

( ) = ( * )( ) = ( -′) ( ′) ∈ (8)
where ( ′) = ∑ ( ′) , Ω = { ′ ∈ ℤ : | ′| ≤ 1}, ( ′) is the Kronecker function, and is number of site ∈ Ω.

For impregnation in porous media, several boundary conditions related to the physical exchange made on the concerned interface are distinguished by: the case where liquid can/ cannot flow out of the boundary (i.e., the constant air pressure and the draining condition/ the undrained condition, respectively). The choice of the boundary conditions depends on the interpretation of the physical phenomena. In the SGP model, the convolution of the values on the boundary surface with the smoothing function is then employed to model the boundary conditions.

THE VALIDATION OF THE SGP METHOD

One-dimensional capillary rising tests were performed on porous samples [START_REF] De Bilbao | Identification of Transport Properties of Refractories: Intrinsic permeability and Capillary Pressure Curve[END_REF]. The amount of impregnated liquid into the porous medium was measured by weighing the loss of liquid in the bath as can be seen in Figure 2. Porous material and liquid for the tests are reported in Table 1. To validate the SGP method, a Finite Element (F.E.) model has been implemented in the software "ASTER" [START_REF] De Bilbao | Identification of Transport Properties of Refractories: Intrinsic permeability and Capillary Pressure Curve[END_REF]. To fit the capillary pressure curve, two well-known models (c.f. Brook's model and van Genuchten's model) [START_REF] Brooks | Hydraulic Properties of Porous Media[END_REF][START_REF] Van Genuchten | A closed-form for predicting the hydraulic conductivity of unsaturated soils[END_REF] are employed where pressure reference or and empirical parameter or need to be determined. To apply the algorithm of the SGP model on both tests, the values of input data are reported in Table 3.

Comparisons between numerical results and experimental data

For the tests 1 and 2, because there is no experimental data for the capillary pressure profiles, the mass gain curve obtained through the SGP model is compared with that obtained experimentally [START_REF] De Bilbao | Identification of Transport Properties of Refractories: Intrinsic permeability and Capillary Pressure Curve[END_REF] (Figure 3). The Figure 3 shows that numerical results from the SGP algorithm fit well with experimental data.

Comparisons between F.E and SGP results

For the tests 1 and 2, the comparisons of the capillary pressure profiles between the SGP model and the FEM model are shown in Figure 4. The Figure 4 shows that F.E. results fit well with SGP results too. Concerning the comparisons between F.E. and SGP results for the test 2, the CPC from Brook's model is employed, which is assumed to be a PDF of the Laplace distribution, i.e. = 1 for the SGP model. Thus, the capillary pressure profiles obtained by the SGP model always remain the form of this PDF. That's why the F.E. and SGP solutions are not exactly the same (as can be seen in the Figure 4 (b)).

The computational costs for the two tests are also reported in Table 4. By comparisons of the CPU time, the Table 4 shows that the computational time of the SGP model is lower than the one of F.E.M. Finally, the objectives of this study are entirely satisfied: (i) the solutions of the SGP method are free from spurious oscillations and (ii) the computational time is reduced.

CONCLUSIONS

In this study, a novel algorithm based on the gradient percolation theory, named Selforganized Gradient Percolation (SGP) algorithm, is proposed. Here, the gradient percolation theory is used to compute the local saturation. The link between the mathematics and the physical issue (i.e., the non-reactive impregnation phenomenon) is done thanks to the initialization of the algorithm by an analytic solution of the Richard's equation. Furthermore, the convolution operator allowed to ensure the spatial continuity of the wetting fluid and to check the boundary conditions of the problem. The comparisons with experimental measurements show a good agreement for the quasi one-dimensional case of vertical impregnation. The comparisons with F.E.M exhibits a lower computational cost for the SGP method. First results are promising with respect to possible generalization to 3D case.

Figure 1 :

 1 Figure 1: A simplified representation of the Eq. (3)

Figure 2 :

 2 Figure 2: The non-reactive impregnation test

Figure 3 :

 3 Figure 3: Measured and simulated mass gain curve: (a) for the test 1 and (b) for the test 2.

Figure 4 :

 4 The evolution of the simulated capillary pressure profiles: (a) for the test 1 and (b) for the test 2.

Table 1 :

 1 Materials and liquid for the tests

	Test	Porous sample	Liquid
	1	Alumina 99%	Glycerine
	2	Alumina 99%	Oil

Table 2 :

 2 Values of physical chracteristics used of the porous materials and the liquid

			Values
	Properties			Units
		Test 1	Test 2
	Initial porosity	0.2	0.19
	Viscosity of the fluid,	1.02	0.35	•
	Mass density of fluid,	1260	892	•
	Intrinsic permeability,	9.5 × 10	7.44 × 10

Table 3 :

 3 Values of input data used for the SGP algorithm

	Test	( ⁄ )	or	or	(Pa)
	1	1,16 • 10	0,55	1100
	2	1,9 • 10	2,62	4870

Table 4 :

 4 Computational time of the SGP model and the F.E model.

	Test	CPU time (seconds)
	SGP model	F.E model
	6.52		25.5
	1.1		12.1