Pierre Gloaguen
email: pierre.gloaguen@agroparistech.fr

Laetitia Chapel

Chloé Friguet

Romain Tavenard

Scalable clustering of segmented trajectories within a continuous time framework Application to maritime traffic data

Keywords: AIS data, Continuous time, Hierarchical Dirichlet process, Scalability, Stochastic variational inference, Trajectory clustering

In the context of the surveillance of the maritime traffic, a major challenge is the automatic identification of traffic flows from a set of observed trajectories, in order to derive good management measures or to detect abnormal or illegal behaviours for example. In this paper, we propose a new modelling framework to cluster sequences of a large amount of trajectories recorded at potentially irregular frequencies. The model is specified within a continuous time framework, being robust to irregular sampling in records and accounting for possible heterogeneous movement patterns within a single trajectory. It partitions a trajectory into sub-trajectories, or movement modes, allowing a clustering of both individuals' movement patterns and trajectories. The clustering is performed using non parametric Bayesian methods, namely the hierarchical Dirichlet process, and considers a stochastic variational inference to estimate the model's parameters, hence providing a scalable method in an easy-todistribute framework. Performance is assessed on both simulated data and on our motivational large trajectory dataset from the Automatic Identification System (AIS), used to monitor the world maritime traffic: the clusters represent significant, atomic motion-patterns, making the model informative for stakeholders.

Introduction

For the last thirty years, the tracking of movements of individuals or objects has been eased by the increasing development of several tracking devices, such as the Global Positioning System (GPS) that allows recording geographical positions through time. Such devices have led to large movement databases that store trajectory data, in which each point represents a position in space at a given time. These databases contain a great deal of knowledge and require analysis [START_REF] Demšar | Analysis and visualisation of movement: an interdisciplinary review[END_REF] such as the extraction of movement patterns across diverse trajectories in order to perform human mobility or traffic monitoring [START_REF] Li | Traffic density-based discovery of hot routes in road networks[END_REF][START_REF] Zheng | Understanding mobility based on gps data[END_REF], motion prediction [START_REF] Sung | Trajectory clustering for motion prediction[END_REF], human action recognition [START_REF] Yao | Learning universal multiview dictionary for human action recognition[END_REF] or detection of abnormal behavior in maritime routes (Vespe and Mazzarella, 2016).

A striking example of large movement database comes from the worldwide maritime traffic surveillance. This monitoring relies on several sources of data, in a rising context of maritime big data [START_REF] Garnier | Exploiting the potential of the future "maritime big data[END_REF]. Among these sources lies the Automatic Identification System (AIS), which automatically collects messages from vessels around the world, at a high frequency. AIS data basically consist in GPS-like data, together with the instantaneous speed and heading, and some vessel specific static information [START_REF] Clazzer | On the characterization of ais traffic at the satellite[END_REF]. An example of such data is presented in Figure 1 (left), considering 6 months of AIS data of vessels steaming in the Ushant traffic separation scheme in Brittany, west of France [START_REF] Fablet | Next Step for Big Data Infrastructure and Analytics for the Surveillance of the Maritime Traffic from AIS & Sentinel Satellite Data Streams[END_REF]. These data are characterized by their diversity as they ( 1) are collected at different frequencies (2) have different lengths (3) are not necessarily regularly sampled (4) represent very different behaviors, but (5) share common trends or similar subparts (called hereafter movement modes).

Fig. 1 (Left) AIS data in the Ushant traffic separation scheme (Brittany, France), gathered during 6 months (7M. GPS obs.). Source: CLS, Brest, France and ANR SESAME [START_REF] Fablet | Next Step for Big Data Infrastructure and Analytics for the Surveillance of the Maritime Traffic from AIS & Sentinel Satellite Data Streams[END_REF] (Right) Three trajectories extracted from the dataset, sharing common movement modes (characterized by the directions West-South-West and South-South-West). The green (+ shaped points) trajectory, however, subsequently adopts another movement mode (South-South-East). The red (triangles) and blue (crosses) trajectories belong to the same cluster, which is different from the one of the green (+ shaped) trajectory.

One major challenge in this context is the extraction of movement patterns emerging from the observed data, considering trajectories that share similar movement modes (see Figure 1 (right)). This issue can be restated from a machine learning point of view as a large-scale clustering task involving the definition of clustering methods that can handle such complex data while being efficient on large databases, and that are able to both cluster trajectories as a whole and detect common sub-trajectories.

Regarding the problem of clustering entire trajectories, existing methods in the literature follow two main trends:

-The first class of methods is a density-based clustering framework that relies on the definition of typical motions in the trajectories [START_REF] Ester | A density-based algorithm for discovering clusters in large spatial databases with noise[END_REF][START_REF] Rinzivillo | Visually driven analysis of movement data by progressive clustering[END_REF][START_REF] Yang | Tad: A trajectory clustering algorithm based on spatial-temporal density analysis[END_REF]. Similarities between full trajectories are considered to take into account local details in the segments. -An other category of methods relies on distance computation between trajectories, clustering together trajectories that are similar. The most classical distances used in this context are Dynamic Time Warping [START_REF] Sakoe | Dynamic programming algorithm optimization for spoken word recognition[END_REF] or LCSS [START_REF] Vlachos | Discovering similar multidimensional trajectories[END_REF], among others. -The last class of methods is a rule-based clustering framework, that relies on preconceived decision rules to build each cluster (Vespe and Mazzarella, 2016). The main advantage of these methods is their simplicity, but they suppose the prior knowledge of decision rules, and afterwards make strong independence assumption on observed points, that might not hold in practice.

In a more parametric approach, one can aim at modelling the heterogeneous nature of movement data [START_REF] Nathan | An emerging movement ecology paradigm[END_REF] to define a proper normality model, whose structure is made of clusters.

In this context, trajectories would belong to a cluster, and, moreover, a trajectory could be the bringing together of different moving patterns or movement modes.

In this context, different trajectories may share similar movement modes in a given portion. To learn common sub-trajectories in the clustering process, [START_REF] Lee | Trajectory clustering: a partition-and-group framework[END_REF] propose a partition-and-group framework based on a dedicated distance. Wang et al. (2011b) propose the definition of a normality model for trajectories through an unsupervised clustering framework inspired by topic models [START_REF] Steyvers | Probabilistic topic models[END_REF]. Initially developed in the field of text analysis, topic models aim at classifying documents depending on their main topics. The clustering approach then estimates the different topics in each document as well as clusters documents together regarding their topics. In the transposition of topic models to trajectory analysis, Wang et al. (2011b) consider trajectories as documents and quantized GPS observations (position and velocity) as words. It is assumed that each point of the trajectory belongs to a semantic region (a topic), called hereafter a movement mode. A movement mode is therefore a specific distribution to be estimated from trajectory data. Trajectories are then considered as a mixture of movement modes, and a trajectory cluster is a set of trajectories with a common mixture distribution. Estimation of both movement modes and trajectory clusters is performed on a discretized space within a Bayesian framework using Gibbs sampling. If this approach still assumes independence of GPS observations (conditionally to their semantic region), it takes into account different heterogeneity levels for movement data. However, it requires a quantization of the space whose influence over results is not investigated, and relies on Gibbs sampling for inference, which is known to be hardly scalable [START_REF] Blei | Variational inference: A review for statisticians[END_REF].

Nevertheless, none of the aforementionned methods are able to deal with large databases, neither consider the continuous nature of the data. For an in-depth review of spatial trajectory clustering algorithms, one can refer to [START_REF] Yuan | A review of moving object trajectory clustering algorithms[END_REF] or (Zheng, 2015, Section 6).

Given the challenges mentionned above, we propose in this paper a novel model, based on topic models, to cluster large and complex trajectory datasets in a scalable framework.

Following these works, the approach proposed here distinguishes yet from existing works in the definition of GPS observations distributions as it is now assumed that (1) the distribution lives in a continuous space, avoiding the quantization of the (position/velocity) space, and (2) within a movement mode (or a topic), observations are not independent. Each movement mode is modeled through the bivariate velocity as a continuous time Gaussian process, namely the Ornstein Uhlenbeck Process (OUP) [START_REF] Uhlenbeck | On the theory of the brownian motion[END_REF], that has been proposed in movement ecology as a flexible framework to model the fundamental units of movement [START_REF] Gurarie | Correlated velocity models as a fundamental unit of animal movement: synthesis and applications[END_REF]. Thanks to its continuous time formulation, the OUP allows handling irregularly sampled data. Moreover, as a Gaussian process, it accounts for observation autocorrelation. The unsupervised clustering is made in a non parametric Bayesian context using a Hierarchical Dirichlet Process framework [START_REF] Teh | Hierarchical dirichlet processes[END_REF]. To perform movement mode estimation in a scalable approach, we use asymptotic properties of the OUP together with stochastic variational inference [START_REF] Hoffman | Stochastic variational inference[END_REF]. Contrary to MCMC methods that obtain samples from the posterior distribution, the variational approach solves an optimization problem, allowing the use of distributed computation and stochastic optimization to scale and speed up inference [START_REF] Blei | Variational inference: A review for statisticians[END_REF].

The remainder of this paper is organized as follows: in Section 2, the hierarchical parametric framework that models trajectory data is fully specified, and the proposed scalable approach to estimate model parameters from data is described in Section 3. The method is evaluated in Section 4, both on simulated data sets, to assess inference performance, and on a new dataset of current AIS data that is released with the paper, to show the interest of our scalable approach in a realistic context.

Movement model

In this section, we define a parametric framework to model trajectory data, i.e. sequences of geographical positions recorded through time.

In the following, the fundamental modelling of movement is first described, considering the individual's velocity process. Velocity is observed directly (as in the AIS context) or is at least computed from the successive geographical locations of the individuals. The positions and velocities of an individual are modeled as two continuous time processes in R 2 , denoted by (X t ) t≥0 and (V t ) t≥0 . Following a common paradigm, we assume that a moving individual's trajectory might be the collection of heterogeneous patterns [START_REF] Nathan | An emerging movement ecology paradigm[END_REF], namely the movement modes. Different movement modes along a trajectory refer to different ways of moving in terms of velocity distribution, reflecting different behaviors, activities, or routes. It is assumed that the set of possible movement modes of tracked individuals during their trajectories is countable, and that a given movement mode can be adopted by several individuals.

Then, the modelling framework aims to account for two levels of heterogeneity possibly present in trajectory data: (1) heterogeneity of an individual's movement within a single trajectory, and (2) heterogeneity between observed trajectories of several individuals.

Fundamental unit of movement Similarly to [START_REF] Gurarie | Correlated velocity models as a fundamental unit of animal movement: synthesis and applications[END_REF], a movement mode is assumed to be characterized by a specific correlated velocity model, defined in a continuous-time framework. Formally, during a time segment [τ 1 ; τ 2 ], if an individual adopts the movement mode k, then its velocity process (V t ) τ 1 ≤t≤τ 2 is assumed to be the solution of the following Stochastic Differential Equation (SDE) (see Øksendal (2003) for a formal definition of SDEs):

dV t = -Γ k (V t -µ k ) dt + Σ k dW t , τ 1 ≤ t ≤ τ 2 , V τ 1 = v τ 1 (1)
where:

µ k ∈ R 2 is the asymptotic mean velocity of the k-th movement mode; -Γ k is a 2 × 2 real-valued matrix, and is an autocorrelation parameter, modelling the recall force of the process towards the mean µ k ; -Σ k is a 2 × 2 real-valued matrix, and is a diffusion term, modelling the variability of the process around the mean µ k ; -W t is a bivariate standard Brownian motion, modelling the stochasticity of the individual's velocity process; v τ 1 is the initial condition of the SDE: the individual's velocity at time τ 1 .

The solution to Eq. ( 1) is a well known continuous time stochastic process, the Ornstein Uhlenbeck Process (OUP) [START_REF] Uhlenbeck | On the theory of the brownian motion[END_REF], also known as the mean reverting process. This mean reverting property is controlled by the recall parameter Γ k and makes the OUP suitable to describe movement modes of individuals, that are often characterized by a mean velocity regime, which is reached by rapid and brutal accelerations (or decelerations, see Figure 2). The OUP satisfies the three following properties:

1. (V t ) τ 1 ≤t≤τ 2 is a continuous time Markov process; 2. (V t ) τ 1 ≤t≤τ 2 is a Gaussian process, i.e. the law of the random variable

V t |V τ 1 = v τ 1
is Gaussian, with explicit mean and covariance (given in Appendix A); 3. if both eigenvalues of Γ k are positive, then (V t ) t≥τ 1 is an asymptotically stationary Gaussian Process, i.e., the law of V t |V τ 1 = v τ 1 satisfies:

V t |V τ 1 t→∞ -→ N (µ k ,Λ k ) , (2)
where Λ k is the asymptotic covariance matrix of the OUP (given in Appendix A). As discussed in [START_REF] Gurarie | Correlated velocity models as a fundamental unit of animal movement: synthesis and applications[END_REF], the OUP offers a flexible framework to model wide range of autocorrelated velocity processes, such as highly directed movements or rotational ones. Knowing the individual's velocity process, its position process, starting at position x τ 1 , is obtained by deterministic integration over time:

X t = X τ 1 + t τ 1 V s ds.
The resulting process (X t ) τ 1 ≤t≤τ 2 is known as an Integrated Ornstein Ulhenbeck Process (IOUP), which remains a Gaussian Process.

Within trajectories heterogeneity To account for heterogeneous ways of moving during a single trajectory, we assume that an individual's trajectory in a time segment [0, T ] is a sequence of IOUPs. Formally, there exists a sequence of times τ 0 = 0 < τ 1 < • • • < τ L = T such that the trajectory is a sequence of L successive IOUPs, each IOUP being defined over a time segment [τ l , τ l+1 ]. A simulated example of such a sequence is shown on Figure 3.

Therefore, a trajectory is characterized by (1) the movement modes adopted within the trajectory, (2) the time spent by the tracked individual within each movement mode and (3) the transitions from one movement mode to another. Here, we do not impose specific modelling for the last two points. Our first main goal is then to estimate the different movement modes within trajectories from a given dataset.

Between trajectories heterogeneity In addition to this segmentation problem, we assume that the dataset is composed of different (unknown) clusters of (entire) trajectories. A cluster is characterized by similar distributions over movement modes. In other words, two trajectories belong to the same cluster if they are composed of the same movement modes in similar proportions. This clustering problem is illustrated in Figure 4. Our second main goal is then to cluster together such trajectories.

To summarize, the overall modelling framework for the two level clustering is depicted on Figure 5.

Scalable two step clustering

The following three problems must be solved:

1. Characterizing different movement modes present in the dataset; 2. For each GPS observation, estimating in which movement modes it belongs; 3. Clustering together trajectories that have the same distribution over movement modes.

In the previous section, we defined a realistic framework to model movement data, introducing a continuous time autocorrelated model over velocity process. The resulting framework therefore consists in a continuous time regime switching diffusion model. Estimating both parameters and movement mode assignation of data is in this case a statistical challenge, even in the case of a single trajectory [START_REF] Patterson | Statistical modelling of individual animal movement: an overview of key methods and a discussion of practical challenges[END_REF]. For instance, [START_REF] Blackwell | Exact bayesian inference for animal movement in continuous time[END_REF] use thinning Poisson process and

F j Z j 1 Z j 2 . . . Z j L j V j 0 V j 1 . . . V j i V j i+1
. . .

Hidden movement modes dynamics

Observed Ornstein Uhlenbeck dynamics Hidden trajectory cluster j = 1, . . . , J MCMC methods to perform parameter estimation from movement data, resulting in highly time demanding approach. A big challenge for estimation methods therefore remains their scalability, especially in the dual clustering context we target here.

In order to perform scalable parameter inference and clustering of both trajectories and GPS observations (into movement modes), we adopt a pragmatic two-step approach that takes advantage of the inherent properties of the OUP and are described with more details hereafter:

Step-1 A first dual clustering is performed based on a simpler independent Gaussian mixture model, in order to estimate potential movement modes and trajectory clusters: it allows getting rid of within mode autocorrelation in the inference, and therefore eases the computations. The Gaussian hypothesis in this case is rather natural, as the OUP stationary distribution is Gaussian.

Step-2 Among the estimated movement modes, only those meeting a temporal consistency constraint are kept. Parameters of these consistent movement modes are then estimated, and used to reassign observations that were assigned to inconsistent movement modes. It ensures that only trajectory segments for which this stationary distribution was reached are kept to estimate movement modes.

Finally, we discuss the algorithmic complexity of the proposed two-steps procedure, justifying that the procedure remains tractable for large datasets while using a continuous time framework ine the movement modes modelling.

3.1

Step 1: simplified model based on an independence assumption

Hierarchical model

The data set is composed of J independent trajectories that we aim to cluster into C groups. This number of clusters is unknown and has to be inferred from data. Each trajectory j is a sequence of n j velocities (V j i ) i=1,...,n j . For a trajectory j, we denote F j = c if it belongs to the cluster c (1 ≤ c ≤ C). We suppose that:

P(F j = c) = w c
where w c is an unknown parameter to be estimated and we denote w = (w 1 , . . . , w C ).

In this first step, it is assummed that each cluster c is a Gaussian mixture distribution having K components and mixture weights π c = (π c 1 , . . . , π c K ). Again, this number K is unknown and has to be inferred from data. Each mixture component models a movement mode. For each V j i , we denote Z j i = k if the velocity belongs to the kth movement mode. Conditionnally to the cluster assignation of the trajectory j, we therefore have:

P Z j i = k|F j = c = π c k
We denote π = π 1 , . . . , π C the set of all unknown movement modes weights to be estimated.

In this simplified model, conditionnally to their movement mode, all the velocities are mutually independent (as assumed in Wang et al., 2011b). This strong hypothesis eases the computation and ensures the scalability of the method as a first step, but will be relaxed in the second step. We therefore suppose that:

V j i | Z j i = k ∼ N µ k ,Λ -1 k , 1 ≤ k ≤ K where µ k ∈ R 2 (resp. Λ k ∈ M 2×2
) is the mean velocity (resp. the precision matrix) of the k-th movement mode. The set of unknown movement parameters

{(µ k ,Λ k )} k=1•••K is denoted T .
The movement parameter set T , the trajectory clusters weights w and, for each cluster c, the movement mode weights π c are the global variables (reusing the terminology used in Hoffman et al. ( 2013)) of the model that have to be estimated. Moreover, for each trajectory j, the cluster label F j has to be estimated. Finally, for each observation V j i , the movement mode label Z j i has to be estimated. Sets F = F j 1≤ j≤J and Z = Z j i 1≤ j≤J, 1≤i≤n j are the local variables of the model to be estimated. The resulting hierarchical structure of this model is depicted on Figure 6. The next section describes the scalable inference approach to estimate both local and global variables from the data.

Bayesian estimation of the parameters using stochastic variational inference

In the first step, the inference of the local and global variables is made within a Bayesian context, considering unknown parameters as random variables. The infer-

w π T Global variables F j Z j i
Local variables

V j i

Cluster level

Observations

Trajectory level j = 1, . . . , J i = 1, . . . , n j ence therefore aims at obtaining the joint posterior distribution of local and global variables, denoted by p(T , w, π, F, Z|V ).

Together with the prior specification, a classical problem in Bayesian clustering is to set the number of clusters. In the following, we adopt a Bayesian non parametric approach by considering a possibly infinite set of movement modes and clusters in the data.

First, for a given movement k, we consider as prior distribution for (µ k ,Λ k ) a Gaussian-Whishart having hyperparameters m 0 , ρ 0 , γ 0 and W 0 . This distribution is a common prior for Gaussian parameters as it provides nice conjugacy properties that we'll use in our scalable purpose. The hyperparameters m 0 and W 0 are location hyperparameters for the mean and the precision matrix, while ρ 0 and γ 0 are concentration parameters, setting the amount of information carried by the prior.

As we do not specify a priori the number of possible movement modes, we assume that this number is infinite.

Each trajectory cluster c is therefore an infinite mixture distribution whose components are Gaussian distributions. The mixture weights prior is given by a stick breaking distribution, denoted by GEM(1, β ), where β is a concentration hyperparameter (see [START_REF] Sethuraman | A constructive definition of dirichlet priors[END_REF] as well as Appendix B.2 for details).

Finally, the amount of such clusters is itself infinite, and the prior distribution over cluster weights is again a stick breaking distribution GEM(1, α) where α is a concentration hyperparameter. The full prior specification is then given by the following equations:

               µ k ,Λ -1 k i.i.d. ∼ GW (m 0 , ρ 0 , γ 0 ,W 0 ) k ≥ 1 w ∼ GEM(α) π c i.i.d. ∼ GEM(β ) c ≥ 1 F j |w i.i.d. ∼ w 1 ≤ j ≤ J Z j i |π, F j i.i.d. ∼ π F j 1 ≤ j ≤ J, 1 ≤ i ≤ n j
The model depicted in Figure 6 together with these priors is known as a Hierarchical Dirichlet Process [START_REF] Teh | Hierarchical dirichlet processes[END_REF]. In this Bayesian context, the expression of the posterior density p(T , w, π, F, Z|V ) has no known analytical form. A classical way to obtain Bayesian estimators of quantities of interest would be alternatively to obtain samples from this posterior distribution. This can be done using MCMC methods, that can be implemented for Dirichlet processes using a Gibbs sampler [START_REF] Neal | Markov chain sampling methods for dirichlet process mixture models[END_REF], and were used in Wang et al. (2011b) for the inference in their discrete-space framework. However, it is known that the Gibbs sampler does not scale properly [START_REF] Blei | Variational inference: A review for statisticians[END_REF], and therefore could not be used in the context of this paper.

A more scalable alternative to this simulation-based approach is the use of variational inference (see [START_REF] Blei | Variational inference: A review for statisticians[END_REF] for a recent review). Variational methods reduce the inference to an optimization problem by minimizing a divergence (typically the Kullback-Leibler divergence [START_REF] Kullback | On information and sufficiency[END_REF]) between the target posterior distribution and the members of a simpler family of distributions, the variational family. An appropriate and common variational family Q is the one which satisfies the mean field assumption, i.e. the set of q distributions that can be fully factorized. Our variational family is therefore such that:

q(T , w, π, F, Z) = q(w) ∏ c q(π c ) ∏ k q(µ k ,Λ k ) ∏ j q(F j ) ∏ i, j q(Z j i ). (3)
The mean field inference problem reduces to find q * such that:

q * = argmax q∈Q E q [log p(T , w, π, F, Z,V )] -E q [log q(T , w, π, F, Z)] , (4)
where E q [•] denotes the expectation with respect to the p.d.f. q. The right hand side of ( 4) is known as the evidence lower bound (ELBO), and can be computed for appropriate families of distributions Q, and thus can be maximized.

In addition to this mean field property, the variational approximation of the posterior distribution is restricted to finite sets of parameters. Formally, the posterior distributions of the infinite cluster (resp. movement mode) weights are approximated by a distribution on a finite set of weights having C max (resp. K max ) elements. Here, C max and K max are variational parameters given by the user. This variational approximation is known as the truncated stick breaking distribution [START_REF] Hoffman | Stochastic variational inference[END_REF] (see Appendix B for details).

A known algorithm to compute q * is the coordinate ascent variational inference (CAVI, [START_REF] Bishop | Pattern Recognition and Machine Learning[END_REF]). This iterative algorithm starts from an initial guess q (0) for the optimal variational distribution, and successively updates each of its components by supposing the others known, and computes an expectation with respect to their distribution (Algorithm 1). For the model presented in Section 3.1 and the chosen prior distributions, all needed expectations can be computed explicitly (see Appendix B for details).

As described in [START_REF] Hoffman | Stochastic variational inference[END_REF], this algorithm can be seen as a gradient ascent, and is therefore suitable for stochastic approximations, resulting in the stochastic variational inference (SVI).

The estimation algorithm therefore reduces to:

1. sample uniformly a batch from the data set, 2. compute the expectations (as given in Appendix B) using only this batch, 3. update variational distributions using these expectations, as described in [START_REF] Hoffman | Stochastic variational inference[END_REF].

Algorithm 1 Coordinate ascent variational inference algorithm 1: Denote p = p(V, w, π, θ , F, Z) 2: Set q (0) = q (0)

w (•)q (0) π (•)q (0) T (•)q (0) F (•)q (0) Z 3: while convergence is not reached do 4: Local variables update 5: log q (i+1) Z (•) = E q (i) w q (i) π q (i) T q (i) F [log p] 6: log q (i+1) F (•) = E q (i) w q (i) π q (i) T q (i+1) Z [log p] 7: Global variables update 8: log q (i+1) w (•) = E q (i) π q (i) T q (i+1) F q (i+1) Z [log p] 9: log q (i+1) π (•) = E q (i+1) w q (i) T q (i+1) F q (i+1) Z [log p] 10: log q (i+1) T (•) = E q (i+1) w q (i+1) π q (i+1) F q (i+1) Z [log p] 11: At last iteration M, set q * = q (M) w (•)q (M) π (•)q (M) T (•)q (M) F (•)q (M) Z
This procedure can be performed online, therefore with no need of storing the data (Wang et al., 2011a). Performances of these stochastic methods are widely discussed in Wang et al. (2011a) and [START_REF] Hoffman | Stochastic variational inference[END_REF].

It is worth noting here that all the needed computations can naturally be distributed (thanks to the independance simplification), as they are essentially a sum over simple operations involving single observations. Therefore, the SVI algorithm proposed here for the simplified model is widely scalable, unlike Gibbs sampling procedures.

As the optimization is done in a high dimensional space, and the algorithm only guarantees convergence to a local optimum, it is crucial to initialize the algorithm from different starting points, to ensure a good exploration of the space1 . Again, different runs of the algorithm can be distributed.

Step 2: Estimation of OUP parameters from clustering outputs

The previous section described the first step of our two-step approach to perform trajectory clustering and inference of movement modes. In this section, we define how the parameters of the movement modes, defined by Ornstein Uhlenbeck processes, are re-estimated from this first step. This first step gives as an output the set of optimal variational distributions, namely:

q * T , the posterior distribution of {(µ k ,Λ k )} 1≤k≤K max , the parameters (under the Gaussian mixture assumption) of the K max possible movement modes in the data set; q * w , the posterior distribution of the weights of the C max possible trajectory clusters in the data set; q * π = {q * π c } 1≤c≤C max the posterior distribution of the weights of the K max possible movement modes in each trajectory cluster c in the data set; for each mode-coherent sequence of consecutive observations Reallocate the sequence to the most likely movement mode

V j i 1 , • • • ,V j im do 5: if t im -t i 1 ≥ δ then 6: Add sequence V j i 1 , • • • ,V j im to set A[Z j i 1 ] 7: else 8: Add sequence V j i 1 , • • • ,V j im to
-q * F = q * F j 1≤ j≤J
for each trajectory j, the posterior probability of being in each possible trajectory cluster;

-q * Z = q * Z j i
1≤ j≤J,1≤i≤n j for each observation i of trajectory j, the posterior probability of being in each possible movement mode.

From the last two distributions, one can get estimators of clusters and movement modes present in the data. A classical estimator would be the maximum a posteriori (MAP), i.e. the cluster (resp. movement mode) label giving the maximum weight of the posterior multinomial distribution q * F j (resp. q * Z j i

).

In order to estimate the OUP parameters of the k-th movement mode, as defined in Section 2, a filter is applied on movement mode sequences based on their temporal consistency. Formally, for a trajectory j, let us consider a sequence of m successively recorded velocities V j i 1 , . . . ,V j i m , that were all estimated to belong to a same movement mode k. This sequence is said to be temporally consistent if its length t i m -t i 1 is larger than δ , a user chosen parameter, representing the minimal time lag for a movement mode. All consistent sequences estimated in a same movement mode are considered as independent realizations of OUP with common parameters. From the Markov and Gaussian properties of the OUP, the likelihood related to this data set can be easily maximized to obtain estimates of parameters µ k , Γ k and Σ k (see Appendix A). If an estimated movement mode k (i.e. a movement mode containing at least one observation) has no consistent sequence, this movement mode is considered as inconsistent. Finally, to refine the movement modes classification, observation sequences belonging to inconsistent movement modes are reallocated to the consistent movement mode whose parameters maximize their likelihood. The overall procedure is summarized in Algorithm 2.

Note that this second step depends on a parameter δ , which value is surely data dependent, but has a clear interpretation. Therefore, depending on the context, its value can easily be discussed with relevant experts or stakeholders.

The resulting consistent movement mode concept allows one to (1) have a good estimation of OUP parameters within a movement mode (as a consistent sequence will often be related to a large amount of points) and (2) filter out "noise" movement modes gathering few observations in a temporally inconsistent manner.

Algorithmic complexity

The goal of the two-step inference scheme is to make the approach tractable while still using fine-grained continuous time models inside movement modes.

Time complexity for our first step is at worse O(J ×t max × K max + J ×C max ) where t max is the maximal number of time steps for a trajectory in the data set. It should be pointed out that our approach relies on a variational approximation of the stick breaking distribution, that imposes to set two parameters K max and C max , as the maximal number of movement modes and trajectory clusters present in the data. These numbers should increase with n. It is worth noting that for our prior specification, a Dirichlet process of parameter α, the expected number of classes in a n data set is α log n (Ghosal and Van der Vaart, 2017, Proposition 4.8). The time complexity for this first step is hence O(β n log n + αJ log J) where n = ∑ j n j is the total number of points in the data set.

Then, estimating OU parameters for a given movement mode in the second step of our process is linear in the number of observations assigned to this movement mode. Since this has to be done for all movement modes, complexity for this second step is linear in the total number of observations in the dataset.

Overall, computational complexity of the inference step is then quasilinear in n and, as discussed above, parts of the computations involved can be distributed. We show in Section 4 that these two properties allow tackling large scale datasets of GPS data in reasonable time.

Experiments

Experiments on simulated data

To ensure that the inference proposed in Section 3 suits the model defined in Section 2, a numerical experiment is performed on simulated data.

Simulation set up.

A data set of 40 trajectories containing overall 8,000 observations is simulated, according to a model with 2 trajectory clusters, the two clusters being composed of respectively 2 and 3 movement modes.

Within each movement mode, velocities are drawn from an OUP whose parameters are movement mode dependent. All trajectories start at velocity (0, 0) and spend 50% of the whole time in the first movement mode before switching to a second movement mode. Trajectories of the first cluster spend then the remaining time in the second movement mode. The trajectories in the second cluster spend 12.5% in the second movement mode before switching to the third movement mode for the remaining time (37.5% of the overall time). Simulated data are shown in Figure 7.

Movement mode estimation: results after the first step. A first clustering is performed using the variational inference approach within a hierarchical Dirichlet process, as described in Section 3.1. The variational approximations of the stick breaking distributions are made with truncated breaking distribution with C max = 5 and K max = 12. Optimization is made independently from 200 starting points chosen randomly, performing 1 000 iterations of the SVI, with a decreasing gradient step along iterations. The best estimate is chosen as the final point maximizing the ELBO defined in Eq. ( 4). The estimated movement mode for a velocity V j i is computed as the maximum a posteriori of the weights distribution q * Z j i

. It results that among the K max = 12 possible movement modes, only 5 of them contain at least one observation. Overall, 88.3% of the points are attributed to the correct movement mode (Table 1). However, two extra movement modes are estimated, corresponding to the transition phases towards the light blue movement mode and leaving from it (see Figure 8).

It is worth noting here that this first phase is a continuous space version of the work of Wang et al. (2011b). This first phase of the presented approach, however, already improves on this baseline since it does not require to decide on a space quantization grid.

Movement mode estimation: results after the second step. Including the second step described in Section 3.2, movement modes for inconsistent sequences are re-estimated. Using δ -consistent sequences, OUP parameters are estimated and used to reallocate inconsistent ones. This results in a good reallocation of the problematic transitory phases having now an overall 96.1% good classification rate, as shown in Figure 8 and Table 1.

One can see here that coming back to the OU property enables correct classification of transitory phases of the movement. This is due to the recall force in the OU that enables attaching any observed velocity that is not close to the mean velocity in the movement mode as soon as there is a trend for velocities to move toward that asymptotic mean.

Trajectory clusters estimation.

As expected on such a simple example, the trajectory cluster assignation is 100% right. It is worth noting here that it would not be the case on trickier examples (not shown here), for instance, when the clusters are distinguished by the order of movement modes sequence, which cannot be captured by the HDP for Gaussian mixtures used here. Fig. 8 Estimated movement modes (in the bivariate velocity space) after one step (left) and two steps (middle). The ground truth is shown on the right. After one step, 5 movement modes are present in data, after step 2, only three are remaining. The related contingency table for good labelling is shown on Table 1 4.2 Experiments on maritime traffic data Dataset. We now validate our clustering approach by analysing real data. The considered dataset records 6 months of AIS data of vessels steaming in the area of the Ushant traffic separation scheme (in Brittany, West of France). This is an area with one of the highest maritime traffic density in the world, with a clear separation scheme of two navigation lanes. Different kinds of vessels are sailing in the area, from cargos and tankers with high velocity and straight routes to sailboats or fishing vessels with lower speed and different sailing directions. As such, the area is highly monitored to avoid collision or grounding, and a better analysis and understanding of the different ship behaviors is of prime importance.

The whole trajectory dataset is shown on Figure 1 and is available online2 . It consists in 18,603 trajectories, gathering at all more than 7 millions GPS observations. Only trajectories having more than 30 points are kept, time lag between two consecutive observations ranges between 5 seconds and 15 hours, with 95% of time lags below 3 minutes.

Two step movement modes estimation. Model fit is performed considering a maximum of 90 movement modes and 30 trajectory clusters, using non-informative priors. 200 runs of 50,000 iterations for SVI are performed independently, each run taking approximately 8 hours with SVI computations parallelized on a 8-core machine. The run leading to the highest ELBO is kept as the estimate.

The SVI leads to the estimation of 81 different movement modes containing (in the sense of maximum a posteriori probability) at least one observation, with 50 movement modes containing 95% of observations. The shapes (mean and covariance) of these movement modes are shown on Figure 9. One can note that a lot of movement modes lead to a same steaming direction, but at different speeds. In this context, choosing more informative prior or adding a penalty on movement modes covariance parameters in the ELBO optimization could lead to a model with less movement modes.

Trajectory clusters estimation. In addition to the movement modes estimated above, 29 (non-empty) trajectory clusters are identified. Three emblematic clusters are represented on Figure 10, plotting 50 randomly chosen trajectories with the same colour code as in Figure 9. The first two subfigures present the largest clusters (gathering respectively 18% and 12% of trajectories) and are composed of two preponderant movement modes: they represent the two main marine roads of the Ushant traffic separation scheme, the first one entering the English Channel, the second one exiting it. The third cluster gathers 5% of trajectories, and mostly puts together trips that are performed at low speed (below 5 knots), involving both sailing boats and fishing vessels.

Comparison with other clustering method. Our algorithm can be seen as an improved (continuous-space) version of the method presented in Wang et al. (2011b). As such, in order to validate the results of the clustering, we rather rely on a k-means algorithm using Dynamic Time Warping (DTW, [START_REF] Sakoe | Dynamic programming algorithm optimization for spoken word recognition[END_REF] as a distance between trajectories considering them as sequences of 2D velocity vectors (using DTW Barycenter Averaging method from [START_REF] Petitjean | A global averaging method for dynamic time warping[END_REF] for centroid estimation). It is worth mentionning here that this method has a complexity at worse O(JC max × t 2 max ) (where t max is the maximal number of points for a trajectory in the data set). Thus, running this algorithm on the full dataset was too costly and we subsampled one tenth of the original dataset to feed the clustering. We used the implementation from the python library tslearn [START_REF] Tavenard | tslearn: A machine learning toolkit dedicated to time-series data[END_REF]. We also run our algorithm on the same dataset and we obtain a Normalized Mutual Information between the clusterings equal to 0.63. The resulting five main clusters are presented in Figure 2 in the Appendix, together with a confusion matrix between those five clusters. As no ground truth is available, we cannot determine which solution fits better but we note by a visual inspection of the clusters that the DTW-based clustering fails at differenciating ships that follow the same traffic lane (i.e. identifying sub-routes in the main traffic lanes). One should also note that, by assigning movement modes to the observations, our methods provides extra explanation for cluster assignments, compared to k-means that solely relies on inertia minimization. For instance, we show on Figure 11 typical trajectories from two significant clusters (respectively 4th and 5th in terms of size). This example illustrates the relevance of our modelling framework in this context, where these two clusters differ in the initial route (one coming from southwest, the other from the south), thus differing in one movement mode, but then share two movement modes, for the rest of the route. Interest of the post-filtering step. The post filtering step was performed whith a δ parameter of 2 hours. In the present case, the main advantage of the second step, which consists in a re-estimate of inconsistent sequences, is to avoid the estimation of transitory movement modes. Figure 12 shows a trajectory belonging to the central cluster of Figure 10. The transitory phase of the vessel, corresponding to its change of direction, is first assigned a specific movement mode at step 1 (as would be done by the baseline from Wang et al. (2011b)). As this transition is brief and the corresponding observations can be considered consistent with the OUP attached with the blue movement mode, these observations are reassigned to that movement mode.

Conclusion

In this work, we have defined a generic framework for the clustering of large trajectory data sets.

In our modelling, a trajectory is described by both the successive positions and velocities. Each trajectory is seen as a succession of movement modes in which a Ornstein Uhlenbeck process is used to model the velocity. A cluster of trajectories is then characterized by a its distribution over movement modes.

This framework is then specified in continuous time and space, which makes its formulation insensitive to GPS sampling, relaxing assumptions of previous models of the literature that are based on space quantization.

Inference is done in a two-step scalable approach using stochastic variational inference for a conjugate hierarchical Dirichlet process. This framework is easy to distribute as demonstrated experimentally, and has a quasilinear complexity in the number of observations. Besides, we provide a dataset of several millions of observations in the AIS context to both validate our model and allow future competitive methods to compare on a real-world large-scale trajectory dataset. Note however that the proposed framework is generic and suits, to our opinion, to a wide range of trajectory data sets.

Future works include extension of our framework to model sojourn time in movement modes and/or transitions between modes. This would make anomaly detection possible, both at the observation level (abnormal observation given the estimated movement mode) and at the trajectory level (abnormal trajectory as an unlikely sequence of movement modes). Such an extended model could also be used as a fully generative model.

A Characteristics of the Ornstein Uhlenbeck process

Notations:

-I d is the d × d identity matrix; -For a matrix A, A T is the transposed matrix of A; -For two 2 × 2 matrices A and B, A ⊕ B denotes the Kronecker sum of A and B, defined by : A ⊕ B = A ⊗ I 2 + I 2 ⊗ B, where ⊗ denotes the Kronecker product. Note that in our case, A ⊕ B is a 4 × 4 matrix;

-For a matrix A = a 1 a 3 a 2 a 4 , vec (A) denotes the vectorization of A, i.e. the vector vec

(A) =     a 1 a 2 a 3 a 4     .
Transition density Suppose that the process (V t ) τ 1 ≤t≤τ 2 is solution of equation ( 1), then,

∀τ 1 ≤ t ≤ τ 2 , V t |V τ 1 = v τ 1 ∼ N (m vτ 1 k,∆ ,Λ k,∆ )
where ∆ = tτ 1 and

m vτ 1 k,∆ = e -Γ k ∆ v τ 1 + (I 2 -e -Γ k ∆ )µ k vec Λ k,∆ = (Γ k ⊕ Γ k ) -1 I 4 -e -(Γ k ⊕Γ k )∆ vec Σ k Σ T k Likelihood Let v be a sequence of observations v t 1 = v τ 1 , v t 2 , . . . , v tn = v τ 2 , discrete time observations
of a OUP starting at V τ 1 and such that V τ 1 is a random variable with p.d.f. (possibly depending on µ k ,Γ k and Σ k ) p τ 1 (•). Then the likelihood of the observed sequence is given by:

L (µ k ,Γ k , Σ k ; v) = p τ 1 (v τ 1 ) n-1 ∏ i=1 φ k,∆ i (v t i+1 ),
where ∆ i = t i+1t i and φ k,∆ i (v t i+1 ) is the p.d.f. of a Gaussian distribution with mean m vt i k,∆ i and covariance Λ k,∆ i Stationary distribution From above, one can see that if both e -Γ k ∆ and e -(Γ k ⊕Γ k )∆ vanishes to 0 when ∆ increases, then the process (V t ) t≥τ 1 is asymptotically stationary, and

V t |V τ 1 = v τ 1 distrib. -→ N (µ k ,Λ k ) where vec (Λ k ) = (Γ k ⊕ Γ k ) -1 vec Σ k Σ T k .
The vanishing condition is satisfied when both eigenvalues of Γ k are positives.

B Computation for variational inference

B.1 Prior specification

Overall, the variables prior distributions are defined as follows:

               µ k ,Λ -1 k i.i.d. ∼ GW (m 0 , ρ 0 , γ 0 ,W 0 ) k ≥ 1, w ∼ GEM(α) π c i.i.d. ∼ GEM(β ) c ≥ 1 F j |w i.i.d. ∼ Mult(w) 1 ≤ j ≤ J Z j i |π, F j i.i.d. ∼ Mult(π F j ) 1 ≤ j ≤ J, 1 ≤ i ≤ n j
where GW (•) denotes the Gaussian Wishart distribution, depending on 4 hyperparameters, GEM(•) denotes the stick breaking distribution depending on 1 hyperparameter, and Mult(•) denotes the multinomial distribution (with weights as parameters). The choice of these distributions is convenient as it leads to a conjugate framework for conditional distributions, i.e., the conditional distribution of a hidden variable given the observations and the other hidden variables.

B.2 About the stick breaking construction

To build a prior for an infinite sequence of weights summing to one, we use the stick breaking construction. The topic weights π c are build using an hyperparameter α:

ν c i.i.d ∼ Beta(1, α) c = 1, 2, . . . π c = ν c c-1 ∏ i=1 (1 -ν i ) Init: π 1 = ν 1
Weights π c k are sampled similarlrly, with an hyperparameter β :

η c k i.i.d ∼ Beta(1, β ) k = 1, 2, . . . π c k = η c k k-1 ∏ i=1 (1 -η c i ) Init: π c 1 = η c 1
In the following, we rather deal with the sequence of ν and η c , as there are i.i.d. samples. We have, for each c and each k:

p(ν c ) = (1 -ν c ) α-1 p(η c k ) = (1 -η c k ) β -1

B.3 Likelihood

Keeping the notation of the main text, the likelihood of an observation V j i is given by:

p(V j i |Z j i , µ,Λ ) = ∞ ∏ k=1 ϕ(V j i |µ k ,Λ -1 k ) z j i,k , which results, for the complete data set V p(V|Z, µ,Λ ) = M ∏ j=1 n j ∏ i=1 ∞ ∏ k=1 ϕ(V j i |µ k ,Λ -1 k ) z j i,k
The distribution of latent allocation vectors (or local variables) F j and Z i j is given by an (infinite) multinomial distributions depending on weights w and π:

p(Z j i |F j , π) = ∞ ∏ c=1 ∞ ∏ k=1 (π c k ) z j i,k f j,c (5)
p(F j |w) = ∞ ∏ c=1 w f j,c c (6)
B.4 Variational approximations of the posterior distributions Let L = {Z j i }, {F j }, η, ν, µ,Λ be the set of (local and global) hidden variables. For any hidden variable U in L the corresponding q * (U) is given by ln

q * (U) = E L \U [ln p(V, Z, F, ν, η, µ,Λ )] + constant (7)
The complete joint distribution can be split in simpler terms:

p(V, Z, F, ν, η, µ,Λ ) = p(V|Z, µ,Λ )p(Z|F, η)p(F|ν)p(µ,Λ )p(η)p(ν) (8)

B.4.1 Optimal variational distribution for ν

The variational approximation of the trajectory cluster is a truncated stick breaking distribution. It is therefore defined by C max random variables, such that the first C max -1 are beta distributed, and the last one equals to one almost surely. A simple computation using equations ( 7) and ( 8) shows that for the C max -1 first terms, the variational distribution of ν c is given by:

q * (ν c ) ∼ Beta 1 + J ∑ j=1 E q F [1 F j =c ], α + J ∑ j=1 E q F [1 F j >c ] ,
where 1 denotes the indicator function.

B.4.2 Optimal variational distribution for η

Similarly to the clustering weights, a truncated stick breaking distribution is used. For the cluster c and the movement mode 1 ≤ k ≤ K max -1, we have:

q * (η c k ) ∼ Beta 1 + J ∑ j=1 E q F j [ f j,c ] n j ∑ i=1 E q Z j i [1 Z j i =k ], β + J ∑ j=1 E q F j [ f j,c ] n j ∑ i=1 E q Z j i [1 Z j i >k ]
where we have written F j = f j,1 , . . . , f j,Cmax .

B.4.3 Optimal variational distribution for (µ,Λ )

This computation is done in [START_REF] Bishop | Pattern Recognition and Machine Learning[END_REF]. For each movement mode 1 ≤ k ≤ K max , the optimal variational distribution is a G W (m k , ρ k , γ k ,W k ) with:

ρ k = ρ 0 + N k m k = ρ 0 m 0 + N k Vk ρ k γ k = γ 0 + N k W k = W -1 0 + N k S k + ρ 0 N k ρ 0 + N k ( Vk ) -m 0 )( Vk ) -m 0 )
where

N k = J ∑ j=1 n j ∑ i=1 E q Z j i [Z j i,k ] Vk = 1 N k J ∑ j=1 n j ∑ i=1 E q Z j i [z j i,k ]V j i S k = 1 N k J ∑ j=1 n j ∑ i=1 E q Z j i [z j i,k ](V j i -Vk )(V j i -Vk ) T

B.4.4 Optimal variational distribution for F j

The computation results in q * (F j ) being a multinomial distribution on the set [[1,C max ]] with the c-th weight proportional to:

E qν [log(π c )] + n ∑ i=1 Kmax ∑ k=1 E q Z j i [z j i,k ]E qη [log π c k ]
Here, one can see that

E qν [log(π c )] = E νc [log(ν c )] + c-1 ∑ i=1 E ν i [log(1 -ν i )]
This expectation can therefore be computed using the two following properties of the beta distribution. If X ∼ Beta(α 1 , α 2 ), then 1.

(1 -X) ∼ Beta(α 2 , α 1 ); 2. E[log X] = ψ(α 1 )ψ(α 1 + α 2 ), where ψ is the digamma function.

Of course, the same remarks hold for E qη [log π c k ].

B.4.5 Optimal variational distribution for Z j i

The computation results in q * (Z j i ) being a multinomial distribution on the set [[1, K max ]] with the k-th weight proportional to:

Cmax ∑ c=1 E q F j [ f j,c ]E qη [log π c k ] + 1 2 E q Λ k [log(detΛ k )] - 1 2 E q µ k ,Λ k (V j i -µ k ) T Λ k (V j i -µ k ) -log 2π

C Simulation design for experiments

We consider a framework with K = 12 movement modes, and C = 10 trajectory clusters. Each trajectory cluster c is characterized by: -Its initial movement mode density (the law of the first movement mode chosen during the cluster) denoted by π c 0 = (π c 0,1 , . . . , π c 0,K ); -Its transition densities between movement modes (the laws for the next movement mode after having been in a cluster k) denoted by π c 0 = (π c k,1 , . . . , π c k,K ) for all 1 ≤ k ≤ K; A movement mode k is characterized by:

-Its movement parameters µ k ,Γ k , Σ k , i.e. the OUP parameters; -Its sojourn time distribution, i.e. the distribution of time spent in this movement mode, denoted by d k .

In our example, d k is the p.d.f. of gamma distribution G (∆ k , 1 2 ) where ∆ k is a movement mode-specific parameter.

For each trajectory 1 ≤ j ≤ 500, we independently repeat the following procedure:

-Final time Sample3 a final time T j ; -Number of observations Sample4 a number of observed points n j ; -Observation times Sample independently υ 1 , . . . , υ n j points from a uniform U[0, T j ], and sort them to set times as t 1 := υ (1) , . . .t n j = υ (n j ) ; -Trajectory cluster Sample a trajectory cluster c with known probabilities w = (w 1 , . . . , w 10 ); -First movement mode Sample k 0 with probabilities π c 0 ; -First duration Set τ 0 = 0; Sample δ 0 from d k 0 , set τ 1 = τ 0 + δ 0 ; -First OUP sampling Set an initial value V 0 . Sample an OUP (V t ) τ 0 ≤t≤τ 1 starting at v 0 at all observation times between τ 0 and τ 1 and at τ 1 ; -Set i = 1; -While τ i < T j -Movement mode Sample k i with probabilities π c k i-1 ; -Duration Sample δ i from d k 0 , set τ i+1 = τ i + δ i ; -OUP sampling Sample a OUP (V t ) τ i ≤t≤τ i+1 starting at V τ i at all observation times between τ i and τ i+1 , and at τ i+1 ; -Set i = i + 1;

The resulting process is a continuous time hidden semi-Markov model [START_REF] Yu | Hidden semi-markov models[END_REF] whose emission densities are autoccorelated OU processes.

D Additional experimental results

Table 2 represents the clusters obtained by our algorithm and those of a DTW-based k-means and Table 3 compares the results. The clusters are ordered by the number of trajectories they contain. In Table 2, one may notice that the DTW-based k-means algorithm puts in the 2 first clusters trajectories that enter the traffic scheme from the South. One also notices a correspondance between the clusters 1,4 (our) → 1,2 (DTW), 2 (our) → 3 (DTW), 3 (our) → 4 (DTW), 5 (our) → 5 (DTW).

Fig. 2

 2 Fig.2Five realizations of bivariate Ornstein Ulhenbeck Processes of parameters µ, Γ and Σ , solution of the SDE eq.(1), starting at position v 0 = (0, 0) (top right corner).

Fig. 3 (

 3 Fig. 3 (Left) Simulation of a sequence of Ornstein Uhlenbeck Processes defining a trajectory. The red triangle marks the starting velocity. The trajectory is made of two movement modes (identified by their colors and line styles), characterized by two different sets of parameters for the velocity process from Eq. (1). (Right) Bivariate velocity processes, together with the instant τ 1 at which parameters change (vertical lines).

Fig. 4

 4 Fig. 4 Simulation of three different trajectories with the same model as in Figure 3. The left panel shows the position processes, the right panel shows the velocity processes. All processes start from (0, 0). (Left) Two trajectories (in plain lines) are made of the same movement modes so they belong to the same cluster. The other trajectory (dashed line) has a different distribution over movement modes so it belongs to another cluster. (Right) The trajectories are plotted on the bivariate velocity space, showing the three possible movement modes. Two of them (red and blue) are shared by the two trajectory clusters.

Fig. 5

 5 Fig.5Modelling framework for the two level clustering of our movement data set (consisting in J trajectories). Equations for the Ornstein-Uhlenbeck dynamics are given in section (2). Hidden movement mode dynamics are not fully specified, as a simplified independant model will be used for inference. Classical modelling framework would be a Markov chain or a semi Markov chain. Notations used in this Figure are detailed in Section 3.

Fig. 6

 6 Fig. 6 Graphical representation of the hierarchical structure of the simplified model of Section 3.1.

Algorithm 2

 2 Estimation of OUP parameters and final movement mode assignation Require: δ : a user-input minimal sequence length 1: ∀k, A[k] = {} Sets of Assigned sequences for each movement mode k 2: U = {} Set of Unassigned sequences 3: for each trajectory j do Filter out inconsistent sequences 4:

Fig. 7

 7 Fig. 7 Simulated trajectory data. Three movement modes are shared by two clusters. Each cluster corresponds to two different routes. The black triangle denotes the beginning of each trajectory. The first two movement modes are present in both clusters.

Fig. 9

 9 Fig.9Movement modes estimated distributions (on the velocity space). Each ellipsoid represents a movement mode distribution, the dot being the expected mean, and the contour the 50% centered quantile. The transparency shows the estimated relative weight. The dotted lines show the 0 values.

Fig. 10

 10 Fig. 10 50 trajectories of three estimated trajectory clusters, see the text for details. Colors correspond to those of Fig. 9.

Fig. 11

 11 Fig. 11 Typical trajectories from two estimated clusters. Trajectories starts at the south, and ends at the North East. In the first part of the trajectories, clusters are differentiated by the movement mode (South West to North East for on the left, South East to North West on the right). Then, all trajectories follow the same route.

Fig. 12

 12 Fig. 12 Example of result of the proposed second step in movement mode estimation. Both the velocity (left) and position (right) processes are shown (units are of no importance here). The red cross shows the first point of both processes. Colors correspond to different movement modes. The purple phase in the middle of the trajectory is inconsistent (as defined in the text). It is a transitory phase in the vessel's velocity process. After the second step, it is assigned as the beginning of the second movement mode.

Table 2 Table 3

 23 First five clusters for our and DTW-based kmeans algorithm. Blue color indicates the beginning of the trajectories, yellow color the end. Number of trajectories by cluster. The Normalized Mutual Information on those five clusters is equal to

Table 1

 1 Contingency table (counting the number of points) between true movement modes (in row) and estimated movement modes (in columns, the label assignment was made a posteriori), after the first step.

	True MM		Est. MM: Step 1			Est. MM: Step 2
		1	2	3	4	5	1	2	3
	1	3406	0	0 634	0	4026	12	2
	2	203 2506	0	0 31	213 2517	10
	3	0	74 1186	0	0	0	74 1186

In this purpose also, it should be pointed out that any stochastic approach is better than the fully deterministic CAVI algorithm, as discussed in[START_REF] Hoos | Stochastic local search: Foundations and applications[END_REF], for instance.

https://github.com/rtavenar/ushant_ais

In our example, it was sampled uniformly in the continuous interval[80, 500]

In our example, it was sampled uniformly in the discrete interval[50, 450]

Acknowledgements Authors would like to thank CLS (Collecte Localisation Satellites) and Erwan Guegueniat for providing the raw data that allowed building the AIS dataset used in this paper.

This work has been supported by DGA and ANR through (i) the ANR/Astrid SESAME project (ref: ANR-16-ASTR-0026) and (ii) the ANR MATS project (ref: ANR-18-CE23-0006).


Marcia Maia
email: marcia.maia@univ-brest.fr

Susanna Sichel

Anne Briais

Daniele Brunelli

Marco Ligi

Nicolas Ferreira

Thomas Campos

Berengere Mougel

Isa Brehme

Christophe C Hémond

Maia Marcia

Sichel Susanna

Briais Anne

Brunelli Daniele

Ligi Marco

Campos Thomas

Mougel Berengere

Hemond Christophe

Motoki Akihisa

Moura Denise

Scalabrin Carla

Pessanha Ivo

Alves Eliane

Ayres Arthur

Oliveira Pedro

Extreme mantle uplift and exhumation along a transpressive transform fault

rocks exposed there were not volcanic and postulated that the mechanism for their 50 formation was different from that building other oceanic islands [START_REF] Darwin | Geological Observations on the Volcanic Islands and Parts of South 218 America Visited During the Voyage of H.M.S. 'Beagle[END_REF] . The islets, formed by 51 variably serpentinised and mylonitised peridotites [START_REF] Melson | Paul Rocks: a 221 high--temperature, mantle--derived intrusion[END_REF][START_REF] Campos | Holocene tectonic uplift of the St Peter and St Paul Rocks (Equatorial Atlantic) consistent 224 with emplacement by extrusion[END_REF] , are currently uplifting at rates of 52 1.5 mm/yr [START_REF] Campos | Holocene tectonic uplift of the St Peter and St Paul Rocks (Equatorial Atlantic) consistent 224 with emplacement by extrusion[END_REF] . Previous work suggested that the exposure of such large volume of 53 ultramafic rocks resulted from an abnormally cold upper mantle 13--15 or cold 54 lithosphere [START_REF] Brunelli | Asthenospheric percolation of alkaline melts beneath the St. 233 Paul region (Central Atlantic Ocean)[END_REF] in the Equatorial Atlantic and that the islets were either part of an 55 extensional flexural ridge, as observed in other transform faults [START_REF] Tucholke | Fracture Zone[END_REF][START_REF] Hekinian | Submersible observations of Equatorial Atlantic mantle: the St Paul 226 Fracture Zone region[END_REF] or linked to 56 compression [START_REF] Palmiotto | Nonvolcanic tectonic islands in ancient and modern oceans[END_REF] . Our data reveal that the islets are part of a major uplift of the lithospheric 57 mantle due to a 10 Myr long period of transpression at the transform boundary [START_REF] Maia | Complex ridge--transform evolution and mantle exhumation at the St. 235 Paul fracture zone system, Equatorial Atlantic[END_REF] . 58 Similar push--ups exist along continental strike--slip faults [START_REF] Mann | Global catalogue, classification and tectonic origins of restraining and 239 releasing bends on active and ancient strike--slip fault systems[END_REF][START_REF] Dair | San Andreas fault geometry through the San Gorgonio Pass, 243 California[END_REF] but this is the first time 59 one has been fully mapped in the oceanic lithosphere near a present--day spreading 60 centre. Understanding the processes responsible for its formation will shed new light on 61 the behaviour of large transform boundaries and their response to changes in both local 62 and far--field stresses as well as on mechanisms leading to mantle exhumation at 63 spreading centres. 64 Four transform faults and three intra--transform ridge segments, cumulating a 630 km 66 offset, form the complex St. Paul transform system (Fig. 1 & Supplementary Fig. S1). 67 Tectonic patterns and lithology reveal that the St. Peter & Paul's islets are part of a 68 major, 200 km--long, 30 km--wide submarine shear zone that accommodates 69 transpressive stresses along the northern transform fault of the St. Paul system [START_REF] Maia | Complex ridge--transform evolution and mantle exhumation at the St. 235 Paul fracture zone system, Equatorial Atlantic[END_REF] (Fig. 70 1). Most dredged samples are breccias of mylonitic peridotite and ultramylonites that 71 underwent various degrees of serpentinisation and deformation (Supplementary Fig. 72 S2). The morphology and the distribution of the deformation in the transform fault allow 73 three tectonic domains in the shear zone to be defined: the Western Transfer Zone 74 (WTZ), the Central Transpressive Zone (CTZ) and the Eastern Shear Zone (ESZ) (Fig. 2). 75 76 An important segmentation of the transform fault with associated push--up ridges 77 resulting from dextral transpression characterizes the WTZ (Fig. 2a). Eastwards, a series 78 of left--stepping restraining bends and offsets in the transform, associated with a large 79 topography, form the distinct tectonic pattern of the CTZ. Thrust faults (Figs. 1, 2b &3) 80 mark the base of its south flank, overlapping the crust formed at the north intra--81 transform segment. These thrust faults, imaged by seismic reflection data (Fig. 3), are 82 associated with the exposure of mylonitised peridotites and deformed sediments at this 83 crustal contact. They form positive flower structures resulting from transpression at the 84 restraining offsets along the shear zone (Fig. 2b). The area is seismically very active and, 85 while most of the events have strike--slip focal mechanisms, a few compressive events [START_REF] Wolfe | Oceanic transform earthquakes with 245 unusual mechanism or locations: relation to fault geometry and state of stress in the 246 adjacent lithosphere[END_REF] 86 near the islets confirm the presence of thrust faults. The ~3500 m--high Atobá Ridge, 87 located at the largest offset of the transform fault, is a major push--up ridge in the centre 88 of the wider transpressive feature, marking the location of the most intense 89 deformation [START_REF] Mcclay | Analog models of restraining stepovers in strike--slip fault 248 systems[END_REF] . The gravity--derived density structure over this portion of the shear zone and the nature of the rocks sampled on both the islets and the flanks of the Atobá ridge 91 (Fig. 4 and Supplementary Figs. S2 &S3) imply the presence of a core of uplifted, 92 relatively unaltered high--density mantle beneath the islets and an anomalously thick 93 low density layer on the ridge flanks, especially the southern one. We suggest that the 94 thickening of the low--density layer derives from a local deepening of the 95 serpentinisation front. Accordingly, the 500°C isotherm, that can be considered the limit 96 of the serpentinisation process [START_REF] Ulmer | Serpentine stability to mantle depths and subduction--250 related magmatism[END_REF] , lies at a depth of 25 km below the Atobá ridge (Fig. 4) 97 i.e. deeper than the inferred maximum depths of the low--density layer, supporting the 98 idea that the low--density material is related to hydration of mantle rocks. In the ESZ 99 (Fig. 2c), the transform fault crosses a deep basin similar to the transform fault basins 100 observed elsewhere in the St. Paul system (Fig. 1 and Supplementary Fig. S1). The basin 101 morphology suggests the absence of significant transpression in this section of the shear 102 zone. Close to the intra--transform segment, the transform becomes a double fault. This 103 reveals that even in the more linear part of the transform domain the motion is not 104 purely strike--slip. 105 106

The oblique segments and the offsets disrupting the main trend of the western 107 transform boundary developed after several episodes of southward propagation of the 108 northeastern segment of the Mid--Atlantic Ridge (MAR) in the last 10 Myr, as revealed by 109 the lengthening of its abyssal hills. The absence of asymmetric faulting at this large--110 offset ridge--transform intersection implies that the robust segment receives enough 111 melt supply to counteract the cold--edge effect of the transform fault 23 (Figs. 1 & 2a). The 112 short northern intra--transform segment, although comparatively less robust, did not 113 retreat to accommodate the propagation of the MAR segment, leading to an overlap of 114 both spreading segments. The morphology of the intra--transform segment, with sub--parallel abyssal hills mapped on both flanks, and the basalt samples recovered on and 116 off axis are consistent with spreading dominated by volcanic processes for at least the 117 last 10 Myr. However, contrary to the MAR, this small segment did not increase in 118 length, reflecting the difference in melt supply between the segments. 119 120

The change in the tectonics of the northern boundary of the St. Paul system started ~11 121 Ma, when a 5° counterclockwise change in the spreading direction between the South 122 America and Nubia plates occurred, inducing extension at large--offset left--stepping 123 transforms in the Central and Equatorial Atlantic [START_REF] Bonatti | Flexural uplift of a lithospheric slab near the Vema transform (Central 208 Atlantic): Timing and mechanisms[END_REF][START_REF] Ligi | Oceanic broad multifault 254 transform plate boundaries[END_REF] . The origin of the St. Paul 124 transpressive shear zone relates to the way the spreading geometry locally adjusted to 125 this plate motion change. We construct an evolutionary scenario describing the 126 responses of the transform to a sequence of far and local stress changes (Fig. 5). The 11 127 Ma counterclockwise rotation locally resulted in the transtensional formation of a 128 flexural transverse ridge along the northern border of the transform fault. This same 129 event resulted in the synchronous growth of the Vema transverse ridge at 11°N [START_REF] Bonatti | Flexural uplift of a lithospheric slab near the Vema transform (Central 208 Atlantic): Timing and mechanisms[END_REF] . 130 Shortly after, the spreading ridge segments started to adjust to the new spreading 131 geometry, with the lengthening of the MAR segment, while the intra--transform segment 132 to the East remained roughly stable. As the MAR segment kept propagating south, the 133 western part of the transform started to segment and this non--uniform adjustment 134 resulted in the formation of restraining bends and offsets, inducing localized 135 transpression on the western portion of the flexural transverse ridge. The resulting 136 structure is bounded by low--angle thrusting along the external faults of the system (Fig. 137 2). The four lengthening events, marked by the steps in the transform boundary, 138 correspond to progressively slower propagation rates, as the ridge geometry adjusted to 139 the local stresses (Fig. 5). The oldest and largest offset corresponds to a fast displacement of the western transform segment from its main trend. This event created 141 the first compressive stresses at the transverse ridge. As the transform fault adjusted 142 southwards during the second event, the increase in the fault segment offset induced 143 higher stresses that begun to form the sigmoidal push--up Atobá ridge. Lengthening of ridge segments due to increased magma supply is well documented 159 along the MAR [START_REF] Gente | Characteristics and evolution of the segmentation of the Mid--Atlantic 256 ridge between 20°N and 24°N during the last 10 million years[END_REF] and may be an important mechanism here. Geochemical data 9, [START_REF] Schilling | Thermal structure of the mantle beneath the equatorial Mid--228 Atlantic Ridge: inference from the spatial variation of dredged basalt glass compositions[END_REF][START_REF] Voyer | The effect of primary 259 versus secondary processes on the volatile content of MORB glasses: An example from 260 the equatorial Mid--Atlantic Ridge (5°N-3°S)[END_REF] show 160 the existence of a mantle composition boundary within the St. Paul system, with an 161 enriched mantle beneath the MAR segment north of the system extending to the 162 northern intra--transform segment. The volcanic morphology of these segments possibly 163 reflects the nature of the mantle, which favours enhanced melting. 164 Continental shear zones where restraining step--overs result in transpressive flower--like 166 structures and push--up ridges, exposing mylonites and ultra--mylonites at the surface 167 are well documented [START_REF] Mann | Global catalogue, classification and tectonic origins of restraining and 239 releasing bends on active and ancient strike--slip fault systems[END_REF][START_REF] Dair | San Andreas fault geometry through the San Gorgonio Pass, 243 California[END_REF] . The observed structural and petrologic patterns at the St Paul 168 northern transform boundary are similar, with multi--segmented and sub--parallel faults 169 away from the step--ups and oblique faults close to them. Associated thrust faults uplift 170 large reliefs and expose deep rocks distributed along bands parallel to the main strike--171 slip fault zone, forming positive flower structures. At the centre of these features, the 172 higher topography is associated with a push--up block, oblique faults and intense 173 deformation [START_REF] Mann | Global catalogue, classification and tectonic origins of restraining and 239 releasing bends on active and ancient strike--slip fault systems[END_REF][START_REF] Mcclay | Analog models of restraining stepovers in strike--slip fault 248 systems[END_REF] . The MAR half--spreading rate (16 mm/yr) [START_REF] Demets | Geologically current plate motions[END_REF] yields convergence rates of 174 32 mm/yr at the step--overs and transpressive segments of the fault. This amount of 175 shortening is sufficient to uplift a 100 km--long block by more than 3500 m, with 176 present--day uplift rates of 1.5 mm/yr [START_REF] Campos | Holocene tectonic uplift of the St Peter and St Paul Rocks (Equatorial Atlantic) consistent 224 with emplacement by extrusion[END_REF] , comparable to those estimated for structures 177 observed along continental strike--slip faults [START_REF] Mann | Global catalogue, classification and tectonic origins of restraining and 239 releasing bends on active and ancient strike--slip fault systems[END_REF][START_REF] Dair | San Andreas fault geometry through the San Gorgonio Pass, 243 California[END_REF] . 178 179

The evolution of the northern St. Paul plate boundary, where the southward migration 180 of both the MAR segment and the transform fault is progressively accommodated 181 through bends, offsets and oblique structures over a wide shear zone, lead to a 182 particular situation in the mid--ocean ridge system, with long--lived regional--scale 183 transpression along the transform. Whereas most exposed mantle rocks near mid--184 ocean ridges result from extensional tectonics along detachment faults, our study 185 reveals that the peridotites and mylonites of the St. Peter & Paul's islets result from a 186 major uplift of the oceanic lithosphere along a transpressive shear zone due to local 187 stresses which have been active during a time much longer than usually necessary for a 188 plate boundary to adjust to a kinematic change. This may also be the case along other 189 major transform faults, although probably at a smaller scale. Here we show that even large--offset transforms, materializing large contrasts in lithosphere thickness, may 191 deform due to the local stresses induced by the increase in melt supply at the spreading 192 ridge. Mantle composition and temperature may therefore trigger the local response of 193 the lithosphere. To understand deformation at large transforms, mantle processes must 194 not be neglected. 195 The gravity data obtained during the cruise [START_REF] Maia | [END_REF] were processed in the conventional way 340 (Eötvos, drift and latitude corrections) to compute the free--air anomalies before being 341 merged into a grid, together with satellite--derived free--air anomaly data [START_REF] Sandwell | New global 417 marine gravity model from CryoSat--2 and Jason--1 reveals buried tectonic structure[END_REF][START_REF] Maia | The support mechanism of the young Foundation 420 Seamounts inferred from bathymetry and gravity[END_REF] to obtain 342 the free--air anomaly used in this study. The density model consists of three layers of 343 constant density: the sediment cover, the crust or altered mantle (or their density 344 proxy) and the normal mantle, respectively with densities of 2400, 2800 and 3300 345 kg/m 3 and 1030 kg/m 3 for the water. Modelled densities do not attempt to reproduce 346 the complex reality inferred from the petrology, but consider densities that may be close 347 to altered and serpentinised peridotites, mylonites and crustal rocks. To construct a 348 sediment thickness grid, we used the values derived from the interpretation of the 349 seismic lines obtained during the cruise as well as available data from older cruises and 350 extrapolated them to the neighbouring areas. The sedimentary infill reaches 300 to 500 351 m in the nearby basins, but most of the area presents only relatively thin sediment cover 352 (less than 100 m). As the seismic lines are spaced further apart than the echosounder 353 sampling, the resolution of the sediment grid is not as good as the bathymetry, so for the 354 model, we interpolated the high--resolution bathymetry [START_REF] Maia | [END_REF] onto a 1 km step grid, 355 compatible with the sediment thickness grid, after projecting all data from geographic 356 degrees to UTM kilometres. From the bathymetry and the sediment thickness, we 357 calculated the basement topography. To calculate the Mantle Bouguer anomaly, we 358 assumed a constant, 6 km--thick "crust" (or its proxy, a 2800 kg/m 3 density layer), using 359 the basement as the top of the layer. The gravity effect of this model was computed in 360 the Fourier domain using the multi--layer method, developed to account for both rapidly 361 varying topography and shallow water depths [START_REF] Maia | The support mechanism of the young Foundation 420 Seamounts inferred from bathymetry and gravity[END_REF] and subtracted from the free--air thickness variations from the assumed 6 km constant "crust", we removed the effect of 364 the cooling of the lithosphere. Several thermal models were tested: a simple age 365 model [START_REF] Parsons | An analysis of the variation of ocean floor bathymetry and 422 heat flow with age[END_REF][START_REF] Stein | A model for the global variation in oceanic depth and heat flow 424 with lithospheric age[END_REF] and passive flow models [START_REF] Morgan | Three--dimensional flow and temperature perturbations 426 due to a transform offset: Effects on oceanic crustal and upper mantle structure[END_REF][START_REF] Ligi | Three--dimensional passive mantle flow 429 beneath mid--ocean ridges: An analytical approach[END_REF] using different values for the thermal 366 conductivity. As differences were only minor between the model results and expressed 367 mostly as a broad regional trend, we chose to keep the age model. Moreover, the passive 368 flow models tend to overemphasize the cold--edge effect, which our bathymetry data 369 suggest is not so marked here, despite the large offset. The 2800 kg/m 3 layer thickness 370 displayed in Figure 4 is obtained by adding the computed thickness variations to the 6 371 km thick "crust". 372 373

Seismic processing 374

The seismic lines were acquired with a 24 channels streamer using two air guns with 375 105 and 85 cubic inches [START_REF] Maia | [END_REF] . All lines were processed using software Sispeed v5.5 376 developed by IFREMER. Common mid--point gathers were stacked and migrated using a 377 simple velocity model with a sound velocity of 1500 m/s assigned for both water and 378 sediments. After processing, interpretation was performed using Kingdom 8.5 software. 379

The quality of the seismic lines in this rough topography area was checked by 380 comparison between the sea bottom reflector derived from the seismic profiles and the 381 topography derived from the high--resolution multibeam bathymetry. Only in a few 382 areas, where the seismic lines closely parallel high scarps, significant differences due to 383 lateral echoes were found. A few seismic lines crossing the CTZ were migrated with a 384 more complex model for the velocity of the sound in the sediments and the basement, 385 accounting for the increase in the sound velocity with depth, in order to better identify 386 the reflectors in the basement. The model consists of a sediment layer with a velocity of 1750 m/s followed by a 300 m thick layer with a velocity of 2000 m/s and then a series 388 of 500 m--thick layers with a velocity increase of 500 m/s each. The last horizon, at a 389 depth of 5800 m from the top of the oceanic crust, seismically separates the crust from 390 the mantle with a velocity 7500 to 8000 m/s. The line shown in Figure 3 was processed 391 with this full migration model. 392 393

Crustal age modelling 394

Since the proximity of the magnetic equator hinders the identification of magnetic 395 anomalies, we analysed the evolution of this part of the St. Paul system through a 396 kinematic reconstruction using the most recently published poles of rotation [START_REF] Mueller | New constraints 431 on the late Cretaceous/Tertiary plate tectonic evolution of the Caribbean[END_REF] . 397 Theoretical crustal ages were calculated using the present day geometry of the mid--398 oceanic ridge and these poles. 399 400

Data availability 401

Gravity and seismic data were acquired during the COLMEIA cruise [START_REF] Maia | [END_REF] . The multibeam 402 bathymetry used in this work came mainly from the COLMEIA cruise [START_REF] Maia | [END_REF] , completed with 403 data from PRIMAR96 [START_REF] Gasperini | New data on the geology of the Romanche F.Z., equatorial Atlantic: 434 PRIMAR--96 cruise report[END_REF] and S7 (R/V Strakhov) [START_REF] Udintsev | Equatorial Segment of the Mid--Atlantic Ridge[END_REF] cruises. Satellite gravity data [START_REF] Sandwell | New global 417 marine gravity model from CryoSat--2 and Jason--1 reveals buried tectonic structure[END_REF]

  The following 144 episode, of smaller amplitude, formed a smaller restraining bend, continuing to sustain 145 the uplift of the Atobá ridge and forming the southern part of the CTZ. The most recent 146 evolution also displays southward deviations of the shear zone resulting in 147 transpressive deformation in the WTZ. The morphology, focal mechanisms, and 148 evidence for uplift in the Atobá Ridge show that transpression is still active in the 149 western part of the shear zone. Such an active deformation implies that the plate 150 boundary did not reach a steady--state geometry and that the changes in the spreading 151 geometry are still being accommodated, although at relatively lower present--day rates. 152 Despite their similar offset, the Vema and the northern St. Paul transform show a 153 different sequence of adjustments both starting with the same kinematic change. The 154 Vema transform reacted with a single short transtensive adjustment while St. Paul 155 records a multiple set of events over a larger time--span. A plausible cause of this 156 different behaviour is the influence of the Sierra Leone hotspot, which may have 157 enhanced the magma supply of the MAR segment just north of St Paul transform 9 . 158

Figure captions 298 299 Figure 1 .Figure 2 .Figure 3 .Figure 4 .Figure 5 .

 29912345 Figure captions 298

  13 14 15 16 17 km low-density layer thickness present day propagation growth of the push-up ridge en-echelon small pressure ridges third propagation step growth of the push-up ridge formation of the graben 5

  used to 404 complement the ship data are available at http://topex.ucsd.edu/grav_outreach/#grid. 405 The global relief dataset 38 ETOPO1 used to draw Figure 1 is available at 406 https://www.ngdc.noaa.gov/mgg/global/global.html. The data that support the findings 407 of this study are available from the corresponding author upon request. 408

Geochem. Geophys. Geosyst. 14, 4698-4717 (2013).

Supplementary information is available in the online version of the paper. Reprints and

Competing financial interests 295

The authors declare no competing financial interests. 296

Code availability 410

The FORTRAN code used to calculate the gravity anomalies [START_REF] Maia | The support mechanism of the young Foundation 420 Seamounts inferred from bathymetry and gravity[END_REF] is available from the 411 corresponding author upon request.


Traditional population genetic studies use genotypic or allelic frequency data obtained from several populations sampled at the same time point. However, temporal population genetics data offers a more powerful way to study complex dynamics, since we can follow the allele frequency changes through time in the population. The simulation part of our ABC framework was conducted with SLiM version 3.1 (Haller and Messer 2017). We simulated a simple scenario where a population of size N was sampled twice at the begin and at the end of a time interval. The proportion of the genome, the number of mutations under selection, and the mean of the gamma distribution that defines the distribution of fitness effects were sampled from prior probability density distribution. We propose the use of ABC-RF to co-estimate demography and selection in temporal population genomics data.

We ran ~50.000 simulations and calculated summary statistics to construct the reference table that was used to row the random-forests.

RESULTS

To test if the grown random-forest can recovers the expected population size N and the genetic load (selection parameter) we ran an additional 100 simulations with and without selection and used them as pseudo-observed data-set (POD).

Selection Neutral

Figure 1 .

 1 Figure 1. Schematic representation of the simulated model: a) one population sampled at two time points and; b) the genomic architecture that it a combination of neutral and beneficial mutation defined by the prior.

  Disentangling the effects of selection and demography is a longstanding difficulty in population genetics. Theoretical works show that selection affecting linked sites may bias inference of demography(Ewing and Jensen 2016; Schrider et al. 2016) and vice-versa. In this way, the ability to jointly make inference about both processes is a necessity not yet fully implemented. One proposed way to co-estimate neutral and selective parameters is with the use of ABC(Li et al. 2012); however, the traditional ABC is computationally expensive. The introduction of random forests in ABC, however, reduces the computational burden, making it possible to study complex dynamics as the presence of linked selection(Pudlo et al. 2016, Raynal et al. 2017).

Figure 2 .

 2 Figure 2. The prior (grey) and posterior (black) distributions of the population size N and the genetic load of one POD. The horizontal red lines are the true value, and the dashed lines represent the credible intervals. Top: the model with selection and, down: the model without selection.

Figure 3 .

 3 Figure 3. The comparison between ABC-RF and moment-based methods to estimate the population size. Dashed line represent the true value.