process regression Kriging noisy data Parameter

Keywords:

estimation Model validation A few GPR topics beyond basics Kernel design Two other points of view Links with other methods Kriging issues

Context

Kriging is most often used in the context of expensive (numerical) experiments (simulators, e.g. PDE solvers):

The experiment can be seen as a function of the input parameters

y = f (x)
where f is a costly to evaluate function.

In the following, we will assume that

x ∈ X : There are d input variables. Usually (but not necessarily) X is R d . y ∈ R: The output is a scalar. But extensions to GP regression with multiple outputs exist.

The fact that f is costly to evaluate changes a lot of things... 1. Representing the function is not possible... The fact that f is costly to evaluate changes a lot of things...

2. Uncertainty propagation is not possible... The fact that f is costly to evaluate changes a lot of things...

3.

Optimisation is also tricky...

x > > f + + + + + +
4. Computing integrals is not possible... 5. Sensitivity analysis is not possible...

Statistical modelling

We know an initial Design of Experiments (DoE) of n points (x i , y i ), y i = f (x i ). What can be said about possible y at any x using probabilities? ⇒ kriging for regression (conditional GP) [START_REF] Krige | A statistical approach to some basic mine valuation problems on the witwatersrand[END_REF][START_REF] Matheron | Principles of geostatistics[END_REF] = a family of surrogates (metamodels) with embedded uncertainty.

x i y i

General bibliography: [START_REF] Rasmussen | Gaussian Processes for Machine Learning[END_REF], Durrande and Le Riche, 2017[START_REF] Le Riche | Introduction to Kriging[END_REF] Random Process (1/2)

Random variable Y random event ω ∈ Ω (e.g., throw a dice) =⇒ get an instance y . Ex:

if dice ≤ 3 , y = 1 if 4 ≥ dice ≤ 5 , y = 2 if dice = 6 , y = 3 Random process Y (x)
A set of RV's indexed by x random event ω ∈ Ω (e.g., weather) =⇒ get a function y (x). Ex:

Le Riche & Durrande GP regression Sept. 2019 Gaussian Process (1/2) Assume Y () is a GP, Y (x) ∼ N (µ(x), k(x, x)) ⇔ ∀X =   x 1 . . . x n   ∈ X n×d , Y (X ) =   Y (x 1 ) . . . Y (x n )   ∼ N (µ(X ), K )
where

K ij = Cov(Y (x i ), Y (x j )) = k(x i , x j ) depends only on the x's y 1 y 2 y 3 Y 1 Y 2 Y 3 y (x) : µ(x i ) K 12 K 23 K 13 Le Riche & Durrande GP regression Sept. 2019
Gaussian Process (2/2)

The distribution of a GP is fully characterised by: its mean function µ(.) defined over X its covariance function (or kernel) k(., .) defined over X × X :

k(x, x ) = Cov(Y (x), Y (x )) ⇒ Example path simulation: Say k(x, x ) = σ 2 exp -(x -x ) 2 /θ 2
and in pseudo-R, build a fine grid X, choose mean function mu(), build the covariance matrix, K[i,j]=k(X[i],X[j]), eigenanalysis, Keig = eigen(K), and sample,

y = mu[X] + Keig$vectors %*% diag(sqrt(Keig$values)) %*% matrix(rnorm(n))
⇒ See also Shiny App: https://github.com/NicolasDurrande/shinyApps

Valid kernels

A kernel satisfies the following properties:

It is symmetric: k(x, x ) = k(x , x) It is positive semi-definite (psd): ∀n ∈ N, ∀x i ∈ D, ∀α ∈ R n , n i=1 n j=1 α i α j k(x i , x j ) ≥ 0
Furthermore any symmetric psd function can be seen as the covariance of a Gaussian process. This equivalence is known as the Loeve theorem.

Popular kernels in 1D

There are a lot of functions that have already been proven psd:

constant k(x, x ) = σ 2 white noise k(x, x ) = σ 2 δ x,x (Kronecker delta function) Brownian k(x, x ) = σ 2 min(x, x ) power-exponential k(x, x ) = σ 2 exp (-|x -x | p /θ) , 0 < p ≤ 2 Matérn 3/2 k(x, x ) = σ 2 (1 + |x -x |) exp (-|x -x |/θ) Matérn 5/2 k(x, x ) = σ 2 1 + |x -x |/θ + 1/3|x -x | 2 /θ 2 × exp (-|x -x |/θ) squared exponential k(x, x ) = σ 2 exp -(x -x ) 2 /θ 2 linear k(x, x ) = σ 2 xx . . .
The parameter σ 2 is called the variance and θ the length-scale.

General factorized form: k(x, x ) = σ 2 r (x, x ), r (, ) the correlation function.

Regularity and covariance function

The regularity and frequency content of the y (x) are controlled by the kernel (and its length-scale) For stationary processes (depend on τ = xx only), the trajectories are p times differentiable (in the mean square sense) if k(τ ) is 2p times differentiable at τ = 0 ⇒ the property of k(τ ) at τ = 0 define the regularity of the process. Examples:

trajectories with squared exponential kernels are infinitely differentiable = very (unrealistically?) smooth. trajectories with Matérn 5/2 and 3/2 kernels are twice and once differentiable. trajectories with power-exponential are not differentiable excepted when p = 2.

Popular multi-dimensional kernels (1/2)

constant k(x, x ) = σ 2 white noise k(x, x ) = σ 2 δ x,x exponential k(x, x ) = σ 2 exp (-||x -x || θ ) Matérn 3/2 k(x, x ) = σ 2 1 + √ 3||x -x || θ exp - √ 3||x -x || θ Matérn 5/2 k(x, x ) = σ 2 1 + √ 5||x -x || θ + 5 3 ||x -x || 2 θ × exp - √ 5||x -x || θ sq. exp. k(x, x ) = σ 2 exp - 1 2 ||x -x || 2 θ where ||x -x || θ = d i=1 (x i -x i ) 2 θ 2 i 1/2 .
A common general recipee: a product of univariate kernels,

k(x, x ) = σ 2 d i=1 r i (x i , x i )
which has d + 1 parameters.

(more on kernel design later)

Le Riche & Durrande GP regression Sept. 2019

Gaussian process regression

Assume we have observed a function f () over a set of points X = (x 1 , . . . , x n ):

0 5 10 15 x -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 2.0 2.5 f
The vector of observations is

F = f (X ) (ie F i = f (x i ) ). Le Riche & Durrande GP regression Sept. 2019
Since f () in unknown, we make the general assumption that it is the sample path of a Gaussian process Y ∼ N (µ(), k(, )):

0 5 10 15 x -4 -3 -2 -1 0 1 2 3 4 Y (here µ(x) = 0) Le Riche & Durrande GP regression Sept. 2019
If we remove all the samples that do not interpolate the observations we obtain: Kriging equations (1/2)

The conditional distribution can be obtained analytically:

By definition, (Y (x), Y (X )) is multivariate normal. Formulas on the conditioning of Gaussian vectors give the distribution of

Y (x)|Y (X ) = F . It is N (m(.), c(., .
)) with :

m(x) = E[Y (x)|Y (X )=F ] = µ(x) + k(x, X )k(X , X ) -1 (F -µ(X )) c(x, x ) = Cov[Y (x), Y (x )|Y (X )=F ] = k(x, x ) -k(x, X )k(X , X ) -1 k(X , x ) Le Riche & Durrande GP regression Sept. 2019 28 / 74
Kriging equations (2/2)

The distribution of Y (x)|Y (X ) = F is N (m(.), c(., .)) with:

m(x) = E[Y (x)|Y (X )=F ] = µ(x) + k(x, X )k(X , X ) -1 (F -µ(X )) c(x, x ) = Cov[Y (x), Y (x )|Y (X )=F ] = k(x, x ) -k(x, X )k(X , X ) -1 k(X , x ) k(X , X ) = [k(x i , x j )] : covariance matrix, Gram matrix in SVM. k(x, X ) = [k(x, x 1 ), . . . , k(x, x n )] : covariance vector, only dependance on x beside µ(x).
It is a Gaussian distribution: gives confidence intervals, can be sampled, this is actually how the previous slides were generated. Bayesian:

Y (x)|Y (X ) = F is the posterior distribution of Y (x) once Y (X ) = F is observed.
A few remarkable properties of GPR models They (can) interpolate the data-points The prediction variance does not depend on the observations The mean predictor does not depend on the variance parameter They (usually) come back to the a priori trend µ(x) when we are far away from the observations.

(proofs left as exercise)

Le Riche & Durrande GP regression Sept. 2019

Changing the kernel has a huge impact on the model:

Gaussian kernel: Exponential kernel:
This is because changing the kernel means changing the prior on f

Gaussian kernel:

Exponential kernel:

There is no kernel that is intrinsically better... it depends!

Gaussian kernel: Exponential kernel:

The kernel has to be chosen according to the prior belief on the behaviour of the function to study: is it continuous, differentiable, how many times? is it stationary ? is it monotonous, bounded? Cf. [START_REF] López-Lopera | Finite-dimensional gaussian approximation with linear inequality constraints[END_REF] ... (more on this in the kernel design section later) Default: constant trend µ (empirical mean or µ from max likelihood [START_REF] Roustant | [END_REF]) and Matérn 5/2 kernel.

Examples of models with observation noise for n(x, x ) = τ 2 δ x,x :

0.0 0.2 0.4 0.6 0.8 1.0

0 1 2 3 x Z(x)|Z(X) + N (X) = F 0.0 0.2 0.4 0.6 0.8 1.0 0 1 2 3 x Z(x)|Z(X) + N (X) = F 0.0 0.2 0.4 0.6 0.8 1.0 -1 0 1 2 3 4 x Z(x)|Z(X) + N (X) = F
The values of τ 2 are respectively 0.001, 0.01 and 0.1.

Kriging with noise kernel (nugget) does not interpolate the data.

A small τ 2 (e.g., 10 -10 ) often used to make the covariance matrix invertible (more on regularization of GPs in [Le Riche et al., 2017]).

Parameter estimation

We have seen previously that the choice of the kernel and its parameters (σ 2 , the θ's, the trend and other parameters) have a great influence on the model.

In order to choose a prior that is suited to the data at hand, we can: minimise the model error maximize the model likelihood We now detail the second approach.

Definition: The likelihood of a distribution with a density p U given observations u 1 , . . . , u p is:

L = p i=1 p U (u i )
The likelihood measures the adequacy between observations and a distribution.

In the GPR context, we often have only one observation of the vector F . The likelihood is then:

L = p Y (X ) (F ) = 1 (2π) n/2 det(k(X , X )) 1/2 × exp - 1 2 (F -µ(X )) k(X , X ) -1 (F -µ(X )) .
It is thus possible to maximise L -or log(L) -with respect to the kernel and model parameters in order to find a well suited prior.

The likelihood in a multi-modal function in θ's and must be optimized with global optimization algorithms.

(more details on likelihood such as concentration in, e.g.

[Le Riche, 2014])

We have seen that given some observations F = f (X ), it is very easy to build lots of models, either by changing the kernel parameters or the kernel itself.

The question is now how to measure the quality of a model to build the best one at the end. Principle: introduce new data and to compare them to the model prediction.

0.0 0.2 0.4 0.6 0.8 1.0

-1 0 1 2 x Z(x)|Z(X) = F Le Riche & Durrande GP regression Sept. 2019
Let X t be the test set and F t = f (X t ) be the associated observations.

The accuracy of the mean can be measured by computing:

Mean Square Error MSE = mean((F t -m(X t )) 2 ) A "normalized" criterion Q 2 = 1 - (F t -m(X t )) 2 (F t -mean(F t )) 2
On the above example we get MSE = 0.038 and Q 2 = 0.95. The predicted distribution can be tested by normalizing the residuals.

According to the model,

F t ∼ N (m(X t ), c(X t , X t )). c(X t , X t ) -1/2 (F t -m(X t
)) should thus be independents N (0, 1):

standardised residuals Density -3 -2 -1 0 1 2 3 0.0 0.1 0.2 0.3 0.4 0.5 0.6 -2 -1 0 1 2 -3 -2 -1 0 1 2 3 Normal Q-Q Plot Theoretical Quantiles Sample Quantiles Le Riche & Durrande GP regression Sept. 2019
When no test set is available, another option is to consider cross validation methods such as leave-one-out.

The steps are:

1. build a model based on all observations except one 2. compute the model error at this point This procedure can be repeated for all the design points in order to get a vector of error. Model to be tested:

0.0 0.2 0.4 0.6 0.8 1.0 -2 -1 0 1 x Z(x)|Z(X) = F Le Riche & Durrande GP regression Sept. 2019 43 / 74
Step 1:

0.0 0.2 0.4 0.6 0.8 1.0 -2 -1 0 1 x Z(x)|Z(X) = F Le Riche & Durrande GP regression Sept. 2019 44 / 74
Step 2:

0.0 0.2 0.4 0.6 0.8 1.0 -2 -1 0 1 x Z(x)|Z(X) = F Le Riche & Durrande GP regression Sept. 2019 45 / 74
Step 3:

0.0 0.2 0.4 0.6 0.8 1.0 -2 -1 0 1 x Z(x)|Z(X) = F Le Riche & Durrande GP regression Sept. 2019
We finally obtain: MSE = 0.24 and Q 2 = 0.34.

We can also look at the residual distribution. For leave-one-out, there is no joint distribution for the residuals so they have to be standardized independently.

standardised residuals Density -2 -1 0 1 2 3 0.0 0.1 0.2 0.3 0.4 -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 -2 -1 0 1 2 Normal Q-Q Plot Theoretical Quantiles Sample Quantiles Le Riche & Durrande GP regression Sept. 2019
Sample code in R (with 6D Hartman function) library(DiceKriging) library(DiceDesign) X <-lhsDesign(n=80,... dimension=6)$design X <-data.frame(X) y <-apply(X, 1, hartman6) mlog <-km(design = X, response = -log(-y)) plot(mlog) Trajectories of Y () and Y 1 () can be obtained from each other: Symmetric kernel: to have Y

Y (x) | Y (X ) = F sampled through Y 1 () with h(x)Y 1 (x) | Y 1 (X ) = F /h(X ) (component-
x 1 x 2 = Y x 2 x 1 , use k x 1 x 2 , x 1 x 2 = k x 1 x 2 , x 1 x 2 + k x 2 x 1 , x 1 x 2

Example

Symmetrical kriging, mean and std. deviation.

Note how the variance is null symmetrically to observations. from [Ginsbourger, 2009] The statistical point of view

Kriging is often introduced as best linear interpolator: linear:

Y (x) = n i=1 λ i (x)Y (x i ) = λ λ λ(x) Y (X ) unbiased: E Y (x) = λ λ λ(x) EY (X ) = EY (x) = µ(x) best: λ λ λ(x) = arg min λ λ λ∈R n E Y (x) -Y (x) 2
This constrained optimization problem is solved in λ λ λ(x) and the kriging equations are recovered from

m(x) = E( Y (x) | Y (X ) = F ) and c(x, x ) = E ( Y (x) -Y (x))( Y (x ) -Y (x ))
but the link with the GP interpretation is typically not discussed,

E(Y (x) | Y (X ) = F ) ? = E( Y (x) | Y (X ) = F ) and E((Y (x) -µ(x)) 2 | Y (X ) = F ) ? = E( Y (x) -Y (x)) 2 .
Le Riche & Durrande GP regression Sept. 2019

⇒ Another way to make kernels (Mercer): choose the φ i (.)'s,

k(x, x ) = N i=1 λ i φ i (x)φ i (x )
Degenerated k(., .) if N < +∞ (the covariance matrix becomes non-invertible beyond N observations)

Proof: k(x, .) ∈ L 2 (X ), k(x, .) = i k(x, .), φ i (.) L 2 φ i (.) = i X k(x, t)φ i (t)dtφ i (.) = i λ i φ i (x)φ i (.)
Alternative definition of the RKHS:

H = {f (.) ∈ L 2 (X ) : f (.) = ∞ i=1 c i φ i (.) and ∞ i=1 c 2 i λ i < +∞}
i.e., impose a sufficiently fast decrease in eigencomponents, a kind of regularization.

Intuition behind the alternative definition of the RKHS, H: 1. construct an orthonormal basis of H

The φ i (.)'s are bi-orthogonal w.r.t. ., . L 2 and ., . H but need to be normalized in H: we use the integral equation of slide 61, which gives an intuition that the φ i ()'s belong to H (think of the integral as a sum). Then,

φ i , φ j H = 1 λ i λ j X k(., t)φ i (t)dt, X k(., t )φ j (t )dt = 1 λ i λ j X X φ i (t)φ j (t )k(t, t )dtdt = 1 λ i λ j X φ i (t) X k(t, t )φ j (t )dt dt = λ j λ i λ j φ i , φ j L 2 = λ j λ i λ j δ ij ⇒ φ i () = λ i φ i () is an orthonormal basis of H 2. f is in the RKHS if its coefficients are a converging series, f () = i c i φ i () = i c i √ λ i φ i () f () H = i c 2 i λ i < +∞ Le Riche & Durrande GP regression Sept. 2019

GPR and complexity control

If complexity is measured as the norm of the function, the Representer Theorem [START_REF] Schölkopf | A generalized representer theorem[END_REF] says that m() is the least complex interpolator:

m(.) = arg min h∈H h 2 H such that h(x i ) = f (x i ) , i = 1, . . . , n Proof: h = h k + h where h k = n i=1 c i k(x i , .
) and h ⊥ h k . Then, f (x i ) = k(x i , .), h k + h = k(x i , .), n j=1 c j k(x j , .) + 0

h ⊥h k = n j=1 c j k(x j , x i ) = k(x i , X )c, c the vector of n c i 's. The problem becomes, min c∈R n ,h ⊥span{k(x i ,.)} c k(X , X )c such that k(X , X )c = F whose solution is h = 0, c = k(X , X ) -1 F , i.e.,

m(x) = h k (x) = k(x, X )k(X , X ) -1 F
A regularization can also be seen in the likelihood in Slide 38 with det(k(X , X )) which must be as small as possible (⇒ large θ's). Other names, (almost) same equations m(x) = k(x, X )k(X , X ) -1 F is ubiquitous Bayesian linear regression: the posterior distribution is identical to the GPR equations under conditions on the kernel, cf.

[Le Riche, 2014] slide 35 and [START_REF] Rasmussen | Gaussian Processes for Machine Learning[END_REF], slide 20 of [Rosić, 2019]. Kalman filter, see slide 21 of [Rosić, 2019]. LS-SVR: same functional form of predictor (sum of kernels centered), but explicit regularization control (C , whereas GPR is implicit in likelihood), no uncertainty. RBF (Radial Basis Functions) [START_REF] Broomhead | Radial basis functions, multi-variable functional interpolation and adaptive networks[END_REF]: same prediction, no uncertainty (hence no likelihood).

  It can summarized by a mean function and 95% confidence intervals.

  x ) = h(x)h(x )k 1 (x, x ) is a valid kernel. Can be seen as k 1 (, )× composition of function and linear kernel. Better, see it as the covariance of Y (x) = h(x)Y 1 (x).

  wise division) Boundary conditions: say you want to impose that all trajectories go through Y (x) = 0 for all (infinite number) x's such that a(x) = 0. Use h(x) = d(a(x))

  infinite. Then, k(x, x ) = φ(x), φ(x )where ., . is the usual scalar product. A canonical distance in the mapped space (defined by the kernel) is,dist k (x, x ) = φ(x)φ(x ) 2 = k(x, x) + k(x , x ) -2k(x, x )

Sept. 2019

Repeat the random event (say 3 times): 3 y (x)'s. They are different, yet bear strong similarities.

Kriging of noisy data

An important special case, noisy data F = f (X ) + ε. Model F with Y (x) + N(x) where N(x) ∼ N (0, n(., .)) independent of Y (x). Then,

The expressions of GPR with noise become (just apply Gaussian vector conditioning with the above) Kernel design: making new from old

Many operations can be applied to psd functions while retaining this property Kernels can be: Summed

Multiplied

On the same space k(x, x

Composed with a function

to create new (non stationary) kernels, increase their dimension. All these transformations can be combined. Examples ...

Sum of kernels over the same space

Example (The Mauna Loa observatory dataset)

This famous dataset compiles the monthly CO 2 concentration in Hawaii since 1958.

Sum of kernels over the same space

We first consider a squared-exponential kernel with a small and a large length-scale: 1950 1960 1970 1980 1990 2000 Sum of kernels over the same space

What happens if we sum both kernels? Sum of kernels over the same space

We can try the following kernel:

The first term is a product of linear kernels. The periodic kernel is Composition with a function

Let k 1 be a kernel over X 1 × X 1 and h be an arbitrary function

This can be seen as a (non-linear) rescaling of the input space.

A way to make non-stationary kernels. The functional point of view (thanks Xavier Bay)

The kernel k(., .) defines a space of functions, a RKHS, H k := span{k(x, .), x ∈ X } with an inner product a ., . H such that there is a linear evaluation functional f (.), k(x, .) H = f (x).

a The inner product between 2 functions is:

), f (.), g (.) H = i,j α i β j k(x i , x j ), which implies the evaluation functional.

Associated to the psd k(., .) are eigenvalues and eigenfunctions

The φ i (.)'s form an orthonormal basis of H w.r.t. the usual scalar product. All this is a generalization of the eigendecomposition of symmetric positive definite matrices to infinite dimensions.

Trajectories can be generated with (Karhunen-Loève),

⇒ in general the trajectories are not in the RKHS:

But the GP mean is in the RKHS. [START_REF] Garland | Aerospace System Analysis and Optimization in uncertainty, chapter Cokriging for multifidelity analysis and optimization[END_REF], [START_REF] Fricker | Multivariate gaussian process emulators with nonseparable covariance structures[END_REF], with gradient as 2nd output [START_REF] Laurent | An overview of gradient-enhanced metamodels with applications[END_REF]. Discrete x variables: Cf. [START_REF] Roustant | Group kernels for gaussian process metamodels with categorical inputs[END_REF], mixed variables and optimization [START_REF] Pelamatti | Efficient global optimization of constrained mixed variable problems[END_REF].

Kriging issues

Too large n: k(X , X ) is n × n and takes O(n 3 ) operations for its inversion ⇒ not directly applicable beyond n = 1000. Solutions: inducing points [START_REF] Hensman | Gaussian processes for big data[END_REF], nested kriging [START_REF] Rullière | Nested kriging predictions for datasets with a large number of observations[END_REF]


Boost-HiC : Computational enhancement of long-range contacts in chromosomal contact maps

Genome-wide chromosomal contact maps are widely used to uncover the 3D organisation of genomes. They rely on the collection of millions of contacting pairs of genomic loci. Contact frequencies at short range are usually well measured in experiments, while there is a lot of missing information about long-range contacts.

We propose to use the sparse information contained in raw contact maps to determine high-confidence contact frequency between all pairs of loci. Our algorithmic procedure, Boost-HiC, enables the detection of Hi-C patterns such as chromosomal compartments at a resolution that would be otherwise only attainable by sequencing a hundred times deeper the experimental Hi-C library. Boost-HiC can also be used to compare contact maps at an improved resolution.

Boost-HiC is available at https://github.com/LeopoldC/Boost-HiC

Introduction

Chromosomal conformation capture (3C) has originally been developed to identify sets of DNA segments in close spatial proximity within a cell nucleus, and thus get an experimental in vivo access to the genome 3D organisation. It relies on the chemical fixation of chromosomal contacts, digestion with a restriction enzyme and a subsequent re-ligation of the cross-linked fragments [START_REF] Dekker | Capturing chromosome conformation[END_REF]. Next-generation sequencing techniques brought this protocol to the whole-genome scale (Hi-C) in cell populations [START_REF] Nynke | Hi-C: a method to study the three-dimensional architecture of genomes[END_REF] and an everyday increasing number of data-sets are now available. These datasets are usually difficult to compare since they can largely vary in terms of quality and sequencing depth. The highest resolution that can be achieved is in theory the size of the restriction fragments but very few datasets reach this resolution in practice. As a stunning example of pushing the experiments to its limits, Rao and colleagues provided the first 1kb-resolution contact map of the human genome, obtained by sequencing 4 billions of read pairs [START_REF] Suhas | A 3D map of the human genome at kilobase resolution reveals principles of chromatin looping[END_REF]. This dataset as been used in several studies since its release and many other data-sets have been produced in other human cell types and other organisms, but generally at a lower resolution [START_REF] Anthony D Schmitt | A compendium of chromatin contact maps reveals spatially active regions in the human genome[END_REF]. In mouse, another very high resolution data-set has recently been produced, covering the in vitro 1/10 differentiation of mouse neurons [START_REF] Bonev | Multiscale 3D genome rewiring during mouse neural development[END_REF]. These maps are often interpreted by determining the position along genomes of 3D structural features such as Topologically Associated Domains (TADs) boundaries, loops and compartments (e.g. in [START_REF] Suhas | A 3D map of the human genome at kilobase resolution reveals principles of chromatin looping[END_REF]). Many algorithms have been developed to delineate these features of the 3D genome organisation but their power is always hindered by the limited resolution of the dataset itself. While in principle the resolution can be improved arbitrarily by increasing the sequencing depth, this also dramatically increases the financial cost of the experiment. Improving the resolution by computational methods is therefore a good option and a methods relying on deep neural networks (HiCPlus) has recently been developed [START_REF] Zhang | Enhancing Hi-C data resolution with deep convolutional neural network HiCPlus[END_REF]. In this method, contact enhancement relies on determining more precisely low contact frequencies from the contact patterns of the genomic neighbours. Adopting a different viewpoint, our guideline is to use instead the path-length on the contact graph [START_REF] Morlot | Network concepts for analyzing 3D genome structure from chromosomal contact maps[END_REF] as a quantitative index indicating how low-confidence contact frequencies are to be reinforced (Figure 1A). The shorter the path on the contact graph between two genomic sites, the higher should be their contact frequency. In this scheme, low confidence long range contacts are determined as a sum of higher confidence contacts. The numerical implementation of this principle, Boost-HiC, is thus expected to dramatically improve the accuracy, reliability and usability of lower-resolution contact maps.

Methods

Hi-C sequence alignment

We used previously published dataset [START_REF] Bonev | Multiscale 3D genome rewiring during mouse neural development[END_REF] from mouse Embryonic Stem (ES) cells and Cortical Neurons (CN) available as GSE96107. Hi-C reads were processed using the mm9 reference genome with HiC-Pro [START_REF] Servant | HiC-Pro: an optimized and flexible pipeline for Hi-C data processing[END_REF]. Bowtie2 was used with default pipeline parameter -very-sensitive -L 30 -score-min L,-0.6,-0.2 -end-to-end -reorder. We analysed replicates for each dataset at a resolution of 10kb separately before merging them at the end of the HiC-Pro pipeline. This merged array defines the raw contact map, M .

Hi-C Filtering

Hi-C maps were filtered in two steps. First, empty bins are removed. In a second step, the distribution of the sum of contacts for each bin is fitted using a Gaussian Kernel Density function [START_REF] Pedregosa | Scikit-learn: Machine learning in Python[END_REF] with parameter kernel='gaussian', bandwidth=2000. Since Gaussian Kernel Density function returns the logarithm of the distribution, we took the exponential of the output. The resulting distribution is then used to identify bins for which the sum of contacts has a value below 5 % and above 95 % of the mean contact number. These bins are removed from the contact map. Filtered datasets were finally normalised using the Sequential Component Normalization [START_REF] Cournac | Normalization of a chromosomal contact map[END_REF] so that the L1-norm of each column/line is equal to 1.

Contact Probability

Contact probability P (s) as a function of the genomic distance s is computed as the mean number of contacts along each secondary diagonal, where s is its distance to the main diagonal. Since there are less contacts between bins which are more distant than between closer ones, contacts are gathered over bins with a size increasing according to a geometric progression (log-binned) of scale factor 1.01. When necessary, we denote P C (s) the contact probability curve associated with the contact map C.

2/10

Boost-HiC Implementation

Boost-HiC algorithm operates in two steps. The normalised contact map C is first transformed element-wise in a distance matrix (d o ) as described in [START_REF] Lesne | 3D genome reconstruction from chromosomal contacts[END_REF] according to the formula : [START_REF] Robert | Algorithm 97: shortest path[END_REF][START_REF] Warshall | A theorem on Boolean matrices[END_REF] is then run on d o to fill the distance matrix d, by defining the distance d ij between any pair of loci i and j as the length of the shortest-path connecting them on the weighted contact graph [START_REF] Morlot | Network concepts for analyzing 3D genome structure from chromosomal contact maps[END_REF]. This complete distance matrix d is then turned element-wise into a filled contact map F o according to the formula :

d o ij = (C ij ) -α A shortest-path algorithm
F o ij = (d ij ) -1/α
The number of rewired contacts (nr) is computed from C and F o by counting the number of elements (ij) taking different values in these two maps. We started with a small value (e.g. 0.05) of the exponent α and then increased this value by step of 0.01 and computed nr. When nr started to increase, we stopped the procedure and kept this value for the exponent α.

We adjusted the contact probability curve associated with F o so that it matches exactly the one computed for C. Each contact F o ij for sites i and j separated by a genomic distance s = |i -j| is adjusted by the ratio between the raw probability against the improved one :

F ′ ij = F o ij P C (|i -j|)/P F o (|i -j|)
The code recapitulating Boost-HiC steps for a given value of alpha is available at https://github.com/LeopoldC/Boost-HiC.

HiCPlus Implementation

HiCplus as been used as described in [START_REF] Zhang | Enhancing Hi-C data resolution with deep convolutional neural network HiCPlus[END_REF]. We modified the original code in order to apply the algorithm of the whole contact map whereas it was originally developped to enhance short range contacts only. The version of the code we used is available at https://github.com/jbmorlot/HiCPlus/.

Downsampling

In order to quantify the efficiency of our algorithm, we wanted to see if it is capable of inferring the same information as a deeply sequenced map when starting from a map with less contacts. To implement this test, we sampled the raw map with a binomial probability. For every couple of bins i and j in the raw map, the number of contacts M ij between these two genomic regions is replaced with a random number generated from a binomial distribution of parameters M ij and k. We downsampled the matrix with 3 different values of k : 0.1, 0.01 and 0.001 corresponding respectively to the downsampling of the whole map at 10 %, 1 % and 0.1 %.

Compartments and TADs

Compartments were obtained as in [START_REF] Nynke | Hi-C: a method to study the three-dimensional architecture of genomes[END_REF]. First, normalised contact maps were transformed to their observed/expected ratio. Correlation matrices of these ratios, also termed correlation maps, were then computed. The sign of the first eigenvector of these correlation maps was used to infer two compartments. A (active) and B (inactive) compartments were then identified by computing the gene density in each compartment and assigning the label A to the highest-density compartment. To infer the position of Topologically Associating Domains (TADs), we computed a "TAD border score" for each 10 kb position p along the genome. This score obtained by summing the contact map elements over a square region of size h with its corner located on the main diagonal, at the poistion p. We next plotted the score along the genome by sliding the square along the diagonal of the contact map [START_REF] Kruse | TADtool: visual parameter identification for TAD-calling algorithms[END_REF]. We chose here a size h of 300kb. Local minima of this score determined the TAD borders [START_REF] Smith | Invariant TAD boundaries constrain cell-type-specific looping interactions between promoters and distal elements around the CFTR locus[END_REF].

3/10

Results and Discussion

Boost-HiC partially restores the resolution of downsampled contact maps

Hi-C datasets are generally processed in the form of contact maps, in which the contact frequency between two loci is determined by normalising the number of times they have been found together in the Hi-C library. At fine resolution (i.e. when 4/10 contacts are aggregated over small genomic bins), this number is frequently equal to zero due to the finite sequencing depth. The issue is currently circumvented by working at a lower resolution, by binning the restriction fragments into larger and less numerous regions. We propose a computational alternative, Boost-HiC, to infer the missing fine-resolution contact frequencies from the knowledge of the measured contacts, and get a complete fine-resolution contact map. The procedure is based on the computation of shortest-paths between any pair of genomic loci on the contact network, a method that we introduced previously for 3D reconstruction [START_REF] Lesne | 3D genome reconstruction from chromosomal contacts[END_REF][START_REF] Morlot | Network concepts for analyzing 3D genome structure from chromosomal contact maps[END_REF]. More precisely, we first compute the dependency of the contact probability on a normalized contact map [START_REF] Cournac | Normalization of a chromosomal contact map[END_REF] C with respect to the genomic distance s between two loci, P (s). We also transform the contact map, with elements C ij , in a distance map with elements

d o ij = (C ij ) -α
where the exponent α is a parameter that will be later optimised. This initial distance map d o defines a weighted contact network, where the link between the loci i and j is given a length d o ij . We then compute an updated distance map using a shortest-path algorithm : the distance d1 ij between any two loci i and j is set to the minimal distance on the graph between nodes i and j, that is, the length of the shortest path relating i and j. In this way we get a full distance map d1 with finite elements d ij , even when the measured contact frequency between the loci i and j vanishes. This distance map is then converted into a contact frequency map F o by inverting the above relation element-wise, namely F o ij = (d ij ) -1/α . Since short-range contact values were often used to infer longer-range contact values, the P (s) curve may change after the procedure. A proper element-wise rescaling of F o restores the original contact probability curve P (s) at this step (Figure 1B). The resulting map F ′ is finally normalised using the Sequential Component Normalization procedure [START_REF] Cournac | Normalization of a chromosomal contact map[END_REF] to give the final 'boosted' contact map F . This map F lies at the same fine resolution than the original one, but vanishing elements are replaced with their inferred values.

To assess the capabilities of our algorithm, we generated 3 downsampled contact maps sampling 10 %, 1 % and 0.1 of the initial contacts from a high-resolution map of region of chromosome 16 in mouse at 10kb resolution (see Methods). We then constructed the corresponding boosted maps (Figure 1 C-F). In order to quantitatively compare the boosted map and our objective (i.e. the original high resolution version), we computed the Spearman correlation between the boosted maps and the high-resolution map for different values of our free parameter α. While some other methods have been developed to compare contact maps from different origin [START_REF] Michael Eg Sauria | HiFive: a tool suite for easy and efficient HiC and 5C data analysis[END_REF][START_REF] Yan | HiC-spector: a matrix library for spectral and reproducibility analysis of Hi-C contact maps[END_REF][START_REF] Yang | HiCRep: assessing the reproducibility of Hi-C data using a stratum-adjusted correlation coefficient[END_REF], this straightforward measure works well when comparing downsampled maps with the original one. We found that the downsampled maps displayed a reduced correlation of 0.51, 0.26 and 0.11 for 10 % ,1 % and 0.1 % downsampling respectively. When α is optimally chosen, this correlation increases to 0.59, 0.51, 0.47 for the boosted maps obtained for the corresponding downsampled maps. In comparison, HiCPlus [START_REF] Zhang | Enhancing Hi-C data resolution with deep convolutional neural network HiCPlus[END_REF], which uses a deep-learning approach where contacts are enhanced using the information contained in the contacts established between the adjacent sites along the genome, achieves correlation values of 0.59, 0.44 and 0.34. This latter result could be explained by the fact that the enhancing procedure HiCPlus was specifically designed for short-range contacts whereas Boost-HiC mostly improves low, long-range contact frequencies.

Optimising the parameter α

The different values of the Spearman correlation obtained for increasing values of our free parameter α showed that unless α is chosen below 0.1 or above 1, the choice of its 5/10 value has a mild effect on the reconstruction efficiency (e.g. it changes from 0.53 to 0.58 for a down-sampling of 10 %, Figure 1 E)). The best α (i.e. yielding to the most accurate reconstruction of the full initial map) is equal to 0.25, 0.5 and 0.6 for different downsampling rates and it is a priori difficult to know which value to use in a practical case. For low α values, only vanishing elements of the sparse original matrix are changed into non-zero values. When α increases, the number of elements of the matrix that are updated in the Boost-HiC procedure increases and non-zero elements are also modified. In order to see whether the optimal value for α depends on the number of modified values, we plotted the number of non-zeros values that are reassigned as a function of the downsampling rate and α value. As expected, as the downsampling rate increases, more and more elements become equal to zero in the downsampled matrix. Interestingly, we showed that for all three downsampling rates, the optimal value of α corresponds to the updating of vanishing (i.e. all the zeros and few non-zeros) elements only (Figure 1 F). Importantly, this can give us a way to choose an optimal α value in real cases, where we have only access to sparse data. We therefore implemented a search of the optimal α, which corresponds to a re-assignment of all the zeros and almost 10 % of the non-zeros elements only.

Since the alpha optimisation step can be long for large matrices, we propose in this case to estimate the optimal alpha only based on the sparsity of the matrix. For low sparsity data (i.e. lower than 0.9) an alpha value of 0.25 works well whereas for high sparsity data (i.e. higher that 0.99) alpha should be chosen close to 0.6. For intermediate sparsity, alpha should ideally be taken between those two bounds.

Boost-HiC enables the precise determination of compartments from low resolution maps

As the contact patterns vanish in down-sampled maps (Figure 1 C-F), the detection of A/B compartments using the state-of-the-art procedure is gradually impaired. At 0.1 % down sampling compartments cannot be identified anymore (Figure 2 A). When the Boost-HiC procedure is applied to the down-sampled maps, the detection of these structural features is restored, showing the ability of our computational strategy to reliably detect 3D features from sparse data (Figure 2 B). Using HiCPlus to enhance the down-sampled maps also resulted in a partial recovery of the compartments (Figure 2 C). All these results are summarised in Figure 2 D), which shows the bins that are attributed to the good compartment, as determined from the high resolution map, for different down-sampling values and enhancement methods. Boost Hi-C robustly allows a better prediction of compartments from low resolution maps. On the other hand, since TADs are determined from short range contacts, their detection in not changed by the Boost-HiC procedure (Figure 2 E). Indeed the TAD border signal is not (or very mildly) modified by the Boost Hi-C procedure. On the other hand, we noted that HiCPlus does change the signal around the diagonal, although the position of TADs borders as minima of the signal are not modified.

Boost Hi-C enables high-resolution comparison between contact maps

We finally assessed the performance of the Boost-HiC procedure for comparing contact maps obtained in different biological conditions. The comparison is usually done by computing the element-wise log-ratio of two such maps. In order to be informative, this comparison can only be done at a resolution for which the signal-to-noise ratio is high. On Figure 3, we compare two contacts maps from two different cell types at different sequencing depths, with and without using Boost-HiC. We compared here 6/10 Even if these maps only represent a sub-region of mouse chromosome 16, many pairwise contacts between 10kb bins are equal to zero in a typical Hi-C experiment, even at the highest coverage (Figure 3 A). When computing the log ratio, these elements of the ratio map will therefore be undefined, as they will be equal to +/-infinity (in green and yellow on Figure 3, right panels). In contrast, when the maps have been complemented with Boost-HiC prior to this differential analysis, these elements will always be filled with a finite value (in red and blue on Figure 3, right panels). As can be expected, this effect is even more pronounced at low resolution (e.g. at 10 % downsampling in Figure 3B)

Conclusion

We present here Boost-HiC, a computational method aiming at enhancing low resolution chromosomal contact maps. We showed that the resulting maps can be efficiently used in order to recover chromosomal compartments at high resolution. We also show that the procedure can also be efficiently used when comparing Hi-C contact maps when looking for differential contacts between two conditions. 8/10

1 Figure 1 :

 11 Figure1: A) Schematic explanation of our shortest-path enhancement algorithm on Hi-C map : the low-confidence long-range contact between sites a and b (white dot) is enhanced using the information on shorter-range higher-confidence contacts (black dots) with intermediary nodes (here c). On the contact graph, the shortest-path between a and b goes through c and the enhanced contact value is related to its length. B) Contact probability curve P (s) of ES-cells from their Hi-C map downsampled at 10 % of contacts. The blue points display the original probability of the downsampled map. Yellow line displays the probability after information enhancement by means of our shortest-path method. The red line displays the improved contact probability at the end of Boost-HiC algorithm. C) to F) panels display the mouse ES-cell contact map (for a region on chromosome 16 located between from 100kb and 29,8Mb) at different downsampling rates (in order) : 100 %, 10 %, 1 % and 0.1 %. Upper parts displays the raw contact map at fixed downsampling. Lower parts shows the map improvement after implementation of Boost-HiC. Colormap in log10 scale[START_REF] Cournac | Normalization of a chromosomal contact map[END_REF]. G) Pearson correlation (in red) between full raw maps and boosted maps for different downsampling rates and value of the parameter α used in the Boost-HiC algorithm. The blue line shows the fraction of rewired contacts between the downsampled map and the boosted map.From top to bottom, we compared the 100 % map with respectively the 10 %, 1 % and 0.1 % map enhanced with Boost-Hi-C.

FFigure 2 :

 2 Figure2: A) Barplot of the first eigenvector of mouse ES-cell correlation map (for a region on chromosome 16 located between 100kb and 29,8Mb) at different downsampling rates (in order from top to bottom) : 100 %, 10 %, 1 %, 0.1 %. B) Similar barplot computed after application of the Boost-HiC algorithm. C) Similar barplot computed after application of HiCPlus algorithm[START_REF] Zhang | Enhancing Hi-C data resolution with deep convolutional neural network HiCPlus[END_REF] D) Fraction of regions attributed to the correct compartment from downsampled maps using different enrichment methods for 10 %, 1 % and 0.1 % downsampling. E) TAD border signal along the genome coordinate (from 100kb to 29,8Mb of chr16). In all plots the blue line shows the score for the raw normalised dataset at 100%. In the top plot red and yellow line correspond to the dataset downsampled at 10% enhanced respectively with Boost-HiC and HiCPlus methods. F) Same plot as E) with Boost-HiC and HiCPlus methods used on the dataset downsampled at 1%.

Figure 3 :

 3 Figure 3: Log 2 -ratio between CN and ES-cell contact maps, before (left panels) and after (right panels) application of the Boost-HiC method. An null value in one dataset will give an infinite value in the log 2 -ratio. A) Log 2 -ratio for maps without downsampling. B) Log 2 -ratio for maps downsampled at 10 %.

Funding

This work has been funded by the programme InFIniTI 2018 of the Mission for Interdisciplinarity of the French Centre National de la Recherche Scientifique (CNRS), grant 238301 (to A.L.) and by the Agence Nationale pour la Recherche (HiResBac ANR-15-CE11-0023-03 to J.M.),


P Cordesse

A Murrone

T Menard

M Massot

Initiation of a validation strategy of reduced-order two-fluid flow models using direct numerical simulations in the context of jet atomization

In industrial applications, developing predictive tools relying on numerical simulations using reduced-order models nourish the need of building a validation strategy. In the context of cryogenic atomization, we propose to build a hierarchy of direct numerical simulation test cases to assess qualitatively and quantitatively diffuse interface models. The present work proposes an initiation of the validation strategy with an air-assisted water atomization using a coaxial injector.

Introduction

Rocket engines reliability is one of the main priorities given to the European Ariane 6 program. Cryogenic combustion chambers encompass several interacting physical phenomena covering a large spectrum of scales. In particular, the primary atomization plays a crucial role in the way the engine works, thus must be thoroughly studied to understand its impact on high frequencies instabilities. In sub-critical condition, the two-phase flow topologies encountered in the chamber vary from two separated phases at the exit of the nozzle to a polydisperse spray of droplets downstream. In between, liquid structure such as ligaments, rings and deformed droplets detach from the liquid core thus the subscale physics and the topology of the flow in this mixed region is very complex. Predictive numerical simulations are mandatory, at least as a complementary tool to experiments to understand the physics and eventually to conceive new combustion chambers and predict instabilities they may generate.

Since direct numerical simulations of these two-phase flows in realistic configurations are still out of reach, computational power needs being too high, reduced-order models are deployed. However great care must be taken on the choices of these models to combine sound mathematics properties with satisfying representativeness.

Two approaches are found in the literature to build reduced-order models for jet atomization: 1. couple two models, a first one for the dispersed flow by employing an element derived from the Kinetic Based Moment Method [START_REF] Sibra | Simulation of reactive polydisperse sprays strongly coupled to unsteady flows in solid rocket motors: Efficient strategy using Eulerian Multi-Fluid methods[END_REF] which treats polydispersed droplets in size, velocity and temperature, and a second one dealing with the separated phase and the mixed region with either a hierarchy of diffuse interface models [START_REF] Drui | Small-scale kinematics of two-phase flows: identifying relaxation processes in separated-and disperse-phase flow models[END_REF] or some level of Large Eddy Simulations (LES) on the interface dynamics [START_REF] Herrmann | A sub-grid surface dynamics model for sub-filter surface tension induced interface dynamics[END_REF]. 2. use a unified model that encompasses any flow topology. Works in this direction are found in [START_REF] Drui | Small-scale kinematics of two-phase flows: identifying relaxation processes in separated-and disperse-phase flow models[END_REF] where a unified model accounting for micro-inertia and micro-viscosity associated to bubble pulsation is proposed.

In the present study, we focus on the first approach. We eventually want to perform a numerical simulation of the primary atomization by using a Kinetic-Based-Moment Methods (KBMM) modelling the dispersed flow as in [START_REF] Sibra | Simulation of reactive polydisperse sprays strongly coupled to unsteady flows in solid rocket motors: Efficient strategy using Eulerian Multi-Fluid methods[END_REF] coupled with a hierarchy of diffuse interface models describing the separated phases and the mixed zone. In order to qualitatively and quantitatively assess the diffuse interface models, we propose a validation strategy based on a hierarchy of direct numerical simulation test cases.

2 Mathematical modelling 2.1 Reduced order methods -diffuse interface model Among the hierarchy of the diffuse interface models, the Baer and Nunziato model first introduced in [START_REF] Baer | A two-phase mixture theory for the Deflagration-to-Detonation Transition (DDT) in reactive granular materials[END_REF] has been generalized in [START_REF] Saurel | A Multiphase Godunov Method for Compressible Multifluid and Multiphase Flows[END_REF] by introducing the interfacial quantities p I and v I , the interfacial pressure and the interfacial velocity respectively, to model the dynamic of the interface. The system of equation takes the form in one dimension:

∂ t q + {∂ q f (q) + N (q)} ∂ x q = r(q) with ∂ q f (q) =   0 0 0 0 ∂ q 2 f 2 (q 2 ) 0 0 0 ∂ q 1 f 1 (q 1 )   , N (q) =   v I 0 0 n 2 0 0 n 1 0 0   (1)
with q ∈ R 7 , the state vector defined by q = (α 2 , q 2 , q 1 ) t ,

q k = (α k ρ k , α k ρ k v k , α k ρ k E k ), f (q) ∈ R 7 the conservative flux f (q) = (0, f 2 (q 2 ), f 1 (q 1 )) t with f k (q k ) = (α k ρ k v k , α k ρ k v 2 k + α k p k , α k (ρ k E k + p k )v k ) t , N (f ) ∈ RR 7×7 the non-conservative terms with n 2 (q) = -n 1 (q) = (0, -p I , -p I v I ) t , α k the volume fraction of phase k ∈ {1, 2}, ρ k the partial density, v k the phase velocity, p k the phase pressure, E k = e k + v 2
k /2 the total energy per unit of mass, e k the internal energy, r(q)/ ∈ R 7 the relaxation terms detailed in the sequel.

Mathematical properties of System (1) have been widely examined in the literature [START_REF] Embid | Mathematical analysis of a two-phase continuum mixture theory[END_REF][START_REF] Coquel | Closure laws for a two-fluid two-pressure model[END_REF][START_REF] Gallouët | Numerical modeling of two-phase flows using the two-fluid two-pressure approach[END_REF]: it is a first-order non-linear system of partial differential equations which include conservative and non-conservative terms; the system is conditionally hyperbolic and admits seven eigenvalues,

v I , {v k , v k ± a k } k=2,1 (2)
with a k the phase sound speed defined by

a 2 k = ∂p k ∂ρ k s k .
The interfacial terms need a closure and the thermodynamics must be postulated. A recent work [START_REF] Cordesse | Supplementary conservative law for non-linear systems of PDEs with nonconservative terms: application to the modelling and analysis of complex fluid flows using computer algebra[END_REF] has extended the theory of first-and second-order non-linear conservative systems [START_REF] Godunov | An Interesting Class of Quasilinear Systems[END_REF][START_REF] Kawashima | On the normal form of the symmetric hyperbolic-parabolic systems associated with the conservation laws[END_REF] to non-conservative systems, and has applied it to the Baer-Nunziato model to obtain a compatible thermodynamics and model closure. All the existing closures in the literature are recovered.

The interfacial velocity v I is commonly chosen in order to be linearly-degenerated which restricts its definition to

v I = βv 1 + (1 -β)v 2 with β ∈ {0, 1, α 1 ρ 1 /ρ} [7].
For the present work, the phases are assumed immiscible such that the postulated entropy H takes the form:

H = - k=1,2 α k ρ k s k (3)
with s k = s k (ρ k , p k ) the phase entropy following a two-parameter equation of state.

The relaxation of the two phases towards an equilibrium state is described by the sources terms r/ in System (1). Only mechanical, and hydrodynamic relaxations are accounted for in the present study and therefore r/ decomposes into:

r = r v u + r p p , with r v u = 0, r v 2 u , r v 1 u t and r p p = p 2 -p 1 p , r p 2 p , r p 1 p t (4)
where

r v 2 = -r v 1 = (0, v 2 -v 1 , v I (v 2 -v 1 )), r p 2 = -r p 1 = (0, 0, p I (p 2 -p 1 )
) and p (resp. u ) is the characteristic time of the mechanical relaxation (resp. hydrodynamical relaxation). When relaxing also the temperatures, one obtains the compressible multi-species Navier-Stokes equations, called also four equation model. The Baer-Nunziato model introduced herein before is thus called the seven equation model.

Incompressible Navier-Stokes equations and VOF for DNS

In order to describe two-phase flows, the incompressible Navier-Stokes equations are introduced in System ( 5)

   ∇ • v = 0 ∂ t v = -(v • ∇) v + 1 ρ (-∇p + ∇ • (2µT) + F) (5)
where v ∈ RR3 is the velocity field, µ the dynamic viscosity, ρ the mixture density, p the equilibrium pressure, T the strain tensor defined by T = 1 2 ∇(v) + ∇(v T ) , F = F V + F ST represents the body force and the surface tension force, F ST = γκδ I n, γ is the surface tension, κ the curvature at the interface, δ I the Dirac function and n the interface outwards normal vector.

Then, the interface is implicitly derived by the zero of a Level Set of a scalar function φ which stands for the signed distance to the interface (|∇φ| = 1, φ > 0 in the liquid and φ < 0 in the gas). Geometrical properties of the interface, normal and curvature, are easily deduced, n = ∇φ/|∇φ| and κ = ∇ • n. The interface motion is captured by the transport of the Level Set function:

∂ t φ + v • ∇ φ = 0 (6)
However, the property of distance function is not generally conserved after the resolution of Equation ( 6), an additional equation have to be solved [START_REF] Sussman | A Level Set Approach for Computing Solutions to Incompressible Two-Phase Flow[END_REF]. As a consequence of solving both equations, the conservation of mass is not guaranteed. To solve this problem, the volume fraction of liquid (VOF -Volume Of Fluid) α l in a volume Ω is introduced in terms of the Level Set function φ(x) as:

α l = 1 Ω Ω H(φ(x))dΩ (7)
and the transport of this quantity in the control volume Ω is given by :

∂ t α l + v • ∇ α l = 0 (8)
The benefits of a VOF formulation coupled with a Level Set function are to conserve mass and to have access to geometrical properties of the interface.

3 Numerical methods

Diffuse interface model

The numerical methods employed to solve System (1) are implemented in the multiphysics computational fluid dynamics software CEDRE [START_REF] Gaillard | Numerical Simulation of Cryogenic Injection in Rocket Engine Combustion Chambers[END_REF] working on general unstructured meshes and organized as a set of solver [START_REF] Gaillard | Numerical Simulation of Cryogenic Injection in Rocket Engine Combustion Chambers[END_REF]. The solver SEQUOIA is in charge of the diffuse interface model.

A Strang splitting technique is applied on a multi-slope HLLC with hybrid limiter solver [START_REF] Furfaro | A simple HLLC-type Riemann solver for compressible non-equilibrium two-phase flows[END_REF][START_REF] Le Touze | Multislope MUSCL method for general unstructured meshes[END_REF] to achieve a time-space second-order accuracy on the discretized equations. The issue of the non-conservative terms encountered when solving System (1) is tackled in [START_REF] Furfaro | A simple HLLC-type Riemann solver for compressible non-equilibrium two-phase flows[END_REF] by assuming (1) the interfacial terms p I and v I to be local constants in the Riemann problem, (2) the volume fraction to vary only across the interfacial contact discontinuity v I . As a result, the non conservative terms in System (1) vanish, v I and p I are determined locally by Discrete Equation Method (DEM) [START_REF] Saurel | A multiphase model with internal degrees of freedom : application to shock-bubble interaction[END_REF] at each time step and stay constant during the update. Thus, the phases are decoupled and the System (1) splits into two conservative sub-systems onto which we apply the multi-slope HLLC with hybrid limiter solver.

Depending on the application, the relaxations are assumed either instantaneous or finite in time. In our case, it is reasonable to assume an instantaneous pressure relaxation but to consider a finite velocity relaxation as in the context of an assisted air-water coaxial injector, the interface dynamic is mainly driven by the shear stress induced by a high velocity gradient at the injection.

To obtain the relaxed pressure, one needs to solve the ODE,

∂ t q = r p (q) p , (9)
with p → ∞ which infers v k remains constant. The problem reduces to apply an iterative procedure as a Newton method to solve a second order equation. Details of the equation can be found in [START_REF] Furfaro | A simple HLLC-type Riemann solver for compressible non-equilibrium two-phase flows[END_REF]. As for the velocities, the following ODE is solved

∂ t v d - A o v d u = 0, ( 10
)
where v d is the slip velocity v d = v 2 -v 1 and superscript o denotes the state before relaxation. A first numerical approach is to fix a remaining slip velocity ratio target at each computational time step ∆t. It defines the characteristic relaxing time:

u A o = ln (X)∆t with X = v d (∆t) v o d and A o = α o 1 ρ o 1 + α o 2 ρ o 2 α o 1 ρ o 1 α o 2 ρ o 2 (11)
An instantaneous velocity relaxation is in pratice also possible and manipulating the ODE leads to a unique relaxed velocity, which is the mass weighted average of the two velocities before relaxing.

DNS numerical methods

As mentioned by Rudman [START_REF] Rudman | A volume-tracking method for incompressible multifluid flows with large density variations[END_REF], it is recommended to solve Navier-Stokes equations in their conservative form to ensure consistance between mass fluxes (equivalent to VOF fluxes) and momentum fluxes. Since a staggered grid is adopted, VOF and velocity do not have the same control volume and it is difficult to obtain this consistency. This is the reason why Rudman introduces a grid two times smaller in each direction for VOF transport. In [START_REF] Vaudor | A consistent mass and momentum flux computation method for two phase flows. Application to atomization process[END_REF], Vaudor developed a method to avoid this finest grid. This method is used here and allows to reduce computation time.

The Navier-Stokes equations are implemented in the code ARCHER [START_REF] Vaudor | A consistent mass and momentum flux computation method for two phase flows. Application to atomization process[END_REF]. A projection method is employed and physical properties (viscosity and density) are expressed in term of both VOF and Level Set. The temporal integration is performed through a second-order Runge-Kutta scheme. The discretization of convective term is achieved with WENO 5 scheme [START_REF] Jiang | Efficient Implementation of Weighted {ENO} Schemes[END_REF]. For the viscosity term, we retain the method presented by Sussman [START_REF] Sussman | A sharp interface method for incompressible two-phase flows[END_REF]. The Ghost-Fluid [START_REF] Fedkiw | A Non-oscillatory Eulerian Approach to Interfaces in Multimaterial Flows (the Ghost Fluid Method)[END_REF] method is employed to take into account surface tension force F ST , treated as a pressure jump.

To compute correctly mass fluxes at the interface (consequently momentum fluxes) and consistence between both flow and interface solvers, VOF and Level set transport (Equations ( 6) and ( 8)) are performed with a CLSVOF algorithm [START_REF] Sussman | A Coupled Level Set and Volume-of-Fluid Method for Computing 3D and Axisymmetric Incompressible Two-Phase Flows[END_REF] at each step of Runge Kutta scheme.

Results and discussions

As part of the validation strategy, the hierarchy of direct numerical simulation test cases starts with an air-assisted water atomization with a coaxial injector and results are presented hereafter.

The injector at Geophisic and industrial flow laboratory (LEGI) has been the subject of several experiments [START_REF] Rehab | Flow regimes of large-velocity-ratio coaxial jets[END_REF][START_REF] Lasheras | Liquid Jet Instability and Atomization in a Coaxial Gas Stream[END_REF][START_REF] Marmottant | Atomisation primaire dans les jets coaxiaux[END_REF][START_REF] Delon | Flapping instability of a liquid jet[END_REF] covering a large range of flow conditions and thus offers experimental data to assess numerical simulations. The present simulations reproduce the experiments conducted by [START_REF] Delon | Flapping instability of a liquid jet[END_REF], whose results of interests for our study are presented in [START_REF] Vaudor | A consistent mass and momentum flux computation method for two phase flows. Application to atomization process[END_REF] along with the numerical results obtained with ARCHER.

Description of the LEGI experiment

Figure 1 renders the simulation configuration which is identical to the experiment of [START_REF] Delon | Flapping instability of a liquid jet[END_REF] and indicates the velocity profiles at the injector outlets measured experimentally. The inside diameter measures d l = 7, 6 mm, the outside diameter d g = 11, 4 mm, the lip length e = 0.2 mm. The gas velocity profile V g given in Equation ( 13) models the boundary layer measured experimentally with δ the boundary layer thickness and V max g the maximum gas velocity.

V l Vg δ rg r l r l +e r x V l = 0.26 m.s -1 (12) V g =                          r -r l -e δ V max g r l + e ≤ r and r ≤ r l + e + δ V max g r l + e + δ ≤ r and r ≤ r g -δ 1 - r -(r g -δ) δ V max g r g -δ ≤
The fluid properties, type, density ρ, capillarity coefficient γ and viscosity coefficient µ of each fluid are summarized in Table 1.

Phase ρ (kg.m -3 ) γ (N.m -1 ) µ (1e

l = ρ L V L d l µ L = 1972, Re g = ρ G V G (d g -d l -2e) µ G = 5854, M = ρ G V 2 G ρ L V 2 L = 11, (14a)
We

l = ρ L V 2 L d l γ = 7, We g = ρ G V 2 G (d g -d l -2e) γ = 36, We = ρ G V 2 G d l γ = 81. ( 14b
)
The Weber numbers are small inferring the surface tension should not play a crucial role in the dynamics of the flow. For the present configuration, there exists a break-up regime map established by [START_REF] Lasheras | Liquid Jet Instability and Atomization in a Coaxial Gas Stream[END_REF], reported in Figure 2. The red dot of coordinates (Re l = 1972, We = 81) locates the investigated flow in the shear w Fig. 2: Break-up regimes in the parameter spaces Re l -We [24] -LEGI breakup zone in Figure 2. We should not expect atomization of the liquid jet and therefore the atomization is not activated in CEDRE.

Numerical set-up

The ARCHER simulations are performed on a Cartesian mesh 1024 × 512 × 512 with a cell size equal to ∆x = 6.68 10 -5 , so a total of 806 millions of faces whereas CEDRE mesh is composed of 278 978 tetrahedral cells. The mesh size is 563 116 faces.

To be comparable to the DNS, two difficulties must be tackled. Firstly, the DNS uses an incompressible solver meaning the acoustic is not solved and thus does not interact with the liquid jet and the density remains constant. Secondly, the boundary conditions imposed on the wall of the DNS box do not let any flow backwards such that no wall effect acts upon the liquid core. Therefore, to eliminate reflected acoustic waves and wall effect on the liquid jet, we have designed an outer box with a coarsening mesh as shown in Figure 3. The smallest cell is located at the lip of the injector and measures ∆x = 2.0 10 -4 m. As for the thermodynamics, CEDRE uses two Stiffened-Gas equations of state and thus the temperature of the phases has been modified to obtain the same initial pressure and density conditions as in Table 1. time, the total CPU is 150 greater for ARCHER than for CEDRE.

Qualitative comparison

We propose to compare the liquid jet obtained with ARCHER and CEDRE at given simulation time. The Level Set function of the DNS permits an exact reconstruction of the interface whereas for the diffuse interface model, the interface lays in the region where the volume fraction varies from α ≈ 0 to α ≈ 1. Consequently, in Figure 4 we have superimposed the solved interface of the DNS to a volume rendering of the liquid volume fraction and a single liquid volume fraction isosurface. We distinguish two regions: close to the injector, on the first half of the DNS box, the diffuse interface model is able to match accordingly to the DNS results. In this region, the mesh used by CEDRE prevents the interface to diffuse too much for the diffuse interface model and the interface of the DNS lays in the volume rendering of the liquid volume fraction of CEDRE.

On the second half of the DNS box, where the DNS shows complex liquid structures such as rings, droplets and ligaments, the diffuse interface model is not able to capture these effects since the mesh is coarsened and the volume fraction alone is not sufficient to describe such complex interface dynamics. Nevertheless, the seven equation model diffusion accords well with the DNS, the red cloud which corresponds to a low liquid volume fraction is limited to the zones where liquid elements of the DNS exists.

This comparison gives an interesting interpretation of the diffuse interface models. The volume fraction is not enough to reconstruct the whole dynamic of the interface but attests the presence or the absence of liquid.

Quantitative comparison

The liquid core is defined as the region of the liquid jet that is always occupied by liquid. To obtain it, a time-averaging of the liquid volume fraction is needed. To compare the liquid core obtained with CEDRE, L lc (α 2 = 0.99) for t > 60 ms.

For t > 60 ms the liquid jet length starts to oscillate around a time averaged value which equals L lc (α 2 = 0.99 = 23 mm. We then computed the time average liquid volume fraction between t sim ∈ [60, 160] ms to obtain the liquid core shown in Figure 6. The liquid core of the diffuse interface model was identified as the isovolume of α l = 1 -2 , where = 1e -6 is the residual volume fraction. The reason for not choosing α l = 1holds to the fact that due to numerical diffusion in a unstructured mesh and mesh interpolation, the volume fraction in single phase region does not stay at the initial value α k = 1 -.

Experimental, DNS and CEDRE values show the same trend. We also could have compared the angle of the spray, but for the diffuse interface model, the result depends highly on the threshold chosen for the liquid volume fraction.

Conclusion

In the present contribution, we started a validation process of reduced-order models in the context of cryogenic propulsion using direct numerical simulations.

The test case proposed here is an air-assisted water atomization using a coaxial injector which also provides experimental results from the LEGI test bench. The comparison showed good agreements and important CPU gains between the seven equation model implemented in the CEDRE code and the DNS results from the ARCHER code.

It showed the limits of diffuse interface models to capture complex liquid structures such as ligaments, rings or deformed droplets and encourages to add a sub-scale description of the interface dynamics through geometric variables such as the interfacial area density, the mean and Gaussian curvatures as proposed in [START_REF] Cordesse | Derivation of a two-phase flow model with twoscale kinematics and surface tension by means of variational calculus[END_REF].

1 ( 13 )Fig. 1 :

 1131 Fig. 1: Injector schematic and velocity profiles.

(a ) Fig. 3 :

 )3 Fig. 3: CEDRE mesh of the configuration

Fig. 4 :

 4 Fig. 4: Instantaneous liquid core comparison. CEDRE: volume rendering of the liquid volume fraction, α l high low, grey isosurface α l = 0.99. ARCHER: Grey isovolume of liquid volume fraction α l = 1

Figure 5

 5 plots the instantaneous liquid jet length L lc over time for several α l .

Fig. 5 :

 5 Fig. 5: Instantaneous liquid jet length L lc over time at isovalue α 2 ∈ [0.91, 1 -2 ], isovalue α 2 = 0.99, α 2 = 1 -2 (bottom) and α 2 = 0.91 (top), L lc (α 2 = 0.99) for t > 60 ms.

ExpeFig. 6 :

 6 Fig. 6: Liquid core comparison -from left to right: experiment, DNS, CEDRE and DNS. CEDRE liquid core is the isovolume at α l = 1 -2 , with α l the time averaged liquid volume fraction and = 1e-6 residual volume fraction.

  -5 P a.s) Physical properties of water and air.Let us define and compute the following flow parameters. Re l is the liquid Reynolds number, Re g is the gas Reynolds number, M the momentum flux ratio, We l the liquid Weber Number, We g the gas Weber Number and We the aerodynamic Weber number. They are defined as follows Re

	Liquid (l) H 2 O	1000	0.0072	1002
	Gas (g)	Air	1.226	0.0072	17.8
		Tab. 1:

  Table (2). Remarkably, when considering the same simulation

		Tab. 2: Simulation costs
		Simu Time (s) N proc Total CPU (h)
	CEDRE	0.160	420	2.14 10 4
	ARCHER	0.500	8192	10 10 6

Acknowledgment

The support of CNES and ONERA through a PhD grant for P.Cordesse, the support of the Hadamard Doctoral School of Mathematics (EDMH), the help of L.H. Dorey (ONERA) and M. Théron (CNES) are gratefully acknowledged.