Loic Le Gratiet
email: loic.le-gratiet@cea.fr

Bayesian analysis of hierarchical multi-fidelity codes

Keywords: surrogate models, co-kriging, multifidelity computer experiment, Bayesian analysis

This paper deals with the Gaussian process based approximation of a code which can be run at different levels of accuracy. This method, which is a particular case of co-kriging, allows us to improve a surrogate model of a complex computer code using fast approximations of it. In particular, we focus on the case of a large number of code levels on the one hand and on a Bayesian approach when we have two levels on the other hand.

The main results of this paper are a new approach to estimate the model parameters which provides a closed form expression for an important parameter of the model (the scale factor), a reduction of the numerical complexity by simplifying the covariance matrix inversion, and a new Bayesian modelling that gives an explicit representation of the joint distribution of the parameters and that is not computationally expensive.

A thermodynamic example is used to illustrate the comparison between 2level and 3-level co-kriging.

Introduction

Large computer codes are widely used in science and engineering to study physical systems since real experiments are often costly and sometimes impossible. Nevertheless, simulations can sometimes be costly and time-consuming as well. In this case, conception based on an exhaustive exploration of the input space of the code is generally impossible under reasonable time constraints. Therefore, a mathematical approximation of the output of the code -also called surrogate or metamodel -is often built with a few simulations to represent the real system. 1 Gaussian Process regression is a particular class of surrogate which makes the assumption that prior beliefs about the code can be modelled by a Gaussian Process. We focus here on this metamodel and on its extension to multiple response models. The reader is refered to [START_REF] Santner | The Design and Analysis of Computer Experiments[END_REF] and [START_REF] Rasmussen | Gaussian Processes for Machine Learning[END_REF] for further detail about Gaussian Process models.

Actually, a computer code can often be run at different levels of complexity and a hierarchy of levels of code can hence be obtained. The aim of this paper is to study the use of several levels of a code to predict the output of a costly computer code.

A first metamodel for multi-level computer codes was built by Kennedy and O'Hagan [START_REF] Kennedy | Predicting the output from a complex computer code when fast approximations are available[END_REF] using a spatially stationary correlation structure. This multi-stage model is a particular case of co-kriging which is a well known geostatistical method. Then, Forrester et al. [START_REF] Forrester | Multi-fidelity optimization via surrogate modelling[END_REF] went into more detail about the estimation of the model parameters. Furthermore, Forrester et al. presented the use of co-kriging for multi-fidelity optimization based on the EGO (Efficient Global Optimization) algorithm created by Jones et al. [START_REF] Jones | Efficient Global Optimization of Expensive Black-Box Functions[END_REF]. A Bayesian approach was also proposed by Qian and Wu [START_REF] Qian | Bayesian Hierarchical Modeling for Integrating Low-accuracy and High-accuracy Experiments[END_REF] which is computationally expensive and does not provide explicit formulas for the joint distribution of the parameters.

This paper presents a new approach to estimate the parameters of the model which is effective in the case of non-spatial stationarity and when many levels of codes are available. In particular, it provides a closed form expression for the estimation of the scale factor which is new and of great practical interest. Furthermore, this approach allows us to consider prior information in the estimation of the parameters. We also address the problem of the inversion of the co-kriging covariance matrix when the number of levels is large. A solution to this problem is provided which shows that the inverse can be easily calculated. Finally, it is known that with a non-Bayesian approach, the variance of the predictive distribution may be underestimated [START_REF] Kennedy | Predicting the output from a complex computer code when fast approximations are available[END_REF]. This paper suggests a Bayesian modelling different from the one presented in [START_REF] Qian | Bayesian Hierarchical Modeling for Integrating Low-accuracy and High-accuracy Experiments[END_REF] which provides an explicit representation of the joint distribution for the parameters and avoids prohibitive implementation.

3 Building a surrogate model based on a hierarchy of s levels of code

Let us assume that we have s levels of code z 1 (x), . . . , z s (x), x ∈ R d , d > 0. For all t = 1, . . . , s the t th scalar output z t (x) is modelled by z t (x) = Z t (x, ω) where Z t (x, ω), ω ∈ Ω is a realization of the Gaussian process Z t (x). We will introduce below a consistent set of hypotheses so that the joint process (Z t (x)) x∈R d ,t=1,...,s is Gaussian given a certain set of parameters. Kenney and O'Hagan [START_REF] Kennedy | Predicting the output from a complex computer code when fast approximations are available[END_REF] suggest an autoregressive model to build a metamodel based on a multi-level computer code. Hence, we have a hierarchy of s levels of code -from the less accurate to the most accurate -and for each level, the conditional distribution of the Gaussian process Z t (x) knowing Z 1 (x), . . . , Z t-1 (x) is entirely determined by Z t-1 (x). Let us introduce here the mathematical formalism that we will use in this paper. Q ⊂ R d is a compact subset of R d called the input space or the domain of interest. For t = 1, . . . , s, D t = {x (t) 1 , . . . , x (t) nt } is the experimental design set at level t containing n t points in Q. Let Z t = Z t (D t ) = (Z t (x (t) 1 ), . . . , Z t (x (t) nt )) T be the random Gaussian vector containing the values of Z t (x) for x ∈ D t where T stands for the transpose. Let Z = (Z T 1 , . . . , Z T s ) T be the Gaussian random vector containing the values of the processes (Z t (x)) t=1,...,s at the points of the design sets (D t ) t=1,...,s . We assume here that the code output is observed without measurement error. The column vector of responses is written z = (z T 1 , . . . , z T s ) T , where z t = (z t (x (t) 1 ), . . . , z t (x (t) nt )) T is the output vector for the level t. If we consider Z s (x), the Gaussian process modelling the most accurate code, we want to determine the predictive distribution of Z s (x 0 ), x 0 ∈ Q given Z = z, i.e. the following conditional distribution:

[Z s (x 0 )|Z = z].
We assume the Markow property introduced in [START_REF] Kennedy | Predicting the output from a complex computer code when fast approximations are available[END_REF]:

Cov(Z t (x), Z t-1 (x ′ )|Z t-1 (x)) = 0 ∀x = x ′ (1)
This means that if Z t-1 (x) is known, then nothing more can be learnt about Z t (x) from any other run of the cheaper code Z t-1 (x ′ ) for x ′ = x. This assumption leads to the following autoregressive model:

Z t (x) = ρ t-1 (x)Z t-1 (x) + δ t (x) t = 2, . . . , s (2)
where δ t (x) is a Gaussian process independent of Z t-1 (x), . . . , Z 1 (x) and ρ t-1 (x) represents a scale factor between Z t (x) and Z t-1 (x). We assume that

ρ t-1 (x) = f ρ t-1 (x) T β ρ t-1 , t = 2, . . . , s, where f ρ t-1 (x) = (f 1 ρ t-1 (x), . . . , f q t-1 ρ t-1 (x))
T is a vector of q t-1 regression functions -generally including the constant function :

x ∈ Q → 1 -and β ρ t-1 ∈ R q t-1 .
Conditioning on parameters σ t , β t and θ t , δ t (x) is assumed to be a Gaussian process with mean f t (x) T β t , where f t (x) is a p t -dimensional vector of regression functions, and with a covariance function of the form

c t (x, x ′ ) = cov(δ t (x), δ t (x ′ )) = σ 2 t r t (x -x ′ ; θ t ),
where σ 2 t is the variance of the Gaussian process and θ t are the hyper parameters of the correlation function r t . Moreover, conditioning on parameters σ 1 , β 1 and θ 1 , the simplest code Z 1 (x) is modelled as a Gaussian process with mean f 1 (x) T β 1 and with covariance function c 1 (x, x ′ ) = σ 2 1 r 1 (xx ′ ; θ 1 ). With this consistent set of hypotheses, the joint process (Z 1 (x), . . . , Z t (x)) x∈Q,t=1,...,s given σ 2 = (σ 2 i ) i=1,...,t , θ = (θ i ) i=1,...,t , β = (β i ) i=1,...,t and β ρ = (β ρ i-1 ) i=2,...,t , is Gaussian with mean:

E[Z t (x)|σ 2 , θ, β, β ρ ] = h ′ t (x) T β (3)
h ′ t (x) T = t-1 i=1 ρ i (x) f T 1 (x), t-1 i=2 ρ i (x) f T 2 (x), . . . , ρ t-1 (x)f T t-1 (x), f T t (x) (4
) and covariance:

cov(Z t (x), Z t (x ′ )|σ 2 , θ, β, β ρ ) = t j=1 σ 2 j   t-1 i=j ρ 2 i (x)   r j (x -x ′ ; θ j ) (5)
For each level t = 2, . . . , s, the experimental design D t is assumed to be such that D t ⊆ D t-1 . Note that this assumption is not necessary but allows us to have closed form expressions for the parameter estimation formulas. Furthermore, we denote by R t (D k , D l ) the correlation matrix between points in

D k and D l , 1 ≤ k, l ≤ s. R t (D k , D l ) is a (n k × n l ) matrix with (i, j) entry given by: [R t (D k , D l )] i,j = r t (x (k) i -x (l) j ; θ t ) 1 ≤ i ≤ n k 1 ≤ j ≤ n l
We will use the notation:

R t (D k ) = R t (D k , D k ). [7] present the case where ∀t ∈ [2, s], ρ t-1 (x) = ρ t-1 is constant.
Here, we will consider the general model presented in equations [START_REF] Dewettinck | Modeling the steady-state thermodynamic operation point of top-spray fluidized bed processing[END_REF]. We will also propose a new approach to estimate the coefficients (β t , β ρ t-1 ) t=2,...,s based on a Bayesian approach, which allows us to get information about their uncertainties. In the following section, we describe the case of 2 levels of code where the scaling coefficient ρ is constant and then we will extend it for s levels in Section 7. The general case in which ρ depends on x is addressed in Appendix A.

Building a model with 2 levels of code

Let us assume that we have 2 levels of code z 2 (x) and z 1 (x). From the previous section we assume that:

Z 2 (x) = ρZ 1 (x) + δ(x), x ∈ Q (Z 1 (x)) x∈Q ⊥ (δ(x)) x∈Q (6)
The goal of this section is to build a surrogate model for Z 2 (x) given the observations Z = z with an uncertainty quantification. The strategy is the following one. In Subsection 4.1 we describe the statistical distribution of the output Z 2 (x 0 ) at a new point x 0 given the parameters (β 1 , β 2 , ρ), (σ 2 1 , σ 2 2 ) and (θ 1 , θ 2 ) and the observations z. In Subsection 4.2 we describe the Bayesian estimation of the parameters (β 1 , β 2 , ρ) and (σ 2 1 , σ 2 2 ) given the observations. As pointed out at the end of Subsection 4.2 the hyper-parameters (θ 1 , θ 2 ) are estimated using a concentrated restricted log-likelihood method.

Conditional distribution of the output

For a point x 0 ∈ Q we determine in this subsection the distribution of

[Z 2 (x 0 )|Z = z, (β 1 , β 2 , ρ), (σ 2 1 , σ 2 2 ), (θ 1 , θ 2 )].
Standard results for normal distributions give that:

[Z 2 (x 0 )|Z = z, (β 1 , β 2 , ρ), (σ 2 1 , σ 2 2 ), (θ 1 , θ 2 )] ∼ N (m Z 2 (x 0 ), s 2 Z 2 (x 0 )) (7)
with mean function:

m Z 2 (x) = h ′ (x) T β + t(x) T V -1 (z -Hβ) (8)
and variance:

s 2 Z 2 (x) = ρ 2 σ 2 1 + σ 2 2 -t(x) T V -1 t(x) (9)
where we have denoted

β = β 1 β 2 , z = z 1 z 2
and where H is defined by:

H =            f T 1 (x (1) 1 ) 0 . . . . . . f T 1 (x (1)
n 1 ) 0 ρf T 1 (x (2) 1 ) f T 2 (x (2) 1
) . . . . . . ρf T 1 (x (2)
n 2 ) f T 2 (x (2)
n 2 )            =         F 1 (D 1 ) 0 ρF 1 (D 2 ) F 2 (D 2 )         with the notation F i (D j ) =    f T i (x (j) n 1 ) . . . f T i (x (j) n j )   . Furthermore, we have h ′ (x) = ρf T 1 (x), f T 2 (x) and: t(x) T = Cov(Z 2 (x), Z) = ρσ 2 1 R 1 ({x}, D 1 ), ρ 2 σ 2 1 R 1 ({x}, D 2 ) + σ 2 2 R 2 ({x}, D 2 ) (10)
The covariance matrix V of the Gaussian vector Z = Z 1 Z 2 can be written :

V = σ 2 1 R 1 (D 1 ) ρσ 2 1 R 1 (D 1 , D 2 ) ρσ 2 1 R 1 (D 2 , D 1 ) ρ 2 σ 2 1 R 1 (D 2 ) + σ 2 2 R 2 (D 2 ) (11)

Bayesian estimation of the parameters with 2 levels of code

In this Subsection, we describe the estimation of the parameters (β 1 , β 2 , ρ, σ 2 1 , σ 2 2 , θ 1 , θ 2 ) for the 2-level model given the observations Z = z. Due to the conditional independence between Z 1 (x) and δ(x), it is possible to estimate separately the parameters (β 1 , σ 2 1 , θ 1 ) and (β 2 , ρ, σ 2 2 , θ 2 ). We first describe the estimations of (β 1 , σ 2 1 ) given θ 1 and (β 2 , σ 2 2 , ρ) given θ 2 , which can be obtained in closed forms. We then describe how to estimate θ 1 and θ 2 .

Firstly, we consider the parameters (β 1 , σ 2 1 , θ 1 ). We choose the following non-informative prior distributions corresponding to the "Jeffreys priors" [START_REF] Jeffreys | Theory of Probability[END_REF]:

p(β 1 |σ 2 1 , θ 1 ) ∝ 1 p(σ 2 1 ) ∝ 1 σ 2 1 ( 12
)
Considering the probability density function of [Z 1 |β 1 , σ 2 1 , θ 1 ] and the Bayes formula, the posterior distribution of [β 1 |z 1 , σ 2 1 , θ 1 ] is :

[β 1 |z 1 , σ 2 1 , θ 1 ] ∼ N p 1 [F T 1 R 1 (D 1 ) -1 F 1 ] -1 [F T 1 R 1 (D 1 ) -1 z 1 ], [F T 1 R 1 (D 1 ) -1 σ 2 1 F 1 ] -1
(13) Then, using the Bayes formula, we obtain that the posterior distribution of

[σ 2 1 |z 1 , θ 1 ] is: [σ 2 1 |z 1 , θ 1 ] ∼ IG(α σ 2 1 |n 1 , Q 1 2 ) (14)
where IG(α, Q) stands for the inverse gamma distribution with density function:

p α,Q (x) = Q α Γ(α) e -Q x x α+1 1 x>0
and the parameters are given by:

α σ 2 1 |n 1 = n 1 -p 1 2 Q 1 = (z 1 -F 1 β1 ) T R 1 (D 1 ) -1 (z 1 -F 1 β1 ) (15)
with β1 = E[β 1 |z 1 , σ 2 1 , θ 1 ] = [F T 1 R 1 (D 1 ) -1 F 1 ] -1 [F T 1 R 1 (D 1 ) -1 z 1 ]
. Bayesian estimation of parameters with non-informative "Jeffreys priors" [START_REF] Jeffreys | Theory of Probability[END_REF] gives the same results as maximum likelihood estimation for the parameter β 1 . For the parameter σ 2 1 , the estimation given by

Q 1 2α σ 2 1 |n 1
is identical to the one obtained with the restricted maximum likelihood method. This method was introduced by Patterson and Thompson [START_REF] Patterson | Recovery of interblock information when block sizes are unequal[END_REF] in order to reduce the bias of the maximum likelihood estimator. Secondly, let us consider the set of parameters (β 2 , ρ, σ 2 2 , θ 2 ). In order to have closed form formulas for the estimation of (β 2 , ρ), we estimate them together. The idea to carry out a joint estimation is proposed for the first time in this paper and we believe it is important. Indeed, if the cheaper code is perfectly known, it can be considered as a regression function and so ρ will be a regression parameter. In this case, it is clear that a separated estimation of β 2 and ρ cannot be optimal. Using similar Jeffrey prior distributions as in [START_REF] Qian | Bayesian Hierarchical Modeling for Integrating Low-accuracy and High-accuracy Experiments[END_REF] and the same methodology as for the estimation of (β 1 , σ2 1 ), we find that:

[(ρ, β 2 )|z 1 , z 2 , σ 2 2 , θ 2 ] ∼ N p 2 +1 [F T R 2 (D 2 ) -1 F ] -1 [F T R 2 (D 2 ) -1 F ], [F T R 2 (D 2 ) -1 σ 2 2 F ] -1 (16) and: [σ 2 2 |z 2 , z 1 , θ 2 ] ∼ IG(α σ 2 2 |n 2 , Q 2 2 ) (17)
where:

α σ 2 2 |n 2 = n 2 -p 2 -1 2 Q 2 = (z 2 -F λ) T R 2 (D 2 ) -1 (z 2 -F λ) (18) with λ = E[(ρ, β 2 )|z 1 , z 2 , σ 2 2 , θ 2 ] = [F T R 2 (D 2 ) -1 F ] -1 [F T R 2 (D 2 ) -1 z 2 ]. The design matrix F is such that F = [ρz 1 (D 2 ) F 2 ]. Furthermore, the estima- tion of σ 2 2 given by Q 2 2α σ 2 2 |n 2
is the same as the restricted maximum likelihood one.

The hyper-parameters θ 1 and θ 2 are found by minimizing the opposite of the concentrated restricted log-likelihoods:

log (|det (R 1 (D 1 )) |) + (n 1 -p 1 )log( σ1 2 ) (19) log (|det (R 2 (D 2 )) |) + (n 2 -p 2 -1)log( σ2 2 ) (20)
These minimizations problems must be numerically solved with a global optimization method. We use an evolutionary method coupled with a BFGS algorithm. The drawback of the maximum likelihood estimation (see [START_REF] Lehmann | Theory of Point Estimation[END_REF]) is that, contrarily to Bayes estimation, we do not have any information about the variance of the estimator. Nevertheless, Bayes estimation of the hyper parameters θ 1 and θ 2 are prohibitive and as noted in [START_REF] Santner | The Design and Analysis of Computer Experiments[END_REF] the choice of the prior distribution is non trivial. Therefore, in this paper, we will always estimate these parameters with a concentrated restricted likelihood method.

5 Bayesian prediction for a code with 2 levels

The aim of a Bayesian prediction is to provide a predictive distribution for Z s (x) integrating the posterior distributions of the parameters and hence taking into account their uncertainty.

A Bayesian prediction for a code with s = 2 levels was suggested in [START_REF] Qian | Bayesian Hierarchical Modeling for Integrating Low-accuracy and High-accuracy Experiments[END_REF]. Nevertheless, we propose here a new Bayesian approach with some significant differences. First, we assume that the adjustment coefficient is a regression function whereas Qian and Wu [START_REF] Qian | Bayesian Hierarchical Modeling for Integrating Low-accuracy and High-accuracy Experiments[END_REF] model it with a Gaussian process. Secondly, we use different prior distributions for the parameter estimation. More specifically, according to the Bayesian estimation of parameters previously presented, we use a joint prior distribution for (β 2 , ρ) conditioned by σ 2 by σ 2 2 . Then, we use a hierarchy between the different parameters. At the lowest level is the regressor parameter β. At the second level is the variance parameter σ 2 which controls the distribution of the parameter β. At the top level is the parameter θ which controls the distribution of the parameters at the bottom levels. It is common to use a hierarchical specification of models for Bayesian prediction as presented in [START_REF] Rasmussen | Gaussian Processes for Machine Learning[END_REF]. This strategy will allow us to obtain explicit formulas for the joint distribution of the parameters and above all, to reduce dramatically the cost of the numerical implementation of the complete Bayesian prediction.

We will also present the case in which we do not have any prior information about the parameters. As described in the previous section, the hyper parameter θ is estimated by minimizing the opposite of the concentrated restricted log-likelihood and it is assumed to be fixed to this estimated value from now on.

Prior distributions and Bayesian estimation of the parameters

Many choices of priors can be made for the Bayesian modelling. Here we study the two following cases:

(I) Priors for each parameter are informative.

(II) Priors for each parameter are non-informative.

For the non-informative case (II), we use the improper distributions corresponding to the "Jeffreys prior" and then the posterior distributions are given in Section 4.2. Note that non-informative distributions are used when we do not have prior knowledge. For the informative case (I), we will consider the following prior distributions:

[β 1 |σ 2 1 ] ∼ N p 1 (b 1 , σ 2 1 V 1 ), [(ρ, β 2 )|z 1 , σ 2 2 ] ∼ N 1+p 2 b λ = b ρ b 2 , σ 2 2 V λ = σ 2 2 V ρ 0 0 V 2 [σ 2 1 ] ∼ IG(α 1 , γ 1 ), [σ 2 2 |z 1 ] ∼ IG(α 2 , γ 2 ) where b 1 ∈ R p 1 , b λ ∈ R 1+p 2 , V 1 is a (p 1 × p 1 ) diagonal matrix, V λ is a ((1 + p 2 ) × (1 + p 2 )) diagonal matrix and α 1 , γ 1 , α 2 , γ 2 > 0.
The forms of the priors are chosen in order to be able to get closed form expressions for the posterior distributions. Note that there are enough free parameters in the priors to allow the user to prescribe their means and variances. From the previous prior definitions, the posterior distributions of the parameters are:

[β 1 |z 1 , σ 2 1 ] ∼ N p 1 (A 1 i ν 1 i , A 1 i ) [(ρ, β 2 )|z 1 , z 2 , σ 2 2 ] ∼ N p 2 +1 (A λ i ν λ i , A λ i ) (21)
where:

A 1 i =    [F T 1 R -1 1 (D 1 ) σ 2 1 F 1 + V -1 1 σ 2 1 ] -1 i = (I) [F T 1 R -1 1 (D 1 ) σ 2 1 F 1 ] -1 i = (II) ν 1 i =    [F T 1 R -1 1 (D 1 ) σ 2 1 z 1 + V -1 1 σ 2 1 b 1 ] i = (I) [F T 1 R -1 1 (D 1 ) σ 2 1 z 1 ] i = (II) (22) A λ i =    [F T R -1 2 (D 2 ) σ 2 2 F + V -1 λ σ 2 2 ] -1 i = (I) [F T R -1 2 (D 2 ) σ 2 2 F ] -1 i = (II) ν λ i =    [F T R -1 2 (D 2 ) σ 2 2 z 2 + V -1 λ σ 2 2 b λ ] i = (I) [F T R -1 2 (D 2 ) σ 2 2 z 2 ] i = (II) (23) and F = [ρz 1 (D 2 ) F 2 ]
. Furthermore, we have:

[σ 2 1 |z 1 ] ∼ IG(α σ 2 1 |n 1 i , Q 1 i 2 ), [σ 2 2 |z 2 , z 1 ] ∼ IG(α σ 2 2 |n 2 i , Q 2 i 2 ) (24)
where:

Q 1 i = γ 1 + (b 1 -β1 ) T (V 1 + [F T 1 R -1 1 (D 1 )F 1 ] -1 ) -1 (b 1 -β1 ) + Q 1 2 i = (I) z T 1 [R -1 1 (D 1 ) -R -1 1 (D 1 )F 1 (F T 1 R -1 1 (D 1 )F 1 ) -1 F T 1 R -1 1 (D 1 )]z 1 i = (II) Q 2 i = γ 2 + (b λ -λ) T (V λ + [F T R -1 2 (D 2 )F ] -1 ) -1 (b λ -λ) + Q 2 2 i = (I) z T 2 [R -1 2 (D 2 ) -R -1 2 (D 2 )F (F T R -1 2 (D 2 )F ) -1 F T R -1 2 (D 2 )]z 2 i = (II) β1 = (F T 1 R -1 1 (D 1 )F 1 ) -1 F T 1 R -1 1 (D 1 )z 1 λ = (F T R -1 2 (D 2 )F ) -1 F T R -1 2 (D 2 )z 2 α σ 2 1 |n 1 i = n 1 2 + α 1 i = (I) n 1 -p 1 2 i = (II) α σ 2 2 |n 2 i = n 2 2 + α 2 i = (I) n 2 -p 2 -1 2 i = (II)
Mixing of informative and non-informative priors are of course possible and easy to implement. As we will discuss in Subsection 5.4 and see in the examples of Section 6, the use of informative priors has minor impact on the mean estimation but may have a strong impact on variance estimation.

Predictive distributions when

β 2 , ρ, σ 2
1 and σ 2 2 are known

As a preliminary step towards the Bayesian prediction carried out in the next subsection, we give here Bayesian prediction in the form of closed form expressions when the parameters

β 2 , ρ, σ 2 1 and σ 2 2 are known. The conditional distribution of [Z 2 (x)|Z = z, β 2 , ρ, σ 2 1 , σ 2 2
] is given by:

[Z 2 (x)|Z = z, β 2 , ρ, σ 2 1 , σ 2 2 ] ∼ N µ i (x), σ 2 i (x) (25)
where:

µ i (x) = h ′ (x) T A 1 i ν 1 i β 2 + t(x) T V -1 z -H A 1 i ν 1 i β 2 σ 2 i (x) = s 2 Z 2 (x) + k 1 A 1 i k T 1
and A 1 i and ν 1 i are defined by ( 22). Note that the estimated variance is augmented by the term k 1 A 1 i k T 1 which quantifies the uncertainty due to the estimation of

β 1 . k 1 is a (1 × p 1 ) vector composed of the p 1 first elements of the (1 × p 1 , 1 × p 2 ) vector k = (k 1 , k 2 ) = h ′ (x) T -t(x) T V -1
H. H is given by (4.1). The existence of closed form formulas is important as it will allow for a fast numerical implementation.

Bayesian prediction

Before performing the Bayesian prediction we note that -thanks to the explicit joint prior distribution for β 2 and ρ, the independence hypotheses and the hierarchical specification of the parameters -conditioning on θ, we have an explicit formula for the following joint density:

p(β 1 , β 2 , ρ, σ 2 1 , σ 2 2 |z 1 , z 2 ) = p(β 1 |σ 2 1 , z 1 )p(β 2 , ρ|σ 2 2 , z 1 , z 2 )p(σ 2 1 |z 1 )p(σ 2 2 |z 1 , z 2 )
(26) This explicit joint density is an original result which contrasts with [START_REF] Qian | Bayesian Hierarchical Modeling for Integrating Low-accuracy and High-accuracy Experiments[END_REF] and which allows us to avoid prohibitive implementation for the Bayesian analysis.

First, we consider the predictive distribution with σ 2 1 and σ 2 2 known. Considering the conditional independence assumption between (δ(x)) x∈Q and

(Z 1 (x)) x∈Q , the probability density function of [Z 2 (x)|Z = z, σ 2 1 , σ 2 2
] can be deduced from the following integral:

p(z 2 (x)|z 1 , z 2 , σ 2 1 , σ 2 2 ) = R 1+p 2 p(z 2 (x)|z 1 , z 2 , β 2 , ρ, σ 2 1 , σ 2 2 )p(ρ, β 2 |z 1 , z 2 , σ 2 2 ) dρdβ 2 (27) where p(z 2 (x)|z 1 , z 2 , β 2 , ρ, σ 2 1 , σ 2
2 ) is given by [START_REF] Kennedy | Predicting the output from a complex computer code when fast approximations are available[END_REF]. This integral has to be numerically evaluated. Since [ρ, β 2 |z 1 , z 2 , σ 2 2 ] has a known normal distribution given by (21), we here use a crude Monte-Carlo algorithm when the dimension of β 2 and ρ is high, or a trapezoidal quadrature method when it is low.

Then, we infer from the parameters σ 2 1 and σ 2 2 . Due to the independence between (δ(x)) x∈Q and (Z 1 (x)) x∈Q , the probability density function of

[Z 2 (x)|Z = z] is: p(z 2 (x)|z 1 , z 2 ) = R 2 p(z 2 (x)|z 1 , z 2 , σ 2 1 , σ 2 2 )p(σ 2 1 |z 1 )p(σ 2 2 |z 1 , z 2 ) dσ 2 1 dσ 2 2 (28)
where p(σ 2 1 |z 1 ) and p(σ 2 2 |z 1 , z 2 ) are given by (24). This integral has also to be numerically evaluated. Since we have a double integration, a quadrature method will be efficient. We use here a trapezoidal numerical integration, defining the region of integration [σ ). We hence have a predictive distribution taking into account the uncertainties due to the parameter estimations.

Discussion about the numerical evaluations of the integrals

We saw in the previous section that we can obtain an analytical prediction when β 2 , ρ, σ 2 1 and σ 2 2 are known. From this analytical formula, we can have a Bayesian prediction with only two nested integrations. One of them can be approximated with a quadrature or a crude Monte Carlo method, which is not too expensive. The other is a double integration approximated with a quadrature method which is efficient and not expensive. Therefore, we do not use any Markov chain Monte Carlo method and we considerably reduce the time and the complexity of the method. This allows us to easily build an accurate Bayesian metamodel. Practically, we use 441 integration points to approximate (28) and 1000 Monte-Carlo particles to approximate (27). Therefore, we have 441000 call to the predictive density function (25).

To avoid a prohibitive implementation, another approach has also been proposed in [START_REF] Cumming | Small Sample Bayesian Designs for Complex High-Dimensional Models Based on Information Gained Using Fast Approximations[END_REF]. They adopt a Bayes linear formulation which requires only the specification of the means, variances, and covariance. See [START_REF] Goldstein | Bayes Linear Statistics: Theory and Methods[END_REF] for further details about the Bayes linear approach. The strength of this method is that its computational cost is low. Nonetheless, since it only focuses on posterior means and covariances, it does not provide the full posterior predictive distribution.

Finally, we highlight the fact that our Bayesian procedure can be used to perform multi-fidelity analysis with more than 2 levels of code whereas the cost of the one presented in [START_REF] Qian | Bayesian Hierarchical Modeling for Integrating Low-accuracy and High-accuracy Experiments[END_REF] is too high to allow for such analysis. We illustrate in Section 8 through an industrial case the great practical importance of using more than 2 levels of code.

Toy examples

We will present in this section some co-kriging metamodels using one-dimensional functions inspired by the example presented in [START_REF] Forrester | Multi-fidelity optimization via surrogate modelling[END_REF]. For the following examples, we will use a non-Bayesian co-kriging modeli.e. the one presented in [START_REF] Kennedy | Predicting the output from a complex computer code when fast approximations are available[END_REF] -but with a Bayesian estimation of the parameters (see Section 4.2) and for the second example we will use a Bayesian co-kriging. Furthermore, the correlation kernels are assumed to be:

r t (x (k) i -x (l) j ; θ t ) = exp - x (k) i -x (l) j 2 θ 2 t where t, k, l = 1, 2 1 ≤ i ≤ n 1 1 ≤ j ≤ n 2 .
Example 1. The aim of this example is to emphasize the effectiveness of the presented Bayesian estimation of the parameters (see Section 4.2). We assume that the expensive code is given by z 2 (x) = (6x-2) 2 sin(12x-4) and the cheaper code by z 1 (x) = 0.5z 2 (x) + 10(x -0.5) -5. The experimental design set of the cheapest code is D 1 = {0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1} and the one of the expensive code is D 2 = {0, 0.4, 0.6, 1}. This example is identical to the one-dimensional demonstration presented in [START_REF] Forrester | Multi-fidelity optimization via surrogate modelling[END_REF]. Figure 1 shows the functions x → z 2 (x) and x → z 1 (x), the training data for z 2 and z 1 , the ordinary kriging using only the expensive data and the co-kriging using expensive and cheap data. To validate the model, the Root-mean-square errors (RMSE) and

Q 2 = 1 -x∈T (mZ 2 (x)-z 2 (x)) 2 x∈T (mZ 2 (x)-z 2) 2 are computed.
The test set T is composed of a regular grid points sampled from 0 to 1 with a grid step equal to 0.01 and z2 is the empirical mean evaluated in T . The estimated RMSE is 5.68 × 10 -2 and the coefficient Q 2 is 99.98%, so we have a prediction error closed to 0. The Bayesian estimation of the parameters of co-kriging are given in Table 1. Furthermore, the estimations of the hyper-parameters (θ 1 , θ 2 ), calculated by maximizing the concentrated log-likelihoods ( 19) and (20), are θ1 = 0.25 and θ2 = 0.80. D 1 being a regular grid with a grid step equal to 0.1 and D 2 being composed of points sampled from 0 to 1, points of the experimental designs are hence strongly correlated which will imply a smooth surrogate model. We see that the Bayesian estimation of parameters is very effective since the estimations of parameters ρ and β 2 are perfect. Nevertheless this example does not highlight the strength of the method since there is a relation between z 2 (x) x∈[0,1] and z 1 (x) x∈[0,1] which exactly corresponds to the equation ( 2) with the error δ 2 that can be written in terms of the regression functions f 2 exactly. Therefore, if the cheap code is well modelled, like in our case, the co-kriging is equivalent to a linear regression. Moreover, the very small value of σ 2 2 illustrates this. Example 2. This example illustrates a case where the non-Bayesian cokriging underestimates the predictive variance whereas the Bayesian one adjusts it. We assume that the expensive code is given by z 2 (x) = (6x -2) 2 sin(12x -4) + sin(10 cos(5x)) and the cheaper code is given by z 1 (x) = 0.5((6x-2) 2 sin(12x-4))+10(x-0.5)-5. Through the term sin(10 cos(5x)), the expensive code has high frequencies which are not captured by the cheap code and the error δ 2 is not a simple linear combination of the regression functions f 2 . Figure 2 shows the results of kriging and co-kriging for these two functions. The estimated RMSE is 1.05 and the coefficient Q 2 is 93.57%, we still have a good prediction. The Bayesian estimations of the parameters are given in Table 2 and we have θ1 = 0.25 and θ2 = 0.07. The values of these parameters have been fixed according the following arguments. As the cheap code is the same as the one of the Example 1, we keep the same estimation for θ 1 . Then, we consider that there are not enough points to carry out a significant estimation of θ 2 . Therefore, we fix the value of θ2 according to the high frequencies introduced by the term sin(10 cos(5x)). Due to the additional term sin(10 cos(5x)), the estimation of the parameter ρ is less effective than in the first example. This highlights the dependence between the estimation of ρ and the mean of δ(x) x∈[0,1] . Furthermore, Figure 3 represents the confidence interval at plus or minus twice the standard deviation of the predictive distribution in the Bayesian and non-Bayesian case. We see that we underestimate the variance of the predictive distribution in the non-Bayesian case. This estimation is well adjusted in the Bayesian case.

The case of s levels of code

The aim of this Section is to perform a multi-level co-kriging with any number of codes. Let us consider s levels of code. The generalization of the previous model is straightforward. Actually, if we note β = (β T 1 , . . . , β T s ) T , ρ = (ρ 1 , . . . , ρ s-1 ), σ 2 = (σ 2 1 , . . . , σ 2 s ) and θ = (θ 1 , . . . , θ s ), we have:

∀x ∈ Q [Z s (x)|Z = z, β, ρ, σ 2 , θ] ∼ N m Zs (x), s 2 Zs (x)
where:

m Zs (x) = h ′ s (x) T β + t s (x) T V -1 s (z -H s β) (29)
and:

s 2 Zs (x) = σ 2 Zs -t s (x) T V -1 s t s (x) (30)
Furthermore, let us denote by R t = R t (D t ) the correlation matrix for D t and ρ s = 0, ∀s ≤ 0. The matrix V s has the form:

V s =   
V (1,1) . . . V (1,s) . . . . . . . . .

V (s,1) . . . V (s,s)    (31)
The s diagonal blocks of size n t × n t are defined by:

V (t,t) = σ 2 t R t (D t ) + σ 2 t-1 ρ 2 t-1 R t-1 (D t ) + • • • + σ 2 1 t-1 i=1 ρ 2 i R 1 (D t ) (32)
and the off-diagonal blocks of size n t × n t ′ are given by:

V (t,t ′ ) = t ′ -1 i=t ρ i V (t,t) (D t , D t ′ ) 1 ≤ t < t ′ ≤ s (33)
The vector t s (x) is such that t s (x) = (t * 1 (x, D 1 ) T , . . . , t * s (x, D s ) T ) T , where:

t * t (x, D t ) T = ρ t-1 t * t-1 (x, D t ) T + s-1 i=t ρ i σ 2 t R t (x, D t ) 1 < t ≤ s (34)
where

s-1 i=s ρ i = 1 and t * 1 (x, D 1 ) T = s-1 i=1 ρ i σ 2 1 R 1 (x, D 1
). If we define:

F k (D l ) =    f T k (x (l) 1 ) . . . f T k (x (l) n l )    1 ≤ k, l ≤ s
then the matrix H s can be written as:

H s =         F 1 (D 1 ) ρ 1 F 1 (D 2 ) F 2 (D 2 ) 0 ρ 1 ρ 2 F 1 (D 3 ) ρ 2 F 2 (D 3 ) . . . . . . . . . s-1 i=1 ρ i F 1 (D s ) s-1 i=2 ρ i F 2 (D s ) . . . F s (D s )         (35)
h ′ s (x) T and var(Z s (x)) = σ 2 Zs are given by the equations ( 3) and (5).

Bayesian estimation of parameters for s levels of code

From the assumptions of conditional independence between (δ t (x)) x∈Q and (Z t-1 (x), . . . , Z 1 (x)) x∈Q , we can extend the Bayesian estimation of the parameters to the case of s levels. Note that we do not assume the independence of β t and ρ t-1 . We can obtain a closed form expression for the estimation of (β t , ρ t-1 ). For all t = 2, . . . , s, we have:

[(ρ t-1 , β t )|z t , z t-1 , θ t , σ 2 t ] ∼ N H T t R -1 t (D t )H t -1 H T t R -1 t (D t )z t , σ 2 t H T t R -1 t (D t )H t -1 (36) where H t = [ρ t-1 z t-1 (D t ) F t (D t )]. Furthermore, if we note λt = E[(ρ t-1 , β t )|z t , z t-1 , θ t , σ 2 t ], then we have: [σ 2 t |z t , z t-1 , θ t ] ∼ IG(α t , Q t 2 ) (37)
where

α t = nt-pt-1 2 and Q t = (z t -H t λt ) T R -1 t (D t )(z t -H t λt ). The REML estimator of σ 2
t is σ2 t = Qt 2αt and we can estimate θ t by minimizing the expression:

log(|det(R t (D t ))|) + (n t -p t -q t-1 )log(σ 2 t ) (38)
The generalization of the Bayesian estimation previously presented is important since it shows that the parameters estimation for a s-levels co-kriging is equivalent to the one for s independent krigings.

Reduction of computational complexity of inverting the covariance matrix V s

V s is an ( s i=1 n i × s i=1 n i ) matrix, its inverse can hence be difficult to process. We present in this Subsection a method to reduce the complexity of the processing of V -1 s . By sorting the experimental design sets such that ∀t = 2, . . . , s, D t-1 = (x

(t-1) 1 , . . . , x (t-1) n t-1 -nt , x (t) 1 , . . . , x (t) nt ) = (D t-1 \ D t , D t )
it can be shown that ∀t = 2, . . . , s the inverse of the matrix V s has the form:

V -1 s =     V -1 s-1 + 0 0 0 ρ 2 s-1 R -1 s σ 2 s - 0 ρ s-1 R -1 s σ 2 s -0 ρ s-1 R -1 s σ 2 s R -1 s σ 2 s     V -1 1 = R -1 1 σ 2 1 (39) with V -1 s-1 an ( s-1 i=1 n i × s-1 i=1 n i ) matrix and R -1
s an (n s × n s ) matrix. This is a very important result since it shows that we can deduce V -1 s from R -1 t , t = 1, . . . , s. Therefore, the complexity of the processing of V -1

s is O( s i=1 n 3 i ) instead of O(( s i=1 n i ) 3
). Furthermore, from the equation (39) and the Bayesian estimation of parameters presented in Section 7.1, we have shown here that building a s-level co-kriging is equivalent to build s independent krigings. We emphasize that, for practical applications, the form (39) for the inverse of V s allows us to perform fine matrix regularization in the case of ill-conditioned problems. Indeed, V s is invertible if and only if the matrices R t , t = 1, . . . , s are invertible. Therefore, if the problem is ill-conditioned, we just have to regularize the matrices R t which are illconditioned too. Then, since (t * 1 (x, D 1 ) T , . . . , t * s-1 (x, D s-1 ) T ) = ρ s-1 t T s-1 (x) it can also be shown that in the equation (29):

t s (x) T V -1 s = ρ s-1 t T s-1 (x)V -1 s-1 -[0 1×( s-1 i=1 n i -ns) , ρ s-1 R s ({x}, D s )R -1 s ],R s ({x}, D s )R -1 s (40) Therefore, t s (x) T V -1 s is independent of σ 2 s . Since t 1 (x) T V -1 1 = R 1 ({x}, D 1 )R -1 1
does not depend on σ 2 1 , by induction, t s (x) T V -1 s is independent of σ 2 i for all 1 ≤ i ≤ s. We have just shown here that the co-kriging mean does not depend on the variance coefficients.

Numerical test on the reduction of computational complexity

In the previous section, we have presented a reduction of complexity for the co-kriging model by expressing the inverse of the matrix V s with the inverses of the matrices R t , t = 1, . . . , s. We present here a numerical test to highlight the gain of CPU time obtained with this method. We focus on the case of 2 levels of code and we consider the Gaussian kernel for the 2 levels:

r(x -x ′ ; θ) = exp - ||x -x ′ || 2 θ 2
The experimental design set for the cheap code is a regular grid composed of n 1 points between 0 and 1 and the experimental design set for the expensive code are the n 2 first points of this grid. We consider the relation n 1 = 4n 2 with n 2 = 50, 60, . . . , 500 and the parameter θ = 5 n 2 (the parameter θ is controlled by n 2 in order to avoid ill-conditioned covariance matrices). The total number of observations is hence n = n 1 + n 2 . Figure 4 compares the CPU time needed to build a co-kriging model with or without reduction complexity.

First, the slope of the two CPU times is close to 3 (the least-squares estimation value is 3.03). The complexity of a matrix inversion being O(n 3 ), with n the size of the matrix, the estimation of the slope highlights the fact that it is the matrix inversion which leads the CPU time. Then, Figure 4 emphasizes that the reduction of complexity is worthwhile. Indeed, we see that the ratio between the two CPU time is approximately a constant equal to 1.93. We are hence close to the theoretical ratio equal to (n 1 + n 2 ) 3 /(n 3 1 + n 3 2 ) ≈ 1.92 which is obtained when we consider that the CPU time is essentially due to the matrix inversion.

Toy example on the complexity reduction

A 3-level co-kriging metamodel is presented in this section to illustrate the gain of CPU which can be obtained with the presented reduction of complexity. We focus on the inversion of the co-kriging matrix V s by comparing the CPU time needed with a direct inversion or by using the formula (39). We assume that the 3 levels of code are given by the followings three dimensional functions: popular in the field of sensitivity analysis [START_REF] Saltelli | Sensitivity analysis, Wiley Series in Probability and Statistics[END_REF]. We consider n 3 = 50 observations for the most accurate code z 3 (x), n 2 = 200 for the intermediate code and z 1 = 400 for the less accurate code. All experimental design sets are randomly sampled from the uniform distribution. As presented in Section 3 we consider nested experimental designs ∀t = 2, . . . , s D t ⊆ D t-1 .

z 1 (x) = sin(x 1 ) (41) z 2 (x) = z 1 (x) + asin(x 2 ) 2 (42) z 3 (x) = z 2 (x) + bx 4 3 sin(x 1 ) (43) with x = (x 1 , x 2 , x 3 ) ∈ [-π, π]
We use a tensorised Matérn-5 2 kernel for the three correlation functions:

r t (x, x ′ ; θ t ) = d i=1 r 1D (x i , x ′ i ; θ t,i ) (44) with r 1D (t, t ′ ; θ) = 1 + √ 5 |t-t ′ | θ + 5 3 (t-t ′ ) 2 θ 2 exp - √ 5 |t-t ′ | θ , t, t ′ ∈ R.
The estimations of the hyper-parameters θ t are presented in Table 3.

Parameter Estimation

θ1

( 0.61 1.99 2.04 ) θ2

( 1.98 0.26 2.48 ) θ3 ( 0.23 0.89 0.21 ) Table 3: Toy example on the complexity reduction. Estimation of the hyperparameters (correlation lengths) for the 3-level co-kriging.

The hyper-parameter estimates show us that z 1 (x) is very smooth in the directions x 2 and x 3 reflecting the fact that it depends only on the first direction x 1 . Similarly, the bias between z 2 (x) and z 1 (x) only depending on the second direction x 2 , it is rougher on this direction and very smooth in the others. Finally, the bias between z 3 (x) and z 2 (x) is rougher in the direction x 3 than in the directions x 1 and x 2 . This is due to the important impact of x 3 on the third level.

The estimation of the variance, scale and regression parameters are given in Table 4.

Parameter Estimation

β 1 0.00 ρ 1 β 2 0.99 2.44 ρ 2 β 3 0.95 0.64 σ 2 1 0.09 σ 2 2
1.66 σ 2 3 6.25 Table 4: Toy example on the complexity reduction. Estimation of the variance, scale and regression parameters for the 3-level co-kriging.

Table 4 shows the efficiency of the suggested method for the parameter estimations since it provides very accurate estimations of ρ 1 and ρ 2 .

To evaluate the accuracy of the co-kriging model, we use a test set of 30,000 points uniformly sampled from the uniform distribution. Then, we compute the coefficient Q 2 with the co-kriging predictions and the responses of z 3 (x) on this set. We obtain for the co-kriging model Q 2 = 83.21%, we hence have a middling accuracy despite the large number of observations used. Nonetheless, we have a significant improvement relatively to the kriging model since with the same kernel and the same experimental design set D 3 we obtain Q 2 = 47.97% which is a very poor accuracy. The hyper-parameter estimation of the kriging model is θ = (0.79, 0.14, 0.29), the variance one is σ 2 = 13.66 and the trend coefficient one is β = 3.89.

Let us now compare the difference of CPU time between the co-kriging building with a crude inversion of the covariance matrix V s and the one with an inversion using the formula presented in Subsection 7.2. The CPU time necessary without the reduction complexity is CPU crude = 0.47 whereas the one necessary with the complexity reduction is CPU light = 0.14. We hence find that the CPU time ratio between the two methods approximately equals 3.36. This is not far from the theoretical ratio which equals 650 3 /(400 3 + 200 3 + 50 3 ) ≈ 3.80. We note that the complexity reduction could be of important practical interest. For example, without it the computational cost of a leave-one-out cross validation procedure will be much more important (the ratio will still be around 3 in our example). The complexity of this procedure being O(n 4 ), the gain of CPU time will be substantially.

Example : Fluidized-Bed Process

This example illustrates the comparison between 2-level and 3-level co-kriging. A 3-level co-kriging method is applied to a physical experiment modelled by a computer code. The experiment, which is the measurement of the temperature of the steady-state thermodynamic operation point for a fluidized-bed process, was presented by [START_REF] Dewettinck | Modeling the steady-state thermodynamic operation point of top-spray fluidized bed processing[END_REF], who developed a computer model named "Topsim" to calculate the measured temperature. The code, developed for a Glatt GPCG-1, fluidized-bed unit in the top-spray configuration, can be run at 3 levels of complexity. We hence have 4 available responses:

1. T exp : the experimental response.

2. T 3 : the most accurate code modelling the experiment.

3. T 2 : a simplified version of T 3 . 4. T 1 : the lowest accurate code modelling the experiment.

The differences between T 1 , T 2 and T 3 are discussed by [START_REF] Dewettinck | Modeling the steady-state thermodynamic operation point of top-spray fluidized bed processing[END_REF]. The aim of this study is to predict the experimental response T exp given the two levels of code T 3 and T 2 . We only focus on a 3-level co-kriging using T 3 and T 2 to predict T exp since 28 responses available for each level is not enough to build a nested experimental design relevant for a 4-level co-kriging. The experimental design set and the responses T 1 , T 2 , T 3 and T exp are given by [START_REF] Qian | Bayesian Hierarchical Modeling for Integrating Low-accuracy and High-accuracy Experiments[END_REF] who have presented a 2-level co-kriging using T exp and T 2 . Furthermore, the responses are parameterized by a 6-dimensional input vector presented by [START_REF] Dewettinck | Modeling the steady-state thermodynamic operation point of top-spray fluidized bed processing[END_REF].

Building the 3-level co-kriging

To build the 3-level co-kriging, we use 10 measures of T exp (measures 1, [START_REF] Fang | Uniform desgin: Theory and application[END_REF][START_REF] Jeffreys | Theory of Probability[END_REF][START_REF] Lehmann | Theory of Point Estimation[END_REF][START_REF] Qian | Bayesian Hierarchical Modeling for Integrating Low-accuracy and High-accuracy Experiments[END_REF][START_REF] Santner | The Design and Analysis of Computer Experiments[END_REF]18,19,20,27 in Table 4 in [START_REF] Qian | Bayesian Hierarchical Modeling for Integrating Low-accuracy and High-accuracy Experiments[END_REF]), 20 simulations of T 3 (runs 1, 2, [START_REF] Fang | Uniform desgin: Theory and application[END_REF][START_REF] Goldstein | Bayes Linear Statistics: Theory and Methods[END_REF][START_REF] Harville | Matrix Algebra from Statistician's Perspective[END_REF][START_REF] Kennedy | Predicting the output from a complex computer code when fast approximations are available[END_REF][START_REF] Jeffreys | Theory of Probability[END_REF][START_REF] Jones | Efficient Global Optimization of Expensive Black-Box Functions[END_REF][START_REF] Lehmann | Theory of Point Estimation[END_REF][START_REF] Patterson | Recovery of interblock information when block sizes are unequal[END_REF][START_REF] Qian | Bayesian Hierarchical Modeling for Integrating Low-accuracy and High-accuracy Experiments[END_REF][START_REF] Rasmussen | Gaussian Processes for Machine Learning[END_REF][START_REF] Santner | The Design and Analysis of Computer Experiments[END_REF][START_REF] Stocki | A method to improve design reliability using optimal Latin hypercube sampling[END_REF]18,19,20,22,24,27) and the 28 simulations of T 2 and the input vector is scaled between 0 and 1. The last 18 measures of T exp are used for validation. The design sets are nested such that D t-1 = (D t-1 \ D t , D t ) for t = 2, 3 and we use a Matern 5 2 kernel for the three covariance functions. The estimations of the hyper-parameters which represent correlation lengths of the three covariance kernels are given in Table 5. The estimations of hyper-parameters in Table 5 show us that the surrogate model will be very smooth in the first four directions. For the fifth direction the Gaussian processes modelling the cheap code T 2 and the bias between T exp and T 3 are very smooth and the one modelling the bias between T 3 and T 2 is close to a regression. Finally, the model is sharper in the sixth direction in particular for the two biases where correlation lengths are around 0.3.

θ1
Furthermore, Table 6 gives the estimation of the variance and regression parameters (see section 7.1). Table 6 shows that the responses have approximately the same scale since the adjustment coefficients are close to 1. Furthermore, we see an important bias between T 3 and T 2 with β 3 = 1.93. Finally, the variance coefficients for the biases indicate that they are possibly much simpler to model than the cheap code T 2 as their estimations are smaller.

Regression coefficient Posterior mean

8.2 3-level co-kriging prediction: predictions when code output is available

The aim of this Section is to show that co-kriging can improve significantly the accuracy of the surrogate model at points where at least one level of responses is available.

The predictions of the 3-level co-kriging are here presented and compared with the predictions obtained with a 2-level co-kriging using only the 10 responses of T exp and the 20 responses of T 3 . The predictions for the 2-level and the 3-level co-krigings vs. the real values (i.e., the measured temperature T exp ) are shown in Figure 5. The 3-level co-kriging gives us the same prediction means as the 2-level co-kriging at the 10 points (points 2, 5, 6, 7, 9, 11, 13, 16, 22, 24) where T 3 is known. These overlapped points mean that T 2 does not influence the surrogate model at these points. This follows from the Markov property introduced in Section 3, which implies that the prediction of T exp is entirely determined by T 3 at these points. We also note that, in general, the 2-level co-kriging predictions -at points where T 3 is unknown -are not accurate and the 3-level co-kriging improves significantly the prediction means compared to the 2-level co-kriging. Table 7 compares the 2-level co-kriging with the 3-level co-kriging and summarizes some results about the quality of the predictions on the 18 validation points. Nonetheless, it is important to notice that, in the 3-level case, the output of the cheapest code T 2 is known at the 18 test points. This means that the results of this subsection show that the 3-level co-kriging prediction is more accurate than the 2-level co-kriging prediction at a point where the cheapest response T 2 is available. In the next subsection we show that the 3-level co-kriging prediction is more accurate than the 2-level one at a point where no response is available. Figure 6 shows the prediction errors of the 2-level co-kriging and the confidence interval at plus or minus twice the prediction standard deviation. The last 10 prediction errors and their confidence intervals are the same as those of the 3-level case since it corresponds to the points where T 3 is known. We see in Figure 6 that the confidence intervals are well predicted. Furthermore, we see a significant difference between the accuracy of the prediction means and their confidence intervals for the point where T 3 is unknown (the 8 first validation points) and for the ones where it is known (the last 10 validation points).

Q 2 RMSE

3-level co-kriging prediction: predictions when code output is not available

In this subsection, we show that a multi-level co-kriging can significantly improve the prediction of a surrogate model at points where no response is available.

We have seen in Section 8.2 that the 3-level co-kriging improves significantly the 2-level co-kriging at points where T 3 is unknown and T 2 has been sampled. Nevertheless, to have a fair comparison between these two co-kriging models, we compare their accuracy by applying a Leave-One-Out Cross-Validation (LOO-CV) procedure at the 10 points where T exp is known. This means that we perform for each of these 10 points the following procedure:

1. The experimental and the two code outputs corresponding to the point are removed from the data set.

2. The 2-level co-kriging method and the 3-level co-kriging method are applied using the truncated data set in order to give a confidence interval for the experimental output at the point.

Figure 7 shows the result of the LOO-CV procedure for the 2-level and 3level co-kriging. We see that the 3-level co-kriging is more accurate than the 2-level one. Indeed, the LOO-CV RMSE for the 2-level co-kriging is equal to 1.88 whereas it is equal to 1.09 for the 3-level co-kriging. This shows that the 3-level co-kriging provides better predictions also at points where no response is available. This highlights the strength of the proposed method and shows that a co-kriging method with more than 2 levels of code can be worthwhile. of the 2-level and 3-level co-kriging. We see that the confidence intervals are accurate and the precision of the 3-level co-kriging is significantly better than the one of the 2-level co-kriging.

Conclusion

We have presented a method for building kriging models using a hierarchy of codes with different levels of accuracy. This method allows us to improve a surrogate model built on a complex code using information from a cheap one. It is particularly useful when the complex code is very expensive. We see in our literature review that the first multi-level metamodel originally suggested is a first-order auto-regressive model built with Gaussian processes. The AR(1) relation between two levels of code is natural and the building of the model is straightforward. Nevertheless, we have highlighted some key issues which makes it difficult to use this model in practical ways. First, important parameters of the model, which are the adjustment coefficients between two successive levels of codes, were numerically estimated. We propose here an analytical estimation of these parameters with a Bayesian method. This method allows us to have information about the uncertainties of the estimations and above all, to easily use the AR(1) model and its generalization to the case of non-spatial stationarity. Furthermore, a strength of the proposed method is that it even works for a code with more than 2 levels since its implementation is such that the estimations of the parameters of a s-level co-kriging is equivalent to the ones of s independent krigings. It is important to highlight that this method is based on a joint estimation between the adjustment coefficient and the mean of the Gaussian process modelling the difference between two successive levels of code.

Second, we have seen that the variance of the predictive distribution of the AR(1) model could be underestimated. A natural approach to improve this estimation is a Bayesian modelling. We propose here a Bayesian co-kriging for 2 levels of code and to avoid computationally expensive implementation, we suggest another model than the one presented. This new model is based on a hierarchical specification of the parameters of the model. This allows us to have a Bayesian model including only two nested integrations without Markov chain Monte Carlo procedure.

Finally, for a non-Bayesian s-level co-kriging, we have proved that building a s-level co-kriging is equivalent to build s independent krigings. This result is very important since it solves one of the most important key issues of the co-kriging which is the inversion of the covariance matrix. A 3-level co-kriging example has been provided to show the efficiency of the presented method.

Acknowledgements

The author thanks his supervisor Josselin Garnier for his valuable guidance. He also thanks Claire Cannamela for her advice and constructive suggestions.

A The case of ρ depending on x

A.1 Building a model with s levels of code Let us consider s levels of code, if we note β = (β T 1 , . . . , β T s ) T , β ρ = (β T ρ 1 , . . . , β T ρ s-1 ) T , σ 2 = (σ 2 1 , . . . , σ 2 s ) and θ = (θ 1 , . . . , θ s ), we have [Z s (x)|Z = z, β, β ρ , σ 2 , θ] ∼ N m Zs (x), s 2 Zs (x) where m Zs (x) and s 2 Zs (x) are defined in equations ( 29) and (30). Let us define the notation l i=k A i = A k ⊙ • • • ⊙ A l where ⊙ represents the matrix element-by-element product. Furthermore, let us denote by ρ t = ρ t (D t ) the vector containing the values of ρ t (x), x ∈ D t . The s diagonal blocks of V s (31) of size n t × n t are defined by:

V (t,t) = σ 2 t R t (D t )+σ 2 t-1 ρ t-1 (D t )ρ T t-1 (D t ) ⊙R t-1 (D t )+• • •+σ 2 1 t-1 i=1 ρ i (D t )ρ T i (D t ) ⊙R 1 (D t )
and the off-diagonal blocks of size n t × n t ′ are given by: A.3 Some important results about the covariance matrix V s By sorting the experimental design sets as in Subsection 7.2, it can be shown that ∀t = 2, . . . , s the inverse of the matrix V s has the form:

V (t,t ′ ) =   1 nt t ′ -1 i=t ρ i (D t ′ ) T   ⊙ V (t,
V -1 s =     V -1 s-1 + 0 0 0 (ρ s-1 (D s )ρ T s-1 (D s )) ⊙ R -1 with V -1 1 = R -1 1 σ 2 1
, V -1 s-1 an ( s-1 i=1 n i × s-1 i=1 n i ) matrix and R -1 s an (n s × n s ) matrix. It can also be shown that:

t s (x) T V -1 s = ρ s-1 (x)t T s-1 (x)V -1 s-1 -[0, ρ T s-1 (D s ) ⊙ R s ({x}, D s )R -1 s ],R s ({x}, D s )R -1 s A.4
Bayesian prediction for a code with 2 levels

The equations for the Bayesian prediction when ρ depends on x can be directly derived from the Section 5 by replacing ρ with β ρ and noting that the design matrix F is such that:

F = [F ρ (D 2 ) ⊙ (z 1 (D 2 )1

Figure 1 :

 1 Figure 1: Toy examples. The co-kriging metamodel is very close to the expensive output z 2 (.) and improves significantly the ordinary kriging metamodel using the small design D 2 .

Figure 2 :

 2 Figure2: Toy examples. The high frequency components of the expensive code are not predicted since they are not captured by the cheap code and the coarse grid used for the expensive code cannot detect them either. Nevertheless, the co-kriging improves the ordinary kriging metamodel since the cheap code allows us to predict the low frequencies of the expensive code accurately.

Figure 3 :

 3 Figure 3: A toy example without any prior information. The thick dotted line represents the prediction mean, the thin dotted lines represent the confidence interval at plus or minus twice the standard deviation in the non-Bayesian case and the dashed lines represent the same confidence interval in the Bayesian case.

3 , a = 7 and b = 1 / 10 .

 3110 We note that the complex function z 3 (x) corresponds to the Ishigami function which is very

Figure 4 :

 4 Figure 4: CPU time comparison between 2 level co-kriging models. The triangles represent the CPU time for the crude co-kriging model and the circles represent the CPU time for the co-kriging model with the complexity reduction. The gain of CPU time with the reduction complexity is approximately a factor equal to 1.93.

Figure 5 :

 5 Figure 5: Predictions of the 2-level and the 3-level co-krigings for the fluidized-bed process. The 3-level co-kriging improves significantly the predictions of the 2-level one.

Figure 6 :

 6 Figure6: Prediction errors of the 2-level co-kriging and confidence intervals at plus or minus twice the standard deviation. We see a significant difference between the accuracy of the predictions means and their confidence intervals for the point where T 3 is unknown (the 8 first validation points) and for the ones where it is known (the last 10 validation points).

Figure 7 :

 7 Figure7: Leave-One-Out Cross-Validation predictive errors and variances of the 2-level and 3-level co-kriging. We see that the confidence intervals are accurate and the precision of the 3-level co-kriging is significantly better than the one of the 2-level co-kriging.

  2 2sup |z 1 , z 2 ) are close to 0. This region essentially contains the support of the function. Furthermore, we create a non-uniform integration grid distributed with a geometric progression.Finally p(z 2 (x)|z 1 , z 2 ) is a predictive density function integrating the posterior distribution of parameters (β 2 , ρ, β 1 , σ 2 1 , σ 2 2

	the equation (24) and such that p(σ 2 1	2 1 inf , σ 2 1sup ] × [σ 2 2 inf , σ 2 2sup ] from

inf |z 1 ), p(σ 2 1sup |z 1 ) p(σ 2 2 inf |z 1 , z 2 ) and p(σ

Table 1 :

 1 A co-kriging example with one-variable functions. Bayesian estimation of parameters.

	Regression Coefficient Estimation Variance Coefficient Estimation
	ρ β 2 β 1	2 (20, -20) -3.49	σ 2 1 σ 2 2	32.75 7.02 × 10 -30

Table 2 :

 2

	Regression Coefficient	Estimation	Variance Coefficient Estimation
	ρ β 2 β 1	1.86 (18.39, -17.00) -3.49	σ 2 1 σ 2 2	32.75.03 0.30

A co-kriging example with one-dimensional functions. Bayesian estimation of parameters.

Table 5 :

 5 Example: fluidized-bed process. Estimation of the hyperparameters (correlation lengths) for the 3-level co-kriging.

	1.790 3.988 1.218 1.790 3.595 0.722
	θ2 1.810 1.842 2.008 1.036 0.001 0.345
	θ3 0.890 0.721 2.008 2.952 1.790 0.241

Table 6 :

 6 Example: fluidized-bed process. Bayesian estimation of the variance and regression parameters for the 3-level co-kriging.

	Posterior Covariance
	σ 2 t

Table 7 :

 7 Example: fluidized-bed process. Comparison between 2-level cokriging and 3-level co-kriging. Predictions are better in the 3-level case and the prediction variance seems well-evaluated since the RMSE and the average standard deviation are close.

	MaxAE

  t) (D t , D t ′ ) 1 ≤ t < t ′ ≤ sThe vector t s (x) in equations (29) and (30) is such that t s (x) = (t * 1 (x, D 1 ) T , . . . , t * s (x, D s ) T ) T , where:t * t (x, D t ) T = ρ T t-1 (D t ) ⊙ t * t-1 (x, D t ) T + R 1 (x, D 1 ). Furthermore, the matrix H s in equations 35 can be written as:(D j ) 1 T p 1 ⊙ F 1 (D j ) ⊙ F 2 (D j ) . . . F j (D j ) 0 . . . . .. Bayesian estimation of parameters for s levels of codeWe can extend the Bayesian estimation of the parameters to the case of ρ depending on x. Note that we do not assume the independence of β t and β ρ t-1 . We have:[(β ρ t-1 , β t )|z t , z t-1 , θ t , σ 2 p tq t-1 2 Q t = (z t -H t λt ) T R -1 t (D t )(z t -H t λt ) λt = E[(β ρ t-1 , β t )|z t , z t-1 , θ t , σ 2 t ]The REML estimator of σ 2 t is σ2 t = Qt 2αt and we can estimate θ t by minimizing the expression: log(|det(R t (D t ))|) + (n tp tq t-1 )log(σ 2 t )

			s-1
				ρ i (x) σ 2 t R t (x, D t )
			i=t
	where 1 < t ≤ s,  1 H s = s-1 i=s ρ i (x) = 1 and t * 1 (x, D 1 ) T = s-1 i=1 ρ i (x) σ 2 . . .    j-1 i=1 ρ i j-1 i=2 ρ i (D j )1 T p 2	. . .	   
				
	A.2 Q t 2	)
	where	α t =	n t

t ] ∼ N H T t R -1 t (D t )H t -1 H T t R -1 t (D t )z t , σ 2 t H T t R -1 t (D t )H t -1

where

H t = [F ρ t-1 (D t ) ⊙ (z t-1 (D t )1 T q t-1 ) F t (D t )]

. Furthermore, we have:

[σ 2 t |z t , z t-1 , θ t ] ∼ IG(α t ,

  T pρ ) F 2 ]Finally, for the Bayesian prediction, we just have to adapt the integral (27) :p(z 2 (x)|z 1 , z 2 , σ 2 1 , σ 2 2 ) = R pρ+p 2 p(z 2 (x)|z 1 , z 2 , β 2 , β ρ , σ 2 1 , σ 2 2 )p(β ρ , β 2 |z 1 , z 2 , σ22 ) dβ ρ dβ 2

whereas in[START_REF] Qian | Bayesian Hierarchical Modeling for Integrating Low-accuracy and High-accuracy Experiments[END_REF] they use separated prior distributions with ρ not conditioned


Pierre Dragicevic

Yuanchun Shi

Visualizing and Manipulating Automatic Document Orientation Methods Using Vector Fields

Keywords: Tabletop, orientation, vector fields. ACM Classification H5.m. Information interfaces and presentation

We introduce and illustrate a design framework whereby tabletop documents are oriented according to vector fields that can be visualized and altered by end users. We explore and illustrate the design space using interactive 2D mockups and show how this approach can potentially combine the advantages of the fully manual and fully automatic document orientation methods previously proposed in the literature.

INTRODUCTION AND RELATED WORK

Orienting documents is an old problem in tabletop systems and many strategies have been explored [START_REF] Hancock | Rotation and Translation Mechanisms for Tabletop Interaction, tabletop[END_REF]. They fall into two categories: explicit orientation methods and automatic orientation methods.

Explicit orientation methods provide widgets or gestures that let users rotate documents. Although simultaneous control of rotation and translation is achievable [START_REF] Kruger | Fluid integration of rotation and translation[END_REF][START_REF] Wu | Multi-finger and whole hand gestural interaction techniques for multi-user tabletop displays[END_REF], having to constantly orient and reorient documents can be tedious and distract users from their primary task. Automatic orientation methods free users from the burden of having to rotate documents [START_REF] Shen | DiamondSpin: an extensible toolkit for around-the-table interaction[END_REF][START_REF] Vernier | Visualization techniques for circular tabletop interfaces[END_REF]: each document is adequately oriented depending on its location, allowing users to solely focus on positioning (Figure 1).

One issue with automatic orientation is that it is rarely optimal. For example, circular orientation -i.e., having documents face away from the center of the table -works well for round tables (Figure 1a) but not for rectangular tables (Figure 2). Although orientation methods can be carefully designed for specific situations, they will always lack the flexibility and the expressiveness provided by explicit rotation, especially when collaborating with other users [START_REF] Kruger | How people use orientation on tables: comprehension, coordination and communication[END_REF].

To sum up, none of the two strategies proposed so far is satisfactory. Explicit rotation can be tedious and distract users from their primary task. Automatic orientation is convenient, but a method that works well in one case will fail as soon as the contexte.g., the shape of the table, the number and location of users, their tasks or collaborative statuschanges. This leads us to propose a middle-ground approach, where users are given control over automatic orientation methods. This approach first requires automatic orientation methods to be reified into concrete entities. We use vector fields, i.e., functions that map locations to 2-D vectors. Vector fields are a popular mathematical construct, with a well-defined algebra. Besides, they are widely used to describe physical phenomena such as fluid motion or gravitational fields, and hence have well-known visual representations.

Permission to make digital or hard copies of all or part of this work for personal or classroom use is granted without fee provided that copies are not made or distributed for profit or commercial advantage and that copies bear this notice and the full citation on the first page. To copy otherwise, or republish, to post on servers or to redistribute to lists, requires prior specific permission and/or a fee. (1)

Since vector fields encode document orientation, we call them orientation fields. The design space of the orientation field framework can be split up into three separate problems: a) generation of orientation fields; b) visualization of orientation fields; and c) user manipulation of orientation fields. We examine them in order, before providing a few general illustrations.

GENERATING ORIENTATION FIELDS

Default orientation fields need to be generated before they can be edited. Existing automatic orientation methods (Figure 1) can be described with vector fields, using either continuous functions or a bilinear interpolation of discrete vector fields. On a tabletop system, several such fields can be provided as default configurations. Users would then select them and possibly modify them according to their needs. Modified fields could in turn be saved and retrieved across working sessions.

It is also possible to automatically produce new orientation fields from existing ones. Two potentially useful operations are filtering and combination.

Filtering: Moving documents can sometimes yield discontinuous changes in orientation (Figure 3a). Since abrupt visual changes are known to confuse users [START_REF] Chang | Animation: from cartoons to the user interface[END_REF], lowpass filtering techniques can be used to eliminate these discontinuous changes. Smooth orientation fields additionally offer more flexibility for orienting documents: the low-pass-filtered closest border method, for example, allows diagonal orientation (Figure 3b).

Combination: New orientation fields can be produced by combining existing orientation fields. Such operations are useful when a consensus needs to be reached, e.g., when several users having their own preferred field happen to sit together around a table, or when several tables are physically assembled [START_REF] Shi | Finger gesture interaction on large tabletop for sharing digital documents among multiple users[END_REF]. Orientation field combination is also useful for smoothly animating discrete changes in orientation fields, e.g., when browsing orientation fields from a repository.

In practice, orientation fields can be combined by computing a linear combination of vectors in either Cartesian or polar coordinates.

When Cartesian coordinates are used, vector length encodes orientation weight or "importance". This allows users to, e.g., define a personal area [START_REF] Scott | Territoriality in collaborative tabletop workspaces[END_REF] whose orientation field will not be affected when other users join the table (very long vectors) while leaving the rest open to other users (null-length vectors).

When polar coordinates are used, all vector orientations have equal importance. Vector length can be set to an arbitrary constant or used to encode additional information such as document size [START_REF] Vernier | Visualization techniques for circular tabletop interfaces[END_REF].

VISUALIZING ORIENTATION FIELDS

Since users will be allowed to modify orientation fields, they need visual feedback on the effects of their actions. They also need the behavior of the orientation field to be always visible to be able to use it efficiently. Although this behavior can be suggested by specific visual designs (Figure 1a), we need to be able to visualize arbitrary fields.

Displaying vector fields is a well-researched topic in scientific visualization, with methods ranging from simple vector grids to more sophisticated textures and streamlines [START_REF] Laidlaw | Comparing 2D Vector Field Visualization Methods: A User Study[END_REF]. We experimented with two simple ones: vector grid and direction-to-hue.

Vector grid. This classic technique shows vectors sampled on a regular grid. Rotating vectors by 90° yields a different visual effect (Figures 3 and4a).

Direction-to-hue. Orientations naturally map to color hues, since hue is cyclic and typically expressed in angular units (Figure 4b). Vector lengths can be either ignored or mapped to other color components such as brightness or saturation.

We found direction-to-hue to be a promising approach. Since hue is preattentive, areas with homogeneous orientations pop out and appear as distinct regions (Figure 4b). As we will further see, this allows to rapidly identify meaningful areas.

MANIPULATING ORIENTATION FIELDS

There are many possible ways to edit orientation fields, and the best methods depend on the hardware available. We explored a technique for the uTable system, which supports direct touch and laser pointing [START_REF] Shi | Finger gesture interaction on large tabletop for sharing digital documents among multiple users[END_REF]. Our technique lets users specify centripetal orientation fields -i.e., all vectors are oriented towards a pointwithin regions of arbitrary shape. The technique requires performing two successive gestures, one for specifying the region to be modified, the other one for indicating the center of attraction (Figure 5a).

The system connects the two ends of the first gesture to the nearest edges of the table and closes the curve. The center of attraction C is computed from the two ends P0 and P1 of the second gesture according to the following formula:

C = k • P0P1. Choosing k > 1 allows C to be off-table,
which can be desirable in some cases. For example, personal documents might be more conveniently oriented if they are facing a point behind the observer (Figure 5c) rather than a point on the table (Figure 5b).

ILLUSTRATIONS

We illustrate the use of orientation fields for visualizing and manipulating document orientation with three simple scenarios. These involve the smoothing filter, the directionto-hue visualization and the centripetal reorientation technique described in the previous sections. We assume a table equipped with a laser pen [START_REF] Shi | Finger gesture interaction on large tabletop for sharing digital documents among multiple users[END_REF] with which users can drag documents anywhere on the table. The default orientation field used by the table is a low-pass-filtered closest border method.

Creating a Personal Space

John enters the room and sits at the table next to his colleague Mike (Figure 6a). He starts opening documents but some of them are inconveniently oriented (Figure 6b), so he creates a personal working area oriented towards himself (Figure 6cd). All his documents are immediately reoriented, and all documents he will subsequently drop within that space will be oriented the same way.

If other people join the table, they will see the space John occupies -the red area -even if there is no document on it. This will help them choosing where to sit and where to create their own personal space.

Creating a Shared Space

John and George are inspecting resumes from prospective graduate students and decide to select the top six (Figure 7a). John drags the resume of his favorite candidate to an empty area so that both George and him can see it, but he notices it is upside-down (Figure 7b). With two rapid gestures he re-orients the whole area, turning it into a shared working space (Figure 7c). Orienting the remaining documents is quite straightforward since it only involves dropping them to the green area (Figure 7d).

Temporary Reorientation

Three people are working on an interactive table they reconfigured to suit their needs (Figure 8). John wants to bring Paul's attention to two of his documents. He reorients them with two rapid gestures. Upon reading them, Paul reorients the documents towards John.

PROTOTYPE IMPLEMENTATION

All the figures except Figure 1 are screenshots of juTable, a Java library for creating interactive tabletop mockups running on desktop computers, with user avatars that can be freely arranged around tables of arbitrary size. Such simulations make it much easier to prototype and communicate new interaction ideas. We are currently porting JuTable to the uTable system [START_REF] Shi | Finger gesture interaction on large tabletop for sharing digital documents among multiple users[END_REF] so that mockups can be directly deployed in a real setting. The orientation field demo and source code can be downloaded at http://www.aviz.fr/jutable.

CONCLUSION AND FUTURE WORK

We proposed a conceptual framework for managing document orientation, where orientation is controlled by vector fields that can be visualized and edited. This strategy combines the advantages of previous methods: like automatic methods, it maps locations to orientations so that users do not have to repeatedly reorient documents. But this mapping can also be opportunistically modified by users, hence approaching the flexibility of manual rotation.

One noteworthy property of the technique is its visibility: color-to-hue visualization makes areas with homogeneous orientation clearly visible, which turns rotations tasks into simple spatial target acquisition tasks and is likely to promote territoriality and group awareness [START_REF] Scott | Territoriality in collaborative tabletop workspaces[END_REF].

The ideas we proposed have been explored and illustrated using mockups, and need to be tested and validated in physically realistic setups. A better understanding on how people use document orientation is also critical to validate the approach [START_REF] Kruger | How people use orientation on tables: comprehension, coordination and communication[END_REF].

The design space can be further explored in several directions. This includes testing more sophisticated visualizations [START_REF] Laidlaw | Comparing 2D Vector Field Visualization Methods: A User Study[END_REF] and exploring physical metaphors such as fluids and forces. For example, users could place magnetic or gravitational objects on the table to modify vector fields, or manipulate conservative vector fields as 3D surfaces. Another open problem is the integration of explicit rotation with automatic rotation techniques.

  TableTop 2009 November 23-25, 2009, Banff, Alberta Copyright 2008 ACM XXXXXXXXXXXX...\$5.00.

Figure 1 .

 1 Figure 1. Automatic document orientation in DiamondSpin, on a round table and on a square table. Image courtesy of [8].

Figure 2 .

 2 Figure 2. With a circular (a) or elliptic (b) orientation scheme, the user sitting on the corner does not have his documents correctly oriented ; Even with the closest border method (c), two of his documents are still inconveniently oriented.

Figure 3 .

 3 Figure 3. The closest border method, without (a) and with (b) low-pass filtering.

Figure 4 .

 4 Figure 4. Visualizing elliptic, closest-border, and smoothed closest border orientation fields using vector grids (a) and direction-to-hue (b). (best viewed in color)

Figure 5 .

 5 Figure 5. Reorienting a region, with k = 1 (the center of attraction is on the table's edge) and k = 4 (the center of attraction is behind the user).

Figure 6 .

 6 Figure 6. Creating a personal working space.

Figure 7 .

 7 Figure 7. Creating a shared working space.