Loic Le Gratiet

Claire Cannamela

Kriging-based sequential design strategies using fast cross-validation techniques with extensions to multi-fidelity computer codes

Keywords: kriging, co-kriging, sequential design, cross-validation, resource allocations, multi-fidelity codes

Kriging-based surrogate models have become very popular during the last decades to approximate a computer code output from few simulations. In practical applications, it is very common to sequentially add new simulations to obtain more accurate approximations. We propose in this paper a method of kriging-based sequential design which combines both the error evaluation providing by the kriging model and the observed errors of a Leave-One-Out cross-validation procedure. This method is proposed in two versions, the first one selects points one at-a-time. The second one allows us to parallelize the simulations and to add several design points at-a-time. Then, we extend these strategies to multi-fidelity co-kriging models which allow us to surrogate a complex code using fast approximations of it. The main advantage of these extensions is that it not only provides the new locations where to perform simulations but also which versions of code have to be simulated (between the complex one or one of its fast approximations).

A real multi-fidelity application is used to illustrate the efficiency of the proposed approaches. In this example, the accurate code is a two-dimensional finite element model and the less accurate one is a one-dimensional approximation of the system.

Introduction

Kriging-based surrogate models have become very popular during the last decades to design and analyze computer experiments [START_REF] Sacks | Design and analysis of computer experiments[END_REF]]. The reader is referred to the books of [Stein, 1999], [START_REF] Santner | The Design and Analysis of Computer Experiments[END_REF] and [START_REF] Rasmussen | Gaussian Processes for Machine Learning[END_REF] for more detail about kriging models. Usually, in real applications, two stages are performed to surrogate a computer code with a kriging model. The first one consists in building a kriging model from simulations coming from an initial experimental design set. Many methods exist to build the initial design set, in order to ensure appropriate space filling properties, the reader is referred to [START_REF] Fang | Design and Modeling for Computer Experiments[END_REF] for a non-exhaustive review of them. The second stage consists in adding simulations sequentially at new design points which complete the initial one. The selection of the new points are usually based on criteria to improve the global accuracy of the kriging model and this will be our goal in this paper. To be complete, we mention that sequential kriging has also been widely used in optimization (see [START_REF] Jones | Efficient global optimization of expensive black-box functions[END_REF]], [START_REF] Picheny | Quantile-based optimization of noisy computer experiments with tunable precision[END_REF]) and to estimate probabilities of failure [START_REF] Bect | Sequential design of computer experiments for the estimation of a probability of failure[END_REF] Kriging models are a powerful tool to enrich an experimental design set since it provides through the kriging variance -also called predictor's Mean Squared Error (MSE) or variance of prediction -an estimation of the model MSE. Kriging literature provides lot of criteria usually based on the kriging variance for sequentially design the experiments [START_REF] Sacks | Design and analysis of computer experiments[END_REF]]. Furthermore, [START_REF] Bates | Experimental design and observation for large systems[END_REF] and [START_REF] Picheny | Adaptative designs of experiments for accurate approximation of a target region[END_REF] propose more efficient criteria by considering the Integrated MSE (IMSE). It consists in integrating the mean value of the MSE integrated over the input parameter space. We note though that the IMSE can be computationally expensive to assess, especially when the dimension increases. Although these criteria are efficient for many cases, they can suffer from an important flaw when the accuracy of the kriging model is not homogeneous over the input parameter space. Indeed, the kriging variance is determined by the distances between prediction and design points but not by the real model errors. To fix this important flaw, we can use the Empirical IMSE suggested in [START_REF] Sacks | Design and analysis of computer experiments[END_REF]] which evaluates the model errors through a test set. Nevertheless, in a complex computer code framework, it could be too expensive to consider an external test set and cross-validation (CV) based criteria are more significant. As an illustration [Kleijnen and van Beers, 2004] and [van Beers and Kleijnen, 2008] combine a bootstrapping and a CV procedure to evaluate the predictor's MSE. Although this method improves the classical approach, it still does not take into account the real model errors. We note that a strength of the method proposed by [Kleijnen and van Beers, 2004] is that it can be applied to others type of surrogate models than the kriging one.

The first focus of this paper is on sequential design to improve the accuracy of a kriging model. In particular, we propose new criteria combining the kriging variance and the Leave-One-Out CV (LOO-CV) errors. The CV errors allow for focusing the new observations on regions where the real model errors are large. Furthermore, thanks to the equations presented in [Dubrule, 1983], the LOO-CV equations are fast to compute and thus the suggested approach is not expensive.

Kriging has recently been extended to allow for the use of coarse versions of a complex computer code to improve the accuracy of its approximation. These so-called multifidelity surrogate models have become of growing interest. A first one was proposed in [START_REF] Craig | Constructing partial prior specifications for models of complex physical systems[END_REF] which is based on a linear regression formulation and was improved by [START_REF] Cumming | Small sample bayesian designs for complex high-dimensional models based on information gained using fast approximations[END_REF] through a Bayes linear formulation. The first multi-fidelity model using co-kriging was suggested by [Kennedy and O'Hagan, 2000]. Then, several works dealing with this model have been developed [Kennedy and O'Hagan, 2001], [START_REF] Higdon | Combining field data and computer simulation for calibration and prediction[END_REF], [START_REF] Reese | Integrated analysis of computer and physical experiments[END_REF], [START_REF] Qian | Bayesian hierarchical modeling for integrating low-accuracy and high-accuracy experiments[END_REF]. Defining sequential design strategies in a multifidelity framework is also of interest and is still an open problem. A method based on nested Latin hypercube designs is suggested in [START_REF] Xiong | Sequential design and analysis of high-accuracy and low-accuracy computer codes[END_REF]. However, it does not allow for adding few simulations (e.g. it cannot perform an one step at-a-time sequential design) and it does not take into account the accuracy of the coarse code versions and the time ratios between two code levels.

The second focus of this paper is on sequential design for co-kriging model. We adapt the new strategies suggested for the kriging model to the multi-fidelity co-kriging one. The strength of the proposed extensions is that they not only provide the new points where to perform new simulations but they also determine which version of code is worth being simulated. These new criteria take into account the computational time ratios between code versions. They are based on a proxy of the IMSE reduction and on an original result giving the contribution of each code on the total variance of the model. We note that sequential design in a multi-fidelity framework has also been applied for optimization purposes [START_REF] Forrester | Multi-fidelity optimization via surrogate modelling[END_REF] and [START_REF] Huang | Sequential kriging optimization using multiple-fidelity evaluations[END_REF].

The paper is organized as follows. First, we introduce the kriging model and present our CV-based sequential design strategies. We illustrate these strategies in tabulated functions. Secondly, we present the proposed co-kriging multi-fidelity model and the extensions of the previous strategies. Finally, we apply the sequential co-kriging approach to a mechanical example.

Kriging models and sequential designs

In this Section, we briefly introduce the kriging model and some of its classical sequential design criteria. Then, we will present our sequential strategies to enhance kriging models considering the region with large LOO-CV errors.

The Kriging model

The kriging model is a widely used method to surrogate the output of a computer code from few simulations [START_REF] Sacks | Design and analysis of computer experiments[END_REF]]. Let us denote by y(x) the output of the code at point x ∈ Q. Here, y(x) is a scalar and Q ⊂ R d . Furthermore, we denote by D = {x 1 , . . . , x n } the experimental design set and y n = y(D) the value of y(x) at points in D.

In a kriging framework, we set that the prior knowledges about the code can be modeled by a Gaussian process Y 0 (x). Usually, we consider a Gaussian process with mean of the form m 0 (x) = f ′ (x)β, with f ′ (x) = (f 1 (x), . . . , f p (x)) and with covariance function k 0 (x, x) = σ 2 r (x, x; θ). Then, the kriging equations are given by the Gaussian process Y 0 (x) conditioned by its known values y n at points in D:

Y n (x) = [Y 0 (x)|Y 0 (D) = y n ] ∼ GP (m n (x), k n (x, x)) (1)
where:

m n (x) = f ′ (x) β + r ′ (x)R -1 (y n -F β) (2)
and:

k n (x, x) = σ 2 r(x, x) -f ′ (x) r ′ (x) 0 F ′ F R -1 f (x) r(x) (3)
where ′ stands for the transpose, F are the values of f ′ (x) at points in D, r(x) is the correlation vector between D and x with respect to the correlation function r(x, x), R is the correlation matrix of D with respect to r(x, x) and β = (F ′ R -1 F) -1 F ′ R -1 y n is the usual least-squares estimates of β. The model parameters σ 2 and θ could be estimated by maximizing their Likelihood [START_REF] Santner | The Design and Analysis of Computer Experiments[END_REF] or with a cross-validation procedure [START_REF] Rasmussen | Gaussian Processes for Machine Learning[END_REF]. Furthermore, the Maximum Likelihood Estimate (MLE) of σ 2 is given by σ2

= (y n -F β) ′ R -1 (y n -F β)/(n -p).
1 point at-a-time Sequential design Now, let us suppose that we want to add a new point x n+1 in D in order to enhance the accuracy of the kriging model. From the kriging variance k n (x, x) -representing the model MSE -some sequential design methods have been derived [START_REF] Sacks | Design and analysis of computer experiments[END_REF], [START_REF] Bates | Experimental design and observation for large systems[END_REF] and [START_REF] Picheny | Adaptative designs of experiments for accurate approximation of a target region[END_REF]. A first one consists in adding x n+1 where the kriging variance is the largest (see [START_REF] Sacks | Design and analysis of computer experiments[END_REF]):

x n+1 = arg max x k n (x, x) (4)
However, as presented in [Kleijnen and van Beers, 2004], its performance is poor. Then, it has been improved with a criterion which consists in adding the new point which leads the most important IMSE reduction (see [START_REF] Bates | Experimental design and observation for large systems[END_REF] and [START_REF] Picheny | Adaptative designs of experiments for accurate approximation of a target region[END_REF]):

x n+1 = arg max x u∈Q k n (u, u) -k n+1 (u, u) du (5)
Here, the covariance kernels k n+1 (u, ũ) corresponds to the one of the Gaussian process Y n (u)

(1) conditioned by a new observation at x. Furthermore, the equation (3) shows that the kriging variance does not depend on the observations if we consider known the parameters σ 2 and θ. Therefore, in that case, k n+1 (u, u) can be computed without new simulations. We denote by MinIMSE this criterion. Finally, we also consider the criterion presented by [Kleijnen and van Beers, 2004] using a Jackknife estimator for the predictor's variance. Its principle is the following one. Let us consider m n,-i (x) the kriging mean built without the i th observation, the jackknife variance is given by:

s 2 jack (x) = 1 n(n -1) n i=1 (ỹ i -ȳ) 2 (6)
where ỹi = nm n (x) -(n -1)m n,-i (x) and ȳ = n i=1 ỹi /n. Then, we consider candidate points coming from a maximin LHS Design [START_REF] Fang | Design and Modeling for Computer Experiments[END_REF] and we add those which maximize the jackknife variance. We denote by KleiCrit this criterion.

q points at-a-time Sequential design

There is a natural way to extend these algorithms when the simulations can be performed simultaneously. Indeed, the covariance kernel k n+1 (x, x) of the Gaussian process Y n (x) conditioned by the new observation at point x n+1 can be computed without knowing y(x n+1 ) when we consider the model parameters σ 2 and θ as known. Then, from k n+1 (x, x), we can find a new point x n+2 where to perform a new simulation using the same criterion as in equation ( 5) and the kernel k n+2 (x, x). Thus, considering the parameters σ 2 and θ as known (they are fixed to their estimated values), we can determine with this procedure q good locations where to perform simulations. We call this method the liar sequential kriging. This idea is also presented in a framework of kriging-based optimization in [Ginsbourger et al., ].

LOO-CV based strategies for kriging sequential design

We present in this subsection new sequential-kriging strategies. The main difference between these new strategies and the previous ones is that they take into account the real model errors through the LOO-CV equations.

The proposed sequential methods are based on the following original proposition. It gives the closed form expressions of the LOO-CV equations where the model parameters β and σ 2 are re-estimated after each removed point. This result is already known when σ 2 is fixed [Dubrule, 1983] and [START_REF] Rasmussen | Gaussian Processes for Machine Learning[END_REF]. Notations: A i,i is the i th element of the main diagonal of A, A i is the i th row of the matrix A, A -i is the matrix A without its i th row, A -i,i is the i th column of A without its i th element, A i,-i = A ′ -i,i and A -i,-i is the matrix A without the i th row and column.

Proposition 1 Let us denote by Y n,-i (x) the Gaussian process Y 0 (x) conditioned by the values y n,-i = y(D) \ y(x i ). Then, the predictive mean of Y n,-i (x) at point x i is given by:

m n,-i (x i ) = y(x i ) -R -1 (y n,-i -F -i β-i ) i / R -1 i,i (7)
where

β-i = (F ′ -i K i F -i ) -1 F ′ -i K i y n,-i and K i = R -1 -i,-i -R -1 -i,i R -1 i,-i / R -1 i,i
. Furthermore, the predictive variance of Y n,-i (x) at point x i is given by:

k n,-i (x i ) = σ 2 -i / R -1 i,i + ς -i (x i ) (8)
where

ς -i (x i ) = R -1 F i / R -1 i,i ′ (F ′ -i K i F -i ) -1 R -1 F i / R -1 i,i and σ 2 -i = y n,-i -F -i β-i ′ K i y n,-i -F -i β-i /(n -p -1)
The previous proposition provides a powerful tool to compute the LOO-CV predictive means and variances. Indeed, several elements of the equations have been already computed during the models construction (e.g. the inverse of the matrix R). Consequently, the LOO-CV equations are fast to compute and can be easily recomputed at each step of the sequential strategy. We note that the originality of this result is the estimation of σ 2 -i . As we use the value of k n,-i (x i ) strongly depending on σ 2 -i in our forthcoming developments, it is important to well estimate it. Now, let us denote by e 2 LOO-CV = ((y

(x i ) -m n,-i (x i )) 2 i=1,...,n
the vector of the LOO-CV squared errors and s 2 LOO-CV = [k n,-i (x i )] i=1,...,n the vector of the LOO-CV variances. Furthermore, let us consider the Voronoi cells (V i ) i=1,...,n associated with the points (x i ) i=1,...,n :

V i = {x ∈ Q, ||x -x i || ≤ ||x -x j ||, ∀j = i}, i, j = 1, . . . , n (9)
In the remainder of this section, we present two strategies to sequentially add simulations which use e 2 LOO-CV , s 2 LOO-CV and V i . The intuitive idea of the suggested criteria is to enhance the predictive variance in the locations where the LOO-CV errors are important.

LOO-CV-based 1 point at-a-time Sequential design

Let us denote by x n+1 the new point that we want to add to D. We consider the point solving the following problem:

x n+1 = arg max x k n (x, x) 1 + n i=1 [e 2 LOO-CV ] i [s 2 LOO-CV ] i I x∈V i (10)
This criterion considers the predictor's MSE k n (x, x) adjusted with the LOO-CV errors and variances. For equivalent k n (x, x), the criterion favors the points close to an experimental design point with large LOO-CV errors. Furthermore, if two points are in the same Voronoi cell, the one with the largest predictor's MSE is considered. Therefore, a sequential strategy with this criterion focus on the regions of Q where the LOO-CV errors are the largest. We note that the standardization with s 2 LOO-CV is important since it is not necessary to enlarge the predictor's MSE in the regions where it is well or over estimated.

LOO-CV-based q points at-a-time Sequential design

We extend here the previous criterion for a q points at-a-time sequential design. First, we emphasize that the liar sequential kriging is not relevant for this new criterion. Indeed, conditioning on model parameters, with a liar method we can compute the kriging variances (k n+i (x, x)) i=1,...,q but not the LOO-CV equations. Therefore, we use another strategy to propose q new locations where to perform the simulations. This approach is proposed in [START_REF] Dubourg | Reliability-based design optimization using kriging surrogates and subset simulation[END_REF] in a different framework. The idea of the suggested method is to select the q best points with respect to the criterion (10) from N candidate points. These N candidate points are chosen with the following algorithm.

1. Generate N MCMC samples with respect to the probability density function proportional to k n (x, x) with a suitable Markov Chain Monte Carlo (MCMC) technique [START_REF] Robert | Monte Carlo statistical methods[END_REF].

2. Extract from these samples N representative points with a N -means clustering technique [MacQueen, 1967].

As presented in [START_REF] Dubourg | Reliability-based design optimization using kriging surrogates and subset simulation[END_REF] the use of this algorithm to select N candidate points in a kriging framework is efficient. Indeed, it allows us to concentrate the points at the modes of the kriging variance. In the proposed strategy, we always take N ≥ q and we choose from the N cluster centers (C i ) i=1...,N the q points where k n,adj (x, x) =

k n (x, x) 1 + n i=1 [e 2 LOO-CV ] i [s 2 LOO-CV ] i
I x∈V i is the largest. For the MCMC procedure, we use a Metropolis-Hastings (M-H) algorithm with a Gaussian jumping distribution. It is centered on the last sample point and has a standard deviation such that the acceptance rate is around 30% (see [START_REF] Robert | Monte Carlo statistical methods[END_REF]). Furthermore, we set N MCMC such that N MCMC ≫ N . For the N -means procedure, we choose the value of N with respect to the following criterion:

max N ≥q min x∈(C i ) i k n (x, x) (11)
where (C i ) i=1...,N are the cluster centers. This criterion prevents from having a cluster center in a region where the kriging variance is close to zero. Furthermore, if the number of clusters is too high, the cluster centers get away from the modes and consequently the value of min x∈(C i ) i=1...,N k n (x, x) decreases. Therefore, this criterion also prevents from having a number of clusters too large.

Sequential design in a multi-fidelity framework

A computer code can often be run at different levels of accuracy. In this case, it can be worth using low-accuracy versions of a code to improve its approximation. Co-kriging models are well suited to build such multi-fidelity surrogate models (see [Kennedy and O'Hagan, 2000] and [START_REF] Qian | Bayesian hierarchical modeling for integrating low-accuracy and high-accuracy experiments[END_REF]). In this section, we present the considered multi-fidelity co-kriging models and we extend the previous sequential design strategies in this framework. We note that, in a multi-fidelity framework, the search for the best locations where to run the code is not the only point of interest. Indeed, once the best locations are determined, we also have to decide which code level is worth being run. This will not only depend on the time-ratios between the code levels but also on the contribution of each code level to the total predictor's MSE.

Multi-fidelity co-kriging models

Let us suppose that we want to surrogate a computer code output y s (x) and that coarse versions of this code (y l (x)) l=1,...,s-1 are available. These codes are sorted by order of fidelity from the less accurate y 1 (x) to the most accurate y s-1 (x). All code levels are modeled by Gaussian processes (Y l (x)) l=1,...,s with respect to the following relationship with l = 2, . . . , s:

   Y l (x) = ρ l-1 Ỹ l-1 (x) + δ l (x) Ỹ l-1 (x) ⊥ δ l (x) Ỹ l-1 (x) ∼ [Y l-1 (x)|Y (l-1) = y (l-1) ] (12)
where

Y (l) = (Y 1 , . . . , Y l ) with Y l = Y l (D l ), y (l) = (y 1 , . . . , y l ) with y l = y l (D l
) and (D l ) l=1,...,s are the experimental design sets at level l with n l points and such that

D s ⊆ D s-1 ⊆ • • • ⊆ D 1 .
We note that the nested property is not necessary to build the model but allows for efficient parameter estimations. Furthermore, conditioning on parameters

β l , σ 2 l and θ l , δ l (x) is a Gaussian process of mean f ′ l (x)β l with f ′ l (x) = (f l 1 (x), . . . , f l p l (x)
) and covariance kernel σ 2 l r l (x, x, θ l ). By convention Y 1 (x) has the same distribution as δ 1 (x). The multi-fidelity co-kriging model at level l = 2, . . . , s is given by the following distribution:

Y l n l (x) = [Y l (x)|Y (l) = y (l) ] ∼ GP(µ l n l (x), k l n l (x, x)) (13)
with:

µ l n l (x) = ρl-1 µ l-1 n l-1 (x) + f l (x) βl + r ′ l (x)R -1 l (y l -F l βl -ρl-1 y l-1 (D l )) (14)
and:

k l n l (x, x) = ρ2 l-1 k l-1 n l-1 (x, x) + σ 2 l r l (x, x) -h ′ l (x) r ′ l (x) 0 H ′ l H l R l -1 h l (x) r l (x) (15
) where F l are the values of f ′ l (x) at points in D l , r l (x) is the correlation vector between D l and x with respect to the correlation function r l (x, x), R l is the correlation matrix of D l with re-

spect to r l (x, x), h ′ l (x) = µ l-1 n l-1 (x) f ′ l (x) and ρl-1 βl = (H ′ l R -1 l H l ) -1 H ′ l R -1 l y l with H l = [y l-1 (D l ) F l ] is the usual least-squares estimates of ρ l-1 β l
. Furthermore, the restricted maximum likelihood estimate of σ 2 l is given by σ2 [START_REF] Santner | The Design and Analysis of Computer Experiments[END_REF].

l = y l -H l ρl-1 βl ′ R -1 l y l -H l ρl-1 βl /(n l - p l -1)
Remark. The important property of this co-kriging model is that its MSE (15) provides through the term ρ 2 l-1 k l-1 n l-1 the contribution of the code level l -1 to the total predictor's MSE. Therefore, it can allow us to determine which code level is worth being simulated at a new location x.

Sequential design for multi-fidelity co-kriging models

The aim of this subsection is to extend the sequential kriging strategies proposed in Subsection 3.2 to the suggested multi-fidelity co-kriging model. These extensions are based on the variance decomposition property presented in Subsection 4.1 in equation ( 15) and on the following extension of Proposition 1:

Proposition 2 For l = 1, . . . , s, let us consider the multi-fidelity model ( 12) with D l ⊆ D l-1 ⊆ • • • ⊆ D 1 . We denote by ξ j the index of D j corresponding to the i th point x l i of D l with 1 ≤ j ≤ l. Then, if we note ε LOO-CV,l (x l i ) the LOO-CV error (i.e. real value minus predicted value) at level l and point x l i , we have:

ε LOO-CV,l (x l i ) = ρl-1 ε LOO-CV,l-1 (x l i ) + R -1 l y l -H l ρl-1 βl ξ l / R -1 l ξ l ,ξ l ( 16
)
where ε LOO-CV,1 (x l i ) is given by Proposition 1,

ρl-1 βl = (H ′ l,-ξ l K l H l,-ξ l ) -1 H ′ l,-ξ l K l y l -ξ l and K l = R -1 l -ξ l ,-ξ l -R -1 l -ξ l ,ξ l R -1 l ξ l ,-ξ l / R -1 l ξ l ,ξ l
. Furthermore, if we note σ 2 LOO-CV,l (x l i ) the variance of the LOO-CV, we have:

σ 2 LOO-CV,l (x l i ) = ρ2 l-1 σ 2 LOO-CV,l-1 (x l i ) + σ 2 l,-ξ l / R -1 l ξ l ,ξ l + ς l ( 17
)
with ς l = u 2 l H ′ l,-ξ l K l H l,-ξ l -1 , u l = R -1 l H l ξs /[R -1 l ] ξ l ,ξ l and: σ 2 l,-ξ l = y l -ξ l -H l,-ξ l ρl-1 βl ′ K s y l -ξ l -H l,-ξ l ρl-1 βl n l -p l -2 (18)
where

H l,-ξ l = [y l -ξ l F l,-ξ l ].
Proposition 2 provides closed form expressions for the LOO-CV errors and variances. From them, the LOO-CV equations are fast to compute and consequently they can be used in a sequential procedure with a low computational cost. Furthermore, since the experimental design sets are nested, we state that during the LOO-CV procedure at level l, the points are removed from all code levels. Finally, from these equations, we can adjust the co-kriging variances k l n l (x, x) l=1,...,s at each level using the same method as presented in equation ( 10).

1 point at-a-time sequential co-kriging. First, let us consider x new the point solving the problem:

x new = arg max x k s ns (x, x) (19)
Therefore, we want to compute a new simulation at point where the predictor's MSE is maximal. Now, let us consider two successive code levels l -1 and l. The question of interest is to estimate which of these two code levels is worth being simulated. First, thanks to the equation ( 15), we can deduce the contribution of each code levels to the predictor's MSE. Let us define the following notation for l = 2, . . . , s:

σ 2 δ l (x) = σ 2 l 1 -h ′ l (x) r ′ l (x) 0 H ′ l H l R l -1 h l (x) r l (x) (20)
and

σ 2 δ 1 (x) = k 1 n 1 (x, x)
. Then, we have:

k l n l (x, x) = l i=1 σ 2 δ i (x) l-1 j=i ρ 2 j ( 21
)
Let us consider that the parameters (θ l ) l=1,...,s define the characteristic length-scales of the kernels ((r l (x, x; θ l )) i=1,...,s (see [START_REF] Rasmussen | Gaussian Processes for Machine Learning[END_REF] p.83). Then, we can approximate the reduction of the IMSE after adding a new point x new at level l with the following formula:

IMSE l red (x new ) = l i=1 σ 2 δ i (x new ) l-1 j=i ρ 2 j d m=1 θ m i (22)
with θ l = (θ 1 l , . . . , θ d l ). Indeed, at each stage, σ 2 δ i (x new ) l-1 j=i ρ 2 j represents the contribution of the bias δ i (x) to the co-kriging variance and d m=1 θ m i represents the volume of influence of x new at level j. This criterion is justify by the fact that the reduction of IMSE l defined by IMSE l = Q σ 2 δ l (x) dx after adding a new point x new has the same order of magnitude than σ 2 δ i (x new ) times the volume of influence d m=1 θ m i of x new . Now, let us consider that the ratio of computational times between the codes y l (x) and y l-1 (x) equals B l/l-1 . It is worth running the code y l-1 (x) if B l/l-1 IMSE l-1 red (x new ) > IMSE l red (x new ), i.e. if the potential uncertainty reduction by running B l/l-1 times y l-1 (x) is greater than the one when we run one simulation on y l (x). From this criterion, we can deduce the following algorithm for an one at-a-time sequential co-kriging model taking into account both the computational ratios between the different code levels and the contribution of each level to the total co-kriging variance.

Algorithm 1 One point at-a-time sequential co-kriging 1: Find x new such that x new = arg max x k s ns (x, x) 2: for l = 2, . . . , s do

3: if σ 2 δ l (x new ) < IMSE l then 4:
Run y l-1 (x new )

5:

end for 6:

else 7: if IMSE l-1 red (x new )/IMSE l red (x new ) > 1/B l/l-1 then 8:
Run y l-1 (x new )

9:

end for 10:

end if 11:

end if 12: end for 13: if (l = s) then 14:

Run y l (x new ) 15: end if Remarks: Algorithm 1 evaluates for two successive code levels l -1 and l, which one is worth being simulated. It starts with the levels one and two, then two and three and so on. When it finds that the level l -1 is more promising than the level l, it stops the loop and simulate x new at code levels y 1 (x), . . . , y l-1 (x). Since the loop is defined from level 1 to level s, it favors simulations at low code levels. Therefore, it will tend to learn the coarse code versions before learning the accurate ones. We note that during the loop of the algorithm 1, the parameters are not re-estimated. In fact, they are re-estimated after adding the new point x new . Moreover, the first test σ 2 δ l (x new ) < IMSE l checks in averaged if the code level l at point x new is worth being run. Then, the test IMSE l-1 red (x new )/IMSE l red (x new ) > 1/B l/l-1 evaluates which code levels between l and l -1 is the most promising. Finally, if we consider that the code level l is more promising than the code level l -1, we confront it to the following code level l + 1. We note that the algorithm 1 is reiterated until a prescribed accuracy is reached or the computational time budget is spent.

1 point at-a-time sequential co-kriging with adjusted predictor's MSE. From Proposition 2, Algorithm 1 and Equation ( 22), we can extend the criterion (10) to the multifidelity co-kriging model. Let us consider the following quantity:

IMSE l red,adj (x new ) = l i=1 σ 2 δ i (x new ) 1 + n i j=1 (εLOO-CV,i(x i j )-ρ i-1 ε LOO-CV,i-1 (x i j )) 2 σ 2 LOO-CV,i (x i j )-ρ 2 i-1 σ 2 LOO-CV,i-1 (x i j ) × l-1 j=i ρ 2 j d m=1 θ m i ( 23
) with the convention ρ0 = 0 and (ρ i ) i=1,...,l is given by Proposition 2. In equation ( 23), the kriging variances σ 2 δ i (x) in equation ( 21) is replaced with the adjusted kriging variance presented in Subsection 3.2. We note that ε LOO-CV,i (x i j )ρi-1 ε LOO-CV,i-1 (x i j )

2 is the part of the LOO-CV squared error explained by the bias δ i (x) and σ 2 LOO-CV,i (x i j ) -ρ2 i-1 σ 2 LOO-CV,i-1 (x i j ) is the corresponding LOO-CV predictive variance. To adapt the adjusted co-kriging variance in a multi-fidelity framework, we just have to replace IMSE l red (x) with IMSE l red,adj (x) in the algorithm 1.

(q i ) i=1,...,s points at-a-time sequential co-kriging. In this paragraph, we propose an extension for the multi-fidelity model of the q points at-a-time sequential design presented in Subsection 3.2. Its principle is the following one. First, we select q l new points for the code y l (x) with the method presented in Subsection 3.2 "LOO-CV based q points at-a-time Sequential design". Then, we consider these points as known for the code y l-1 (x) and we select q l-1 new points for this code with the same method. We note that, as presented in Subsection 3.1, we can use a liar method to compute the new co-kriging variance without simulating y l-1 (x) at the q l new points. Finally, we repeat this procedure for all code levels from y l-2 (x) to y 1 (x). At the end of the procedure, we have l i=j q i new points at level j and we want to find the allocation {q 1 , . . . , q l } leading the largest potential uncertainty reduction and under the constraint of a constant CPU time budget. We note the CPU time budget T = l j=1 l i=j q i t j where (t i ) i=1,...,s represents the CPU times of codes (y i (x)) i=1,...,s . The algorithm 2 presents the suggested q points at-a-time sequential co-kriging.

Algorithm 2 (q i ) i=1,...,s points at-a-time sequential co-kriging 1: Set the budget T > 0 and the allocation {q 1 , . . . , q l } such that l j=1 l i=j q i t j = T 2: Set (N i MCMC ) i=1,...,l for the M-H procedures.

3: Generate N l MCMC samples with respect to k l n l (x, x). 4: Find the N l cluster centers (C l i ) i=1,...,N l such that N l = max N ≥q l min x∈(C l i ) i k l n l (x, x) 5: Select from (C l i ) i=1,.
..,N l the q l points (x l new,i ) i=1,...,q l where k l n l ,adj (x, x) is the largest. 6: for m = l -1, . . . , 1 do

7:

Compute k m nm+ l i=m+1 q i (x, x) with the new points (x j new,i ) i=1,...,q l j=m+1,...,l

8:

Generate N m MCMC samples with respect to k m nm+ l i=m+1 q i (x, x).

9:

Find the N m cluster centers

(C m i ) i=1,...,N m such that N m = max N ≥q m min x∈(C m i ) i k m nm+ l i=m+1 q i (x, x) 10:
Select from (C m i ) i=1,...,N m the q m points (x m new,i ) i=1,...,qm where k m nm+ l i=m+1 q i ,adj (x, x) is the largest. 11: end for In Algorithm 2, k l n l + s i=l+1 q i (x, x) corresponds to the kernel of the Gaussian process Y l n l (x) conditioned by the observations (x j new,i ) i=1,...,q s j=l+1,...,s when the parameters (σ 2 i ) i=1,...,l and (θ i ) i=1,...,l are considered as known (i.e. this corresponds to a liar method). Furthermore, k l n l + s i=l+1 q i ,adj (x, x) corresponds to the predictor's variance k l n l + s i=l+1 q i (x, x) adjusted with the LOO-CV errors and variances:

k l n l + s i=l+1 q i ,adj (x, x) = l i=1 σ 2 δ i + s i=l+1 q i (x) 1 + n i j=1 (εLOO-CV,i(x i j )-ρ i-1 ε LOO-CV,i-1 (x i j )) σ 2 LOO-CV,i (x i j )-ρ 2 i-1 σ 2 LOO-CV,i-1 (x i j ) × l-1 j=i ρ 2 j d m=1 θ m i (24)
where k 1 n 1 + s i=l+1 q i and σ 2 δ i + s i=l+1 q i (x new ) are deduced from the equation ( 15). We note that for the M-H procedures, we use a Gaussian jumping distribution with a standard deviation such that acceptance rate is around 30%. Furthermore, let us consider the following quantity

IMSE red,q = l i=1 r=1,...,q i σ 2 δ i (x i new,r ) l-1 j=i ρ 2 j d m=1 θ m i (25)
We consider the allocation {q 1 , . . . , q l } which solves the following optimization problem:

{q 1 , . . . , q l } = arg max {q 1 ,...,q l } IMSE red,q such that l j=1 l i=j

q i t j = T (26)
i.e. we look for the allocation leading the maximal uncertainty reduction. This optimization problem is very complex to solve and a sub-optimal solution will be often considered. Nevertheless, when the number of code levels and the budget T are low (e.g. s = 2 in our application) an exhaustive exploration of the allocation {q 1 , . . . , q l } can be performed. We are in that case in the presented application . Furthermore, we note that IMSE red,q is a proxy on the IMSE reduction when we add (x m new,i ) i=1,...,qm m=1,...,l at code levels (y m (x)) m=1,...,l .

In practical application, the algorithm 2 is reiterated until we reach a prescribed precision or the computational time budget is exhausted.

Applications

We compare in this Section the MinIMSE, KleiCrit and AdjMMSE criteria on toy examples and on an application concerning a spherical tank under pressure. We present both the cases of 1 point at-a-time and q points at-a-time sequential kriging. Then, we compare on the tank application, the suggested sequential kriging and co-kriging methods with s = 2 levels. The purpose of this section is to emphasize the efficiency of the LOO-CV-based criteria and to highlight that a multi-fidelity analysis can be worthwhile. Finally, for the multi-fidelity sequential co-kriging, we present the allocation of the simulations between the coarse code and the accurate one. We note that for the different examples, we compare the different methods given a prescribed computational time budget.

Comparison between sequential kriging criteria

In this subsection, the 1 point at-a-time sequential kriging criteria (MinIMSE, KleiCrit, Ad-jMMSE) are compared on three tabulated functions: Ackley, 1987]: Xian, 2001]: Michalewicz, 1992]:

• Ackley's function on [-2, 2] 2 [
f (x, y) = -20exp -0.2 x 2 + y 2 2 -exp cos(2πx) + cos(2πy) 2 + 20 + exp(1) • Shubert's function on [-2, 2] 2 [
f (x, y) = 5 k=1 kcos ((k + 1)x + k) 5 k=1 kcos ((k + 1)y + k) • Michalewicz's function on [0, π] 2 [
f (x, y) = -sin (x) sin x 2 π 20 -sin (y) sin y 2 π 20
The comparison is performed on a test set D test composed of n test = 1000 points uniformly spread on the input parameter space and from 50 different initial experimental design sets.

We compare the different methods with respect to the Normalized RMSE:

Norm RMSE = ntest i=1 y real (x i test ) -y pred (x) 2 /n test max x∈Dtest y real (x) -min x∈Dtest y real (x) (27)
where y real (x) is the real value of the output and y pred (x) the predicted one. The initial experimental design sets are LHS designs of 10 points optimized with respect to the S-optimality [Stocki, 2005]. From these designs, 50 sequential krigings are performed and the convergence of the mean and the quantiles of the Normalized RMSE are computed for the three criteria. Furthermore, after each added point, the parameters of the kriging model are re-estimated with a maximum likelihood method. These estimations are performed thanks to the R library 'DiceKriging' [START_REF] Roustant | DiceKriging, DiceOptim: Two R packages for the analysis of computer experiments by kriging-based metamodeling and optimization[END_REF].

Figure 1 illustrates the efficiency of the criterion AdjMMSE. Indeed, for the Shubert's and Michalewicz's functions, we see that the accuracy of the 1 point at-a-time kriging with this criterion is significantly better than the one of the others criteria (both in terms of mean and quantiles of the Normalized RMSE). In fact, these functions have the particularity to have non-stationarity in some areas of the input parameter space. Thus, the errors are more important in these locations and the suggested criterion focus the new points on it. Furthermore, the non-stationarity is particularly important for the Shubert's function. For this reason, the IMSE criterion performed very bad in that case. Indeed, this criterion is efficient for stationary functions (i.e. when the predictor's MSE well predicts the real model errors). In contrast, the jackknife predictor's MSE provided by the criterion KleiCrit manages to catch this non-stationarity and it performs better than the IMSE criterion. Moreover, we see that the performance of the AdjMMSE and IMSE criteria are equivalent for the Ackley's function. We note that the Ackley's function is perfectly stationary and the errors have thus the same order of magnitude over the input parameter space.

These examples illustrate the fact that our criterion is more efficient than the other criteria when the functions are non-stationary and it remains efficient even in the cases where the functions are stationary (its efficiency is equivalent to the one of the IMSE criterion).

Spherical tank under internal pressure example

In this section, we deal with an example about a spherical tank under internal pressure. We are interested in the von Mises stresses on the three points labeled in figure 2. Indeed, we want to prevent from material yielding which occurs when the von Mises stress reaches the critical yield strength.

The system illustrated in figure 2 depends on the following parameters:

• P (M P a) ∈ [30, 50]: the value of the internal pressure.

• R int (mm) ∈ [1500, 2500]: the length of the internal radius of the shell.

• T shell (mm) ∈ [300, 500]: the thickness of the shell.

• T cap (mm) ∈ [100, 300]: the thickness of the cap.

• E shell (GP a) ∈ [63, 77]: the Young's modulus of the shell material.

• E cap (GP a) ∈ [189, 231]: the Young's modulus of the cap material.

• σ y,shell (M P a) ∈ [200,300]: the yield stress of the cap material.

• σ y,cap (M P a) ∈ [400, 800]: the yield stress of the cap material.

The accurate code output y 2 (x) is the value of the von Mises stress provided by an Aster finite elements code (http://www.code-aster.org) modelling the system presented in figure 2. We use the notation x = (P, R int , T shell , T cap , E shell , E cap , σ y,shell , σ y,cap ). We note that the material properties of the shell correspond to high quality aluminum and the ones of the cap corresponds to steel from classical to high quality. Then, the coarse code output y 1 (x) is the value of the von Mises stress given by the 1D simplification of the tank (28) (it corresponds to a perfect spherical tank under pressure, i.e. without cap).

y 1 (x) = 3 2 (R int + T shell ) 3 (R int + T shell ) 3 -R 3 int P (28)
According to equation ( 28), the actual input dimension of y 1 (x) is three (it depends only on P , R int and T shell ) while a sensitivity analysis performed with a Sobol decomposition gives that the accurate code depends essentially on four parameters (P , R int , T shell and T cap ). Furthermore, the response is highly stationary. Therefore, only few points are necessary to well predict the output of the code. For these reasons, we can start the sequential strategies from an initial experimental design set with only 10 points. Thus, for the different comparisons, we use a S-optimal LHS design D 2 of 10 points for the code y 2 (x). For the coarse code y 1 (x), we start with a design D 1 of 20 points. It is created with the following procedure. First, we create a S-Optimal design D1 of 20 points. Second, we remove from D1 the 10 points that are the closest to those of D 2 . Finally, D 1 is the concatenation of D 2 and D1 (this procedure ensures the nested property D 2 ⊂ D 1 ). We note that the CPU time is around 1 minute for the accurate code and 10 -8 seconds for the coarse code. Nevertheless, to be in a more realistic case, we consider that the CPU time ratio between y 2 (x) and y 1 (x) equals B 2/1 = 10. Furthermore, each sequential procedure is performed with 40 different initial design sets. Then, the mean and the quantiles of probabilities 90% and 10% of the empirical Normalized MSE are computed from a test set composed of 1000 points uniformly spread on the input parameter space. Finally, for the M-H procedure, we use a Gaussian jumping distribution such that the acceptance rate is around 30% and we set N MCMC = 50000 (we use 5 000 samples for the the burn-in procedure of the M-H method, see [START_REF] Robert | Monte Carlo statistical methods[END_REF]). For the M-H procedure, we use the package R CRAN mcmc. We note that after each added points, the parameters of the kriging or co-kriging models are re-estimated with a maximum likelihood method.

The remainder of this section is organized as follows. First we compare the MSE of the 1 point at-a-time sequential kriging with the one of the q = 5 points at-a-time one. Second, we compare for a given CPU time budget the sequential kriging and cokriging strategies. In the forthcoming developments, the response i = 1, 2, 3 refers to the value of the von Mises stress at point i on figure 2.

Comparison between sequential kriging criteria

Figure 3 compares the different criteria of the 1 point at-a-time and the q = 5 points at-a-time sequential kriging. We see that the criteria MinIMSE and AdjMMSE give equivalent values for the MSE for the 1 point at-a-time procedure and they perform better than the KleiCrit criterion. There are equivalent since the output y 2 (x) is perfectly stationary. Nevertheless, the criterion AdjMMSE is the most efficient for the q = 5 points at-a-time procedure. Indeed, the 5 points provided by a liar method with the MinIMSE criterion are not necessarily those which maximize the reduction of the IMSE. The method suggested in section 4.2 seems to give better solution.

Comparison between kriging and co-kriging sequential analysis

In this section, we compare the sequential kriging strategy with the sequential co-kriging with respect to the AdjMMSE criterion. Figure 4 gives the convergence of the empirical normalized MSE for the response 1. We see that the sequential co-kriging performs better than the kriging one. Furthermore, at the beginning of the method, the proportion of runs for the accurate code is very low. Indeed, the coarse code and the accurate code are extremely correlated for this response (around 99%) and thus, during the sequential strategy, the bias between the two codes is well estimated. Then, when the coarse code is well approximated, the sequential strategy starts to run the accurate one (for a CPU time around 500).

Figure 5 gives the convergence of the errors for the response 2. For this response, the correlation between the coarse and the accurate code is around 80%. Therefore, the proportion of runs for the accurate code determined by the sequential strategy is more important than in Figure 4. Furthermore, we see that this proportion increases with the CPU time. It means that the sequential co-kriging improves the approximation of the coarse code at the beginning of the procedure and then focuses on the accurate code. As a result, we see that the sequential co-kriging strategy is substantially better than the kriging one.

Figures 4 and5 illustrate the efficiency of the sequential co-kriging when the coarse code bring information on the accurate code. For the response 3, the coarse code is weakly correlated with the accurate code (around 45%). This is due to the fact that the coarse code models the von Mises stress in a perfect spherical tank whereas the response 3 corresponds to the one in the cap. Figure 6 shows that in this case, the sequential co-kriging model manages to determine that the coarse code is not worth being simulated. Indeed, the proportion of runs for the accurate code is very low. Furthermore, it shows that the co-kriging sequential design performs as well as the kriging one when the coarse code is non-informative.

Finally, Figure 7 shows the efficiency of the (q 1 , q 2 ) at-a-time sequential co-kriging. We set in the algorithm 2 that T = q 1 + q 2 + 10q 2 = 120 where the CPU time of the coarse code is 1 and the one of the accurate code is 10. For the the sequential kriging, we use a q = 10 at-a-time sequential procedure. Furthermore, Figure 7 shows that at the beginning of the procedure, the sequential co-kriging focuses on the approximation of the coarse code whereas at the end the number of runs for the accurate code is maximal. We note that the allocation of runs for the accurate code in figure 7 agrees with the proportion of runs given in Figure 5.

The results of the sequential co-kriging on the different responses show that the criterion suggested in section 4.1 performs very well. Indeed, it is always better than the sequential kriging when the coarse code is informative and its performance is equivalent to it when the coarse code is not useful. Furthermore, the different proportions of runs for the accurate code emphasizes that the criterion accurately determines the contribution of each code to the total model error and the optimal run allocation between the accurate and the coarse codes.

Conclusion

This paper deals with sequential strategies for kriging and co-kriging models. The kriging model is used to approximate the output of a complex computer code and the co-kriging one allows to improve this approximation thanks to coarse versions of the code.

First, we have presented classical sequential criteria for the kriging model and we have suggested another criterion based on the Leave-One-Out cross validation errors. This criterion has allowed us to set the new observations at locations where the model error is important. The examples presented in the last section have highlighted the efficiency of the suggested criterion. Indeed, for non-stationary functions, it provides results significantly better than classical criteria and for stationary ones its performance is equivalent to them. We have also emphasized the performance of the suggested criterion on a real application. Furthermore, we show in the application that when the simulations can be performed in parallel, our method has performed better.

Second, we have presented the extension of our criterion to multi-fidelity co-kriging models. We have shown in the application that performing a multi-fidelity sequential co-kriging is worthwhile when the coarse code versions are informative (i.e. highly correlated with the accurate code). Furthermore, a strength of the proposed approach is that it performs as well as a sequential kriging when the coarse code versions are not informative. In fact, the proposed extension takes into account the contribution of each code to the total predictor's mean squared errors and it determines the best run allocation between accurate and coarse code versions given a CPU time budget. Comparison between q = 10 points at-a-time sequential kriging and (q 1 , q 2 ) points ata-time sequential co-kriging. On the left, the bold circles represents the mean of the empirical MSE for the sequential kriging and the bold triangles represent the one of the sequential cokriging. Furthermore, the thin circles and triangles represent the quantiles of probabilities 10% and 90% of the empirical MSE. On the right, the squares represent the median number of runs for the coarse code during the sequential co-kriging and the triangles represent it for the accurate code.

Figure 1 :Figure 2 :

 12 Figure 1: Comparison between 1 point at-a-time sequential kriging criteria on toy examples. The bold triangles represent the mean of the empirical MSE for the AdjMMSE criterion, the bold circles represent it for the MinIMSE criterion and the bold Crosses represent it for the KleiCrit criterion. Furthermore, the thin triangles, circles and crosses represent the quantiles of probabilities 10% and 90% of the empirical MSE.

Figure 3 :

 3 Figure3: Comparison between 1 point at-a-time sequential kriging criteria (on left) and q = 5 points at-a-time sequential kriging criteria (on right) on the spherical tank example. The bold triangles represent the mean of the empirical MSE for the AdjMMSE criterion, the bold circles represent it for the MinIMSE criterion and the bold Crosses represent it for the KleiCrit criterion. Furthermore, the thin triangles, circles and crosses represent the quantiles of probabilities 10% and 90% of the empirical MSE.

Figure 6 :

 6 Figure 6: Comparison between 1 point at-a-time sequential kriging and co-kriging on the response 3 of the spherical tank example with respect to the AdjMMSE criterion (on the left). The thick dashed line represents the mean of the empirical MSE for the sequential kriging and the thick solid line represents it for the sequential co-kriging. The thin lines represent the quantiles of probabilities 10% and 90% of the empirical MSE. Figure on the right represents the proportion of runs allocated to the accurate code.

Figure 7 :

 7 Figure7: Comparison between q = 10 points at-a-time sequential kriging and (q 1 , q 2 ) points ata-time sequential co-kriging. On the left, the bold circles represents the mean of the empirical MSE for the sequential kriging and the bold triangles represent the one of the sequential cokriging. Furthermore, the thin circles and triangles represent the quantiles of probabilities 10% and 90% of the empirical MSE. On the right, the squares represent the median number of runs for the coarse code during the sequential co-kriging and the triangles represent it for the accurate code.

Acknowledgements

The author particularly thanks Professor Josselin Garnier for supervising his work and for his fruitful guidance. He is also grateful to Dr. Gilles Defaux for providing the application.


Armel Menant

Guillaume Bertrand

Christelle Loiselet

Laurent Guillou-Frottier

Laurent Jolivet

Spatial and temporal distribution of Cu-Au-Mo ore deposits along the western Tethyan convergent margin: a link with the 3D subduction dynamics

Along the western Tethyan convergent margin, where Tertiary subduction history is well constrained, porphyry, epithermal and skarn ore deposits show a variable evolution of their spatial distribution.Using different and complementary database on European and Middle East ore deposits, threemetallogenic episodeshave been highlighted: (1) a late Cretaceous -Paleocene phase characterized by a copper mineralization within the Balkan chain and in the Kaçkar mountains (eastern Turkey), (2) an Eocene phase with a few copper ore deposits in eastern Turkey and small Caucasia and (3) an Oligocene -Neogene phase with a more southern distribution along the margin and mainly constituted by epithermal Au systems in the west (Carpathians, Rhodope, Aegean and western Turkey) and by porphyry copper deposits in the east (Zagros).

These changesare suspected to be controlled by complex and evolving subduction dynamics.Using paleogeographic tools, it turned out that, in the eastern Mediterranean area, the late Cretaceous -Paleocene and Oligocene -Neogenemetallogenic episodesarecoeval with a significant decrease of the Africa -Eurasia convergence rate, from about 1.5 to 0.4 cm/yr. Indeed, compressionaltectonics in the volcanic arc domain, associated with a high convergent rate, promote the storage of large volumes of metal-rich magma and the hal-00746488, version 1 -29 Oct 2012

Author manuscript, published in "AGU Fall Meeting, San Francisco : France (2012)" development of an extensive MASH (melting, assimilation, storage and homogenization) zone. When this convergence rate decreases, a stress relaxation occurs in the overriding crust, inducing the ascent of a sufficient flux of this fertile magma and allowing the formation of numerous mineralized systems within the upper crust.

The Au-rich Oligocene -Neogenemetallogenic episode in the eastern Mediterranean region is also correlated with an increase of mantle-derived and/or subduction-modified lithospheric mantle components in magmas. This feature may be a consequence of the emplacement of hot asthenosphere at shallow depth related to (1) the development of a wide back-arc region due to slab retreat such as in the Aegean domain and (2) a slab tear and/or a lithospheric delamination, suspected notably in the Carpathians and western Turkey where alkaline to shoshoniticvolcanism occurs.

As the behavior of the slab and asthenosphere below the upper plate seems to play a key-role in controlling the distribution of ore deposits, it is worth studying the dynamics of the 3D mantle flow related toslab retreat.Thus, 3D numerical models of subduction dynamics with realistic rheologies have been developed.Around the slab edges, the poloidal (i.e. in a vertical plane) and toroidal (i.e. in a horizontal plane) components of the mantle flow in subduction zone appear to depend on the slab rollback to plate velocity ratio. Heat and mass transfersinduced by such 3D mantle flow, promote thermal anomalies in back-arc domain, observed on seismic tomographic models and necessary to produce fertile magmatism. hal-00746488, version 1 -29 Oct 2012

  Spatial and temporal distribution of Cu-Au-Mo ore deposits along the western Tethyan convergent margin: a link with the 3D subduction dynamics A. Menant (1), G. Bertrand (1), C. Loiselet (2), L. Guillou-Frottier (1), L. Jolivet (1) (1) ISTO, Université d'Orléans -BRGM, Orléans, France ; (2) MonashUniversity, Melbourne, VIC, Australia Emplacement conditions of mineralized systems in subduction and post-subduction environments and the sources of metals such as Cu, Mo and Au have been considered in the past. However, despite their importance in exploration strategies at the continental scale, interrelationships between distribution of ore systems and subduction dynamics are still partly unclear.


Xavier Perrin

Francis Colas

Cedric Pradalier

Roland Y Siegwart

Cédric Pradalier

Learning to identify users and predict their destination in a robotic guidance application

User guidance systems are relevant to various applications of the service robotics field, among which: smart GPS navigator, robotic guides for museum or shopping malls or robotic wheel chairs for disabled persons. Such a system aims at helping its user to reach its destination in a fairly complex environment. If we assume the system is used in a fixed environment by multiple users for multiple navigation task over the course of days or weeks, then it is possible to take advantage of the user routine: from the initial navigational choice, users can be identified and their goal can be predicted. As a result of these prediction, the guidance system can bring its user to its destination while requiring less interaction. This property is particularly relevant for assisting disabled person for whom interaction is a long and complex task. In this paper, we implement a user guidance system using a dynamic Bayesian model and a topological representation of the environment. This model is evaluated with respect to the quality of its action prediction in a scenario involving 4 human users, and it is shown that in addition to the user identity, the goals and actions of the user are accurately predicted.

Introduction

Robots are more and more present in the daily life, not only in the industry but also at home as toys or as service robots such as vacuum cleaners. There is also a growing demand in the health-care domain for smart assistive device such as intelligent wheelchairs. The present paper is focused on an intelligent system designed to help the elderly or disabled people in their daily activities. For these people, moving in their houses or passing through doorways may represent challenging tasks. We developed a semi-autonomous robot for improving user mobility while minimizing the required input, i.e. having an interaction process adapted to low throughput devices such as single switches, sip and puff systems, brain machine interfaces, or simple voice recognition. More precisely, at each crossing, the robot proposes a direction of travel to the user who will then either agree or disagree. The better the propositions are, the faster and easier the human-robot interaction is.

In this work, a dynamic Bayesian network (DBN) is used in order to learn the habits of multiple human users of a robotic helper. By habits, we mean the succession of navigational tasks one executes in a known environment: in a retirement home for instance, one resident wakes up in his bedroom, goes to the bathroom, then to the common room where the breakfast is served, and then goes to other rooms for his daily activities. Another resident will have other preferred locations for his own activities. All these accumulated additional information allow the robotic assistant to help its current human user from the first movements of a new travel until the destination. From the system's point of view, the actual identity of its user is not a relevant information. A user is merely defined by his activity pattern.

In the next section, we describe related works in recognition techniques for goal, user, or activity. In section 3, we detail the developed DBN, represented as a graphical model. Section 4 presents experiments in simulation as well as with a real robot and their results. The final discussion appears in section 5.

Related Works

Our aim is to ease the navigation of users, therefore we will consider activity recognition only from the point of view of navigation. In this case, inferring the user's intention requires techniques for plan recognition, which are used in a broad variety of domains, such as motion prediction, speech understanding, video surveillance, and so on. Uncertainty is inherent in plan inference, as the robot does not know in advance the intended destination of a user. Furthermore, many ways can lead to the same destination while one way can lead to several places. Probabilistic reasoning techniques are used in almost all works (review in [START_REF] Carberry | Techniques for Plan Recognition[END_REF]), as they help to express and maintain the beliefs in the possible goal destinations.

In the particular domain of intention recognition for navigation, two aspects were studied recently: local intention recognition (immediate action or location in the vicinity of the wheelchair from uncertain input) and global intention recognition (goal destination from local decision). Our work focuses on this latter issue, assuming that the recognition of the immediate intention is solved. In our test, this will be achieved by an interaction device with low uncertainty, such as a joystick or a reliable speech recognition. Inferring a global intention requires the ability to localize the robot on an available map of the environment. In a discrete environment represented by adjacent cells, Verma and Rao [START_REF] Verma | Goal-Based Imitation as Probabilistic Inference over Graphical Models[END_REF] described a Dynamic Bayesian Network (DBN) in the form of a graphical model composed of a Partially Observable Markov Decision Process (POMDP) enhanced with the notion of goal locations. Relying on reinforcement learning techniques, the agent explores how it should behave in order to reach three possible goals. Based on this acquired knowledge, the system could infer the most probable goal from a set of possibly noisy observations. Taha et al. [START_REF] Taha | POMDP-based Long-term User Intention Prediction for Wheelchair Navigation[END_REF] achieved similar results by fusing the robot position and the possible goals as the new state definition in a common POMDP. Vasquez [START_REF] Vasquez Govea | Incremental Learning for Motion Prediction of Pedestrians and Vehicles[END_REF] described a growing hidden Markov model algorithm dedicated to continuous learning, clustering, and making inference about car motions in a parking lot or people motions in a hallway.

Topological maps were successfully used in the works of Taha [START_REF] Taha | POMDP-based Long-term User Intention Prediction for Wheelchair Navigation[END_REF]. This environment representation is compact and matches the human's natural description of a path better than metric representations [START_REF] Rawlinson | Ways to Tell Robots Where to Go -Directing Autonomous Robots Using Topological Instructions[END_REF] (e.g. "go on the left at the second crossing" instead of "go straight for 100 meters, then turn left"). Many techniques exist for the map construction, e.g. the generalized Voronoi diagram, and applications based on imprecise human drawings have also been reported [START_REF] Setalaphruk | Robot navigation in corridor environments using a sketch floor map[END_REF].

User models are often used in computer applications (e.g. intelligent help) or online search or e-commerce sites (e.g. book recommendation based on other customer's choices) [START_REF] Schickel-Zuber | Inferring User's Preferences using Ontologies[END_REF]. Based on databases containing information from numerous users and some observations from the current user activity, reasoning techniques infer the next user action and try to help him. With the aging population, some research domains are focused on the activity recognition in so-called smart environments, i.e. environments where sensors have been installed in order to monitor the activity of human beings. The accumulated data are used to train algorithms which later serves to determine human activities, supervise the user's condition and medication, or detect anomalies [START_REF] Hristova | Context-aware services for ambient assisted living: A case-study[END_REF][START_REF] Lu | Creating Robust Activity Maps Using Wireless Sensor Network in a Smart Home[END_REF][START_REF] Mori | Anomaly detection algorithm based on life pattern extraction from accumulated pyroelectric sensor data[END_REF][START_REF] Zheng | Human Activity Detection in Smart Home Environment with Self-Adaptive Neural Networks[END_REF]. Some researches try to further determine the attributes (coffee drinker, smoker) of multiple users based on their location [START_REF] Liao | Location-Based Activity Recognition using Relational Markov Networks[END_REF][START_REF] Matsuo | Inferring Long-term User Property based on Users' Location History[END_REF].

As a summary, the studies on activity recognition perform well in capturing the user habits, but they can only monitor ongoing activities, not help the user to perform a task. On the other hand, the studies on intention recognition for facilitating the user's motion do not try to learn the typical user's daily habits. In this article, we propose a system that learns the habits of multiple human users controlling a robot and exploit this knowledge for the navigational control of the robot. In case of an unknown user, the same system first infers the identity of the user from the first movement in the environment before being able to exploit this knowledge for helping its user specifically. The robotic system is given a topological representation of the environment a priori. While incrementally learning the global intentions of a known user, it anticipates the user's destination and proposes better actions, as regular patterns of actions are performed day after day by the same user.

Model description

The graphical model shown in figure 1 represents the dynamic Bayesian network which composes the core of the system. Our DBN is an extension of a Markov Decision Process (MDP), with the agent's State S t and the Action A t (lower part). The state S t is the observed pose of the robot, in terms of position in a given topological map (a node) and orientation. 1 The action A t is among the repertoire: forward, left, right, u-turn, and stop. An action always connects to nodes of the map. The MDP is enhanced by the notion of Goal state G t , indicating the current global goal, and the variable Goal reached R t , a Boolean value indicating whether the agent's current state S t is a goal state or not (P(R

t |G t S t ) = 1 iff G t = S t ). A U t-1 S A G U R G last State Action Goal User goalReached Lastvisitedgoal t-1 A R S G U t t t t t G last G t-1 R t-1 S t-1
Fig. 1 Our dynamic Bayesian model for learning the habits of multiple users.

The term P(A t |S t G t R t ) is specified separately according to whether or not the goal is reached. First, when R t = 0, P(A t |S t G t ) is the action model leading to each goal G t from any state S t as deduced from the topological map (i.e no reinforcement learning technique used). Then, when R t = 1, the probability of staying still P(A t = stop) is higher than the other actions, so that the user can easily stop when it reaches a goal location. The last term of the MDP, the transition model P(S t |S t-1 A t-1 ), is also deduced from the topological map.

With the goal of learning the user habits for multiple users, we added two more variables. First, the user is symbolized with the variable U t . Second, the variable G last , representing the last visited goal, makes possible the learning of the succession of place visits, as will be shown later. Both the goal model

P(G t |G t-1 R t-1 U t-1 G last )
and the user model P(U t |U t-1 G t-1 G last ) result from submodels, where predefined terms related to the persistence of a goal (keeping the same goal or switching to another one), resp. the persistence of a user, are combined with terms learned online representing the knowledge acquired by the system. These submodels are described in details in the appendix. When reaching a goal during an online supervised learning phase, the known user and the particular G t and G last are used for updating the histograms of the probability distributions of the goal and user models.

From the model shown in figure 1, we compute at each time step several probability distributions. After the learning phase, the user is not known any more. Nevertheless, given the known variables S t-1 , A t-1 , and G last , we can infer his identity by computing P(U t |S t-1 A t-1 G last ) using Bayes' rule and the law of total probability:

P(U t |S t-1 A t-1 G last ) (1) = 1 Z ∑ G t-1      P(G t-1 ) ∑ U t-1 [P(U t-1 ) P(U t |U t-1 G t-1 G last )] × ∑ R t-1 [P(R t-1 |G t-1 S t-1 )P(A t-1 |G t-1 S t-1 R t-1 )]     
where the term 1 Z corresponds to a normalization factor. For a guidance robot, we also want to infer the goal the robot should be aiming to or the action it should propose to the unknown user. These two distributions can be computed as follows:

P(G t |S t-1 A t-1 G last ) (2) = 1 Z ∑ G t-1    P(G t-1 ) ∑ R t-1   P(R t-1 |G t-1 S t-1 )P(A t-1 |G t-1 S t-1 R t-1 ) × ∑ U t-1 {P(U t-1 )P(G t |G t-1 R t-1 U t-1 G last )}      (a) (b)
Fig. 2 (a) The differential drive robot used in our experiments. (b) Example of action proposed to the user, who then either agrees of disagrees orally (use of a speech recognition software).

P(A t |S t-1 A t-1 G last ) (3) = 1 Z ∑ St                    P(S t |A t-1 S t-1 ) × ∑ G t-1         P(G t-1 ) × ∑ R t-1            P(R t-1 |G t-1 S t-1 )P(A t-1 |G t-1 S t-1 R t-1 ) × ∑ U t-1     P(U t-1 ) × ∑ Gt    P(G t |U t-1 G t-1 R t-1 G last ) × ∑ Rt P(R t |G t S t ) ×P(A t |G t S t R t )                                             
If the user is known by the robot (e.g. through manual or visual identification), we can use a Dirac distribution for P(U t-1 ) and recompute all the above equations with this additional information, the system giving back the learned user habits.

Experiments

We run experiments of our multi-user guidance robot in our laboratory environment, using a differential-drive robot (fig. 2a). Figure 3 shows that the topological decomposition of the environment is made of goal nodes and connecting nodes. The former ones represent either people's desks (D1-D6) or common rooms like the cafeteria (C), the printer room (P), the robot lab (RL), or the bathroom (B). Four users share the robot, each with a particular desk and typical sequence, all starting from the entrance E. User 1 executes the sequence D3-D4-P-C-D3-D2-D3, user 2 D4-RL-C-D3-D4-RL-D4, user 3 D2-C-D2-D3-D2-P-D2, and user 4 D6-RL-C-D4-D6-RL-D6 (fig. 4a-d). Additionally each user goes to the bathroom at a random point in his sequence. These sequences are repeated 20 times during the learning phase. In order to speed up this process, we used a simulation. However, for the tests on the 21th day, a real robot was used. 2 In a second phase, the real robot was used by each of the user, the robot having no preliminary knowledge about the current user. At each time step, the system computes the probability distribution over the user, the goal, and the action with, respectively, equations 1, 2, and 3. This last distribution is used to propose an action to the user (see fig. 2b), who can confirm or not. In this precise experiment, this was achieved through a speech recognition system. Results of the user recognition are shown in figure 5, which present the evolution of the user probability distribution during a trial of the test phase. Starting from a uniform distribution, the actual user is correctly inferred as soon as a revealing action, mostly leading towards his desk, is executed. In the scenarios, the different user habits shared some common goals, or even some common sequences of goals. These similarities can be seen in figure 5, where the probability in the most probable user decreases and the one(s) for other user(s) increases (e.g. fig 5b, the sequence between goals D3 & D4 is shared by users 1 & 2). We can also notice that when a user goes to the bathroom, his inferred probability decreases slightly. As every user can go to the bathroom at any time, this is not a discriminant observation, the models tending thus slightly to a uniform due to the transition in the user submodel.

For quantifying our model, we compare the model prediction of the user, the goal, and the action with the real values. Based on their distributions computed at each time step, we can check if its maximum is the actual user, intended goal and proper action proposition. We thus introduce a measure µ max , being the mean of the number of times the max of a distribution matches the real value. For example, the formula for the goal is

µ max = 1 N ∑ N t=1 δ ( Ĝt , argmax G t P(G t |S t-1 A t-1 G last ))
, with N the number of steps, Ĝt the intended goal, and δ (a, b) the Kronecker function. 1 Model performance for the user, goal, and action predictions of each of the four single user scenarios (left columns) and the two taking-over scenarios (right columns).

The results are displayed in table 1. As can be seen when comparing with figure 5, the user are properly recognized, the differences being explained by the amount of steps until a revealing action. Concerning the goals and the actions, their probabilities are related to the learned habits of the users and are also influenced by the geographical location of the visited goals. As an example, user 3 visited three goals on the right when leaving his office, while the fourth one was on the left. Accordingly, the probability of going right, resp. of the three right goals were higher than going left to the bathroom. Furthermore, we can notice that the action inference is much better than the goal inference. Indeed, most of the time, the goal confusion is between two nearby goals that share a part of the sequence. As a consequence, the goal could be erroneously inferred during the first steps of the travel, while the action inference will be wrong only at a branching node.

Finally, we want to assess the robustness of our system. To this end, we run scenarios where two users took over the robot when user 2 reaches the cafeteria (C). First, the results with user 3 are displayed in figure 6a. As can be seen, while leaving C, the robot thought user 2 was going to D3, but also increased the probability in user 1, 3, and 4 as they all go in that direction after being in C. At each revealing action, the probability in specific users changes. But when reaching D2, user 3 starts being the most probable one. The increasing probability of user 1 is explained by the shared common goals with user 3, and starts to decline as soon as user 3 comes back to D2. In figure 6b, the results when user 4 takes the robot from user 2 are displayed. As user 4 first goes to the desk of user 2, this sequence is not revealing enough because only the former action differs from the habit of user 2. Then, going to the toilets does not bring any information about the user's identity. Finally, when user 4 reaches his desk, the remaining sequences are specific enough for the system to infer its user correctly. Overall, despite a performance loss in both goal and user recognition, table 1 shows that the actions proposed to the user are still relevant.

Conclusion

In this paper, we introduced a probabilistic model able to learn the daily habits of different users of a guidance robot. When the robot accumulated enough knowledge about a particular user, it is able to actively help its user to reach a goal destination. As shown in the experiments starting with an unknown user, the same model is able to infer the identity of the user after just a few revealing actions have been made, or specific goals have been reached. Having then a high confidence in a particular user, the system is again able to actively help its user. Furthermore, experiments where a different user takes over the guidance robot from an initial user have shown the ability of the system to recognize the switch between the users. We tested our system with a reasonable environment size and user number, but it would be interesting to test its robustness to an increased user number or to complexer user patterns. If leading to pattern overlaps, the system should still be able to help the users as they would go to the same location. However, the extreme case would be if there is a uniform distribution on the goals when leaving a particular one (e.g. the bathroom). For these cases, discretizing the user sequences over different periods of the day or increasing the history length should help.

We have P(G t-1 R t-1 ) which is a uniform distribution. Then, P(G t |G t-1 R t-1 ) defines the persistence of a goal: if a goal is not reached, the probability of having the same goal in mind is much higher than changing, whereas when the user reaches a goal, the other places become equally probable with a preference of staying in the current place. This term is predefined. Finally, the P(U t-1 G last |G t ) distribution is derived from a histogram learned by the system when a user stops at a goal location. From this model, we compute the probability distribution P(G t |G t-1 R t-1 U t-1 G last ), which is then used in the main model as the goal model, as follows:

P(G t |G t-1 R t-1 U t-1 G last ) ∝ P(G t |G t-1 R t-1 )P(U t-1 G last |G t )
The user model is build in a similar manner than the goal model and is the following: P(U t U t-1 G t-1 G last ) = P(U t-1 ) P(U t |U t-1 ) P(G t-1 G last |U t ) Again, P(U t-1 ) is a uniform distribution. P(U t |U t-1 ) describes the persistence of the human user, i.e. keeping the same user is more probable than switching between users. P(G t-1 G last |U t ) is also derived from a histogram, which is also learned when a user reaches a goal location. The probability distribution computed in order to build the user model is P(U t |U t-1 G t-1 G last ):

P(U t |U t-1 G t-1 G last ) ∝ P(U t |U t-1 )P(G t-1 G last |U t )

Fig. 3

 3 Fig. 3 Topological map of our laboratory used for the experiments with an entrance (E), desks (D1-D6), a printer room (P), a cafeteria (C), a robot lab (RL), and toilets (T).

Fig. 4

 4 Fig. 4 Habit description of the four users. The succession of visited places is marked with increasing numbers. Each user visits the bathroom (B) at random among the sequence.

3

 3

Fig. 5

 5 Fig. 5 User probability when (a) user 1, (b) user 2, (c) user 3, (d) user 4 was steering the robot. The dotted lines indicate the steps where a user stopped at a goal.

Fig. 6

 6 Fig. 6 Scenarios where (a) user 3 and (b) user 4 take the robot from user 2 in the cafeteria. The user probability is displayed. The dotted lines indicate the steps where a user stops at a goal. The black region indicates when a new user is steering the robot.

In relatively small and not too dynamic environments as the ones considered for this work, localization can be considered as a solved issue.

In simulation, the robot moves instantaneously to the next node, the inputs to our DBN S t-1 , A t-1 , and G last being identical to when using the real robot.

The software performed well enough to assume a perfect recognition of the user agreement ("Yes") or disagreement ("No").

Acknowledgments

This work has been supported by the Swiss National Science Foundation NCCR 'IM2' and by the EC-contract number BACS FP6-IST-027140. This paper only reflects the authors' view and funding agencies are not liable for any use that may be made of the information contained herein.

Appendix: Submodels for the goal and user models

In this appendix, we first describe the submodel used for building the goal model in the main system, then the submodel for the user model. In the goal model P(G t |G t-1 R t-1 U t-1 G last ), the goal at time t is dependent from a lot of variables issued from the previous time step. Some relations, e.g. between G t and G last , are user-dependent and some others, e.g. between G t and R t , can be predefined. A submodel is a Bayesian model on his own, with specific relation among the variables, which is then used in order to infer the probability distribution over some variables given some other known variables. In the case of the goal model, our submodel is defined as follows: