Fabienne Comte

Adeline Samson

Julien J Stirnemann

HAZARD ESTIMATION FOR CENSORED DATA CONTAMINATED WITH ADDITIVE MEASUREMENT ERROR: APPLICATION TO LENGTH OF PREGNANCY

Keywords: Censored data, Measurement error, Survival function estimation, Hazard rate function estimation, Nonparametric methods, Deconvolution

published or not. The documents may come

Hazard estimation for censored data contaminated with additive measurement error: application to length of pregnancy

Fabienne Comte, Adeline Samson, Julien J.

Introduction

The length of pregnancy is defined as the time between the spontaneous fertilization of the oocyte and spontaneous delivery following labor. However, despite the simplicity of the definition, the estimation of the physiological length of pregnancy is a challenging problem since both the time origin and the time of onset of labor are only partly observed. Several events may occur during pregnancy and interfere with its physiological course prior to spontaneous labor, such as medical induction of labor or cesarean section, intrauterine fetal death or any fetal or maternal medical condition that leads to delivery before its "physiological" time.

Using time-to-event terminology, the time to spontaneous labor is therefore randomly rightcensored. Right-censored data, in its standard presentation, involves independent observations ((X j ∧ C j ), 1 X j ≤C j ) for j = 1, . . . , n where the variable X denotes the true time between the origin and the occurrence of the event of interest and the variable C denotes the true time between the origin and the occurrence of censoring. We classically assume that X and C are independent.

As stated above, the time of fertilization (the time origin) is ill-observed and can never be directly measured, except in the very special case of in vitro fertilization, which we do not consider here, as this case cannot be considered as physiological. In practice, in spontaneous pregnancies, the time of onset of pregnancy is predicted either from the last menstrual period, by simply adding 14 days or, better, by the sonographic measurement of the length of the embryo in the first trimester.

Because it is not directly observed and estimated from biological data, the time of fertilization is known up to an error, the magnitude of which is of several days, regardless of the method.

In this paper, we consider the problem of censored data prone to measurement error. As the time origin of pregnancy is known up to an additive (random) error, both variables X (time between the true onset and the natural childbirth) and C (time between the true onset and any censoring event) are observed up to this additive error.

Let us define the model more precisely. Let ε denote the random error variable assumed to be independent of X and C. We assume that the observations are properly classified as censored or uncensored and that both the censored and uncensored observations are measured with error:

Y j = (X j ∧ C j ) + ε j = (X j + ε j ) ∧ (C j + ε j ), j = 1, . . . , n (1) δ j = 1 X j ≤C j ,
Note that the censoring indicator δ j is unchanged by the measurement error: 1 X j ≤C j = 1 X j +ε j ≤C j +ε j . The purpose of this work is to provide a nonparametric estimator of the hazard function h X of X, based on the observations (Y j , δ j ).

Nonparametric methods have been proposed in the frameworks of both censored data and deconvolution. Regarding censoring, [START_REF] Antoniadis | Density and hazard rate estimation for right-censored data by using wavelet methods[END_REF] consider a wavelet hazard estimator which is not adaptive, whereas [START_REF] Li | On the minimax optimality of wavelet estimators with censored data[END_REF][START_REF] Li | On the block thresholding wavelet estimators with censored data[END_REF] suggests estimators based on wavelet with hard or block thresholding. Estimators based on model selection via penalization have also been proposed: [START_REF] Dohler | Adaptive estimation of hazard functions[END_REF] estimate the log-hazard function using a penalized likelihood-based criterion, Brunel andComte [2005, 2008] consider penalized contrast estimators for both the density and the hazard rate using either the Nelson-Aalen estimator of the cumulative hazard function or the Kaplan-Meier cumulative hazard estimator, [START_REF] Reynaud | Penalized projection estimators of the Aalen multiplicative intensity[END_REF] proposes a penalized projection estimator of the Aalen multiplicative intensity process with adaptive results and minimax rates and [START_REF] Akakpo | Histogram selection for possibly censored data[END_REF] consider a histogram selection for both density and hazard rate estimation.

We can also consider our estimation problem in the setting of deconvolution. Deconvolution has been widely studied in various contexts. We hereby restrain to references with a known density of the noise. Kernel estimators have been proposed by [START_REF] Stefanski | Deconvoluting kernel density estimators[END_REF], [START_REF] Fan | On the optimal rates of convergence for nonparametric deconvolution problems[END_REF], with bandwidth selection strategies [START_REF] Delaigle | Bootstrap bandwidth selection in kernel density estimation from a contaminated sample[END_REF]. Wavelet estimators [START_REF] Pensky | Adaptive wavelet estimator for nonparametric density deconvolution[END_REF]Vidakovic, 1999, Fan and[START_REF] Fan | Wavelet deconvolution. Information Theory[END_REF], and projection methods with model selection [START_REF] Comte | Penalized contrast estimator for adaptive density deconvolution[END_REF] have also been developed. A pointwise estimation method for S X has been proposed by [START_REF] Dattner | On deconvolution of distribution functions[END_REF] when the data are noisy but not censored. In this work, we estimate the hazard rate h X based on the ratio of deconvolution estimators, using the developments of [START_REF] Dattner | On deconvolution of distribution functions[END_REF] in the setting of distribution functions (with no censoring) for the denominator.

The paper is organized as follows: Section 2 studies a quotient estimator of the hazard rate h X . The estimator is illustrated with a simulation study in Section 3 and is compared to results obtained when neither measurement error nor censoring are considered. The motivating application of estimation of length of pregnancy is illustrated by an analysis of real data in Section 4. In Section 5, we discuss the alternative problem where the noise affects only the variable X, a problem that remains open. Proofs are gathered in Appendix.

Notations. We denote f U the density of a variable U . We denote S U (t) = P(U ≥ t) the survival function at point t of a random variable U , h U (t) = f U (t)/S U (t) the hazard ratio at point t and f * U the characteristic function. We denote g * (t) = e itx g(x)dx the Fourier transform of any integrable function g. For a function g : R → R, we denote g 2 = R g 2 (x)dx the L 2 norm. For two integrable and square-integrable functions g and h, we denote g h the convolution product g h(x) = g(x -u)h(u)du. For two real numbers a and b, we denote a ∧ b = min(a, b).

Censored data and measurement noise

2.1. Setting. We observe for j = 1, . . . , n

Y j = (X j ∧ C j ) + ε j , δ j = 1 X j ≤C j .
We assume that the law of the noise is known and that its characteristic function is such that

∀u ∈ R, f * ε (u) = 0.
The following assumption, which is verified by exponential or Gamma distributions for example, will be considered fulfilled throughout this section: Assumption (A1) We assume both X and C to be nonnegative random variables. We also assume E(X) < +∞ and E(C) < +∞.

In this section, we want to estimate the hazard rate h X of X. This hazard rate may be expressed as the following nonstandard quotient, where S X is the survival function:

h X = f X S X = f X S C S X S C = f X S C S X∧C .
The idea is to estimate separately the numerator f X S C and the denominator S X∧C .

2.2. Construction of the estimator for the numerator. It is rather easy to get an estimator of the numerator f X S C , and more precisely of its projection on the space

(2)

S m := {t ∈ L 2 (R), supp(t * ) ⊂ [-πm, πm]}.
For a square-integrable function g, let us denote g m its orthogonal projection on S m , such that g * m (x) = g * (x)1 |x|≤πm . Then (f X S C ) m is estimated by the following deconvolution estimator:

(

3) (f X S C ) m (x) = 1 2πn n j=1 πm -πm e -iux δ j e iuY j f * ε (u) du. Indeed E(δ 1 e iuY 1 ) = E(1 X 1 ≤C 1 e iu(X 1 ∧C 1 ) e iuε 1 ) = E(1 X 1 ≤C 1 e iuX 1 )f * ε (u) = E(S C (X 1 )e iuX 1 )f * ε (u) = (f X S C ) * (u)f * ε (u). Therefore E( (f X S C ) m (x)) = 1 2π πm -πm e -iux (f X S C ) * (u)du := (f X S C ) m (x).
Under integrability conditions, (f X S C ) m (x) tends to f X S C (x) when m tends to infinity by the Fourier inverse formula. Then the risk bound of (f X S C ) m can easily be deduced from [START_REF] Comte | Penalized contrast estimator for adaptive density deconvolution[END_REF]:

E( (f X S C ) m -(f X S C ) 2 ) ≤ f X S C -(f X S C ) m 2 + E(δ 1 ) 2π πm -πm du |f * ε (u)| 2 where the bias term f X S C -(f X S C ) m 2 = 1 2π |u|≥πm |(f X S C ) * (u)| 2 du
is decreasing with m while the variance term, (2π

) -1 E(δ 1 ) πm -πm du/|f * ε (u)| 2
, obviously increases. As a consequence, we wish to choose m so that a good compromise is reached. This is done by estimating each term. The variance is estimated by a quantity pen 1 (m) proportional to the bound on the variance. Since

f X S C -(f X S C ) m 2 = f X S C 2 -(f X S C ) m 2
, the bias term can be replaced by -(f X S C ) m 2 in the search for an optimal m, as it amounts to omit a term which does not depend on m. Therefore, this term is estimated by -(f X S C ) m 2 . This explains why we select m1 as follows.

(4) m1 = arg min

m∈{1,...,m n,1 } (-(f X S C ) m 2 + pen 1 (m)),
where m n,1 is such that m n,1 ≤ n and

(5) pen

1 (m) = κ 1 n 1 n n k=1 δ k log(J 1 (m))J 1 (m), with J 1 (m) = 1 2π πm -πm du |f * ε (u)| 2 .
In pen 1 (m), the constant κ 1 is calibrated from preliminary simulations.

Following [START_REF] Comte | Penalized contrast estimator for adaptive density deconvolution[END_REF], applying Talagrand's Inequality,

E( (f X S C ) m1 -f X S C 2 ) ≤ C inf m∈{1,...,m n,1 } (f X S C ) m -f X S C 2 + pen 1 (m) + C n
for C and C two constants which do not depend on n. In other words, the estimator (f X S C ) m1 realizes the adequate squared-bias/variance trade off, up to the multiplicative factor C and the negligible residual term C /n.

2.3.

Construction of the estimator for the denominator. We now wish to estimate the denominator S X∧C = S X S C . Note that under assumption (A1), the survival functions are squareintegrable over R + (thus over R if they are extended by 0), as opposed to cumulative distribution functions. This is true, for example, for exponential distributions, classically used in survival analysis: the associated survival functions are clearly square integrable. We define for x ≥ 0, the following estimator of S X∧C , as proposed by [START_REF] Dattner | On deconvolution of distribution functions[END_REF]:

(6) (S X∧C ) m (x) = 1 2 + 1 πn n j=1 Re πm 0 1 iu e iu(Y j -x) f * ε (u) du.
While only the pointwise risk of this estimator is studied in [START_REF] Dattner | On deconvolution of distribution functions[END_REF], we want hereby to compute the integrated L 2 -risk of (S X∧C ) m . It is not trivial from (6) why this integrated risk is properly defined. Therefore, before proceeding to the study of the integrated risk we consider an alternate expression of (S X∧C ) m using (1/π) +∞ 0 sin(v)/vdv = 1/2, as follows: 7) is not mandatory because the integral is real, given the symmetry of the domain and the properties of the function under integration.

(S X∧C ) m (x) = Re   1 2πn n j=1 πm -πm e -iux iu e iuY j f * ε (u) -1 du   + ψ m (x) (7) with ψ m (x) = - 1 2iπ |u|≥πm e -iux u du = 1 π +∞ πm sin(ux) u du. Note that ψ * m (u) = -1/(iu)1 |u|≥πm . This implies by Parseval formula that +∞ 0 |ψ m (x)| 2 dx = (1/2π 2 )m -1 . Moreover, the real part Re(•) in formula (
Then, in order to compute the integrated L 2 -risk of (S X∧C ) m , we can see ( 7) as a deconvolution estimator of S * X∧C . First, notice that S * X∧C (u) =

+∞ 0 e iuv S X∧C (v)dv is well defined under assumption (A1) because S X∧C is integrable and square integrable on R + , its support. Then, let us introduce the following estimate of S * X∧C (u): for all u,

(8) Ŝ * X∧C (u) = 1 n 1 iu n j=1 e iuY j f * ε (u) -1 .
Lemma 1. The estimator Ŝ * X∧C given by ( 8) is well defined on R and is an unbiased estimate of S * X∧C (u).

The estimator (S X∧C ) m written as ( 7) can be seen as the Fourier inversion of (8). Here, however, the Fourier inversion is done with a cutoff πm on the first part of the estimator, which is not integrable, and on the whole real line for the non random part which has a known value. This allows us to write

(S X∧C ) m (x) = Re 1 2π πm -πm e -iux Ŝ * X∧C (u)du + ψ m (x) = 1 2π πm -πm e -iux Ŝ * X∧C (u)du + ψ m (x).
We emphasize that, in practice, the ψ m (x) term is a very useful correction of the estimator for x ∈ [0, 1]. We are now able to study the integrated L 2 -risk and prove the following result.

Proposition 1. Let (S X∧C ) m be defined by ( 6). Under assumption (A1), we have

E( (S X∧C ) m -S X∧C 2 ) ≤ 1 π |u|≥πm |S * X∧C (u)| 2 du + 2 π 2 m + 4 πn πm 1 du u 2 |f * ε (u)| 2 + c n where c is a positive constant, c = (E(Y 2 1 )/π) 1 0 du/|f * ε (u)| 2 .
The first two terms in the upper bound are squared bias terms decreasing when m increases, the third is a variance term which increases with m; the last term is a negligible residual. Contrary to the numerator estimator, the bias decreases at a slower rate. This is due to the term 1/(π 2 m) and to

S * X∧C (u) = f * X∧C (u) -1 iu = f * X∧C (u) iu - 1 iu which implies 1 2π |u|≥πm |S * X∧C (u)| 2 du = O 1 m .
This slow bias order is due to the discontinuity in 0 of survival functions for positive random variables, while computing a global risk over R + . Even with a slow bias decrease rate, we could still obtain a satisfactory convergence rate for the estimator. Indeed, the noises we have in mind must also have lower bounded supports. For example, an exponential distribution for ε yields a variance term of order m/n. The optimal value m opt for the cutoff given by the bias-variance compromise is such that m opt = O( √ n) and the resulting rate is O(n -1/2 ), which is good for a nonparametric deconvolution problem.

All these considerations being asymptotic, we propose a finite sample model selection strategy for choosing m. Let us denote by

(9) J 2 (m) := 1 π πm 1 du u 2 |f * ε (u)| 2 , and pen 2 (m) = κ 2 log(n) J 2 (m) n ,
where κ 2 is a constant to be calibrated by simulations. Note that the lower bound of the integral is 1 so that the integral is properly defined. Then, setting

(10) m2 = arg min m∈{1,...,m n,2 } (-(S X∧C ) m 2 + 3 2π 2 m + pen 2 (m)),
for m n,2 such that m n,2 ≤ n and J 2 (m n,2 ) ≤ n, we obtain an adaptive estimator of S X∧C , which is rather simple to implement, compared to the pointwise procedure of [START_REF] Dattner | On deconvolution of distribution functions[END_REF]. We can prove Theorem 1. Let (S X∧C ) m be defined by ( 6) and m2 by (10). Then there exists a numerical constant κ 0 , such that for κ 2 ≥ κ 0 , we have

E( (S X∧C ) m2 -S X∧C 2 ) ≤ inf m∈{1,...,m n,2 } 3 π |u|≥πm |S * X∧C (u)| 2 du + 2 π 2 m + 4pen 2 (m) + C n
where C is a constant depending on f * ε . From the proof we find that κ 0 = 48 suits, but this theoretical value is too large in practice (see Section 3).

Our adaptive procedure (S X∧C ) m2 has the simplicity of choosing a unique global cutoff m2 for m. This is an advantage compared to the pointwise selection procedure described in [START_REF] Dattner | On deconvolution of distribution functions[END_REF], where the selection is repeated for each point. As a counterpart, we obtain a theoretical global rate which is not as good as the pointwise one. This is due to the fact the pointwise strategy avoids the point x = 0 where a discontinuity occurs.

2.4. Construction of the estimator of h X . The two estimators (f X S C ) m1 and (S X∧C ) m2 allow us to build the final estimator of the hazard rate h X as a quotient estimator. To prevent the denominator to get small, a truncation is required when computing the quotient. The estimator of h X (x) is finally

(11) ĥ m1 , m2 (x) = (f X S C ) m1 (x) (S X∧C ) m2 (x) 1 (S X∧C ) m2 (x)≥λ/ √ n
where λ is a constant to be calibrated. Note that, heuristically, the resulting risk of ĥ m1 , m2 is the addition of the risks of the numerator and the denominator, up to a multiplicative constant.

Simulation

3.1. Design of simulation. Simulations are used to evaluate the performances of the estimator. For each design of simulations, 500 datasets are simulated. We consider samples of size n = 400, 1000. Data are simulated with a Laplace noise with variance σ 2 = 2b 2 as follows:

f ε (x) = 1 2b e -|x|/b and f * ε (x) = 1 1 + b 2 x 2 with b = 1/(2 √ 5) or b = 1/( √ 5
). We consider four densities for X.

(1) Mixed Gamma distribution: X = 1/ √ 5.48W with W ∼ 0.4Γ(5, 1) + 0.6Γ(13, 1) (2) Beta distribution: X ∼ B(2, 5)/ √ 0.025 (3) Gaussian distribution: X ∼ |N (5, 1)| (4) Gamma distribution: X ∼ Γ(5, 1)/ √ 5

These densities are normalized with unit variance, thus allowing the ratio 1/σ 2 to represent the signal-to-noise ratio, denoted s2n. We considered signal to noise ratios of s2n = 2.5 and s2n = 10 in our simulations.

The censoring variable C is simulated with an exponential distribution, with parameter chosen to ensure 20% or 40% of censored observations. 3.2. Estimator implementation. We first describe the implementation of the numerator (f X S C ) m1 . The penalty depends on J 1 (m), which is computed by discretization of the integral. Then we compute pen 1 (m) defined by ( 5) with the choice κ 1 = 2, obtained after a set of simulation experiments to calibrate it. We consider m n,1 = argmax(m ∈ N, J 1 (m)/n ≤ 1). Following, we have the final estimation of m1 defined by (4). By plugging (4) into (3) we obtain (f X S C ) m1 which is our estimator for the numerator.

For the implementation of the denominator (S X∧C ) m2 , the penalty depends on J 2 (m), which is computed by discretization of the integral. We take pen 2 (m) as defined by (9) with κ 2 = 5, after a set of simulation experiments to calibrate it. We define m n,2 = argmax(m ∈ N, J 2 (m)/n ≤ 1). Following, we have the final estimation of m2 defined by (10). By plugging this in (6) we obtain our estimator for the denominator (S X∧C ) m2 .

Finally, we estimate h X as a quotient:

ĥ m1 , m2 (x) = (f X S C ) m1 (x) (S X∧C ) m2 (x) 1 (S X∧C ) m2 (x)≥λ/ √ n ,
with the numerical constant λ = 0.1. 1. Note that even based on 500 simulated datasets, a variability remains in the MISE results, that might be due to the quotient estimator. Therefore, only the trends of the MISE should be interpreted.

Results are compared to estimators obtained in the three following cases: 1/ data with no noise and no censoring, 2/ data with no noise but censoring, 3/ data with noise but no censoring. These three cases can be considered as benchmarks for our situation including both noise and censoring. For case 1, h X is estimated as a quotient of a projection estimator of f X with a trigonometric basis and a Kaplan-Meier estimator of S X . For case 2, h X is estimated as a quotient with numerator and denominator adapted from (f X S C ) m1 and (S X∧C ) m (removing the noise 1/|f * ε |). For case 3, h X is estimated as the quotient of the projection estimator of f X with the trigonometric basis and an estimator of S X directly deduced from (S X∧C ) m . Note that trigonometric polynomials are easy to implement but are sometimes subject to bad side-effects.

Table 1 shows that the MISE obtained with the new estimator are close to the MISE obtained with the more standard estimators without noise or without censoring. The results are satisfactory for the four distributions of X. The MISE are reduced when n increases, whatever the censoring level and the signal to noise ratio. Similarly, the MISE decreases when the censoring level decreases, whatever the value of n and the signal to noise ratio.

We also compare the MISE obtained for ĥ m1 , m2 with the MISE obtained on the same noisy and censored data but modeling either only the noise, or only the censoring, or neither the noise nor the censoring. Results are presented in Table 2 for data with 20% of censoring, small noise (s2n = 10) and n = 400, 1000. We see that when the model is misspecified, the MISE increases. This is especially true when censoring is neglected (two last columns). Neglecting the noise increases the MISE in the Gaussian and the Gamma case. For the Mixed Gamma and the Beta distributions, the MISE are of the same order in the first two columns, when censoring is appropriately modeled.

Table 2. MISE×100 of the estimation of h X , compared with the MISE on the same noisy and censored data but assuming in the modeling either only the noise, or only the censoring, or neither the noise nor the censoring. MISE was averaged over 500 samples. Data are simulated with a Laplace noise, and an exponential censoring variable with 20% of censoring, small noise (s2n = 10) and n = 400 or n = 1000.

Application to length of pregnancy

This work was motivated by the problem of estimating the physiological length of pregnancy, i.e. the time between spontaneous oocyte fertilization and spontaneous delivery.

Although many estimates have been reported, usually of around 40 weeks following last menstrual period (i.e. around 38 weeks after conception), they all rely on imperfect dating of the time origin since the precise time of conception remains unknown in spontaneously conceived pregnancies. In practice, the onset of pregnancy may be estimated by adding two weeks to the last menstrual period, by biochemical tests and also by fetal ultrasound, which is in many cases the preferred method [START_REF] Stirnemann | Day-specific probabilities of conception in fertile cycles resulting in spontaneous pregnancies[END_REF]. Prediction of the time origin using ultrasonographic measurement of fetal crown-rump is not exact and one should take into account this prediction error.

In such data, censoring may occur because of medically planned deliveries whenever maternal or fetal conditions require delivery prior to spontaneous labor. For example, induction (or cesarean section) may be offered whenever women reach 41 weeks and is medically indicated at 41 weeks and 3 days in the Obstetrics unit of Paris Necker Hospital, to avoid fetal complications of prolonged pregnancy. Other medical indications for timely delivery include maternal conditions such as preeclampsia or fetal conditions such as growth restriction or fetal malformations. For obvious reasons, termination of pregnancy or intrauterine fetal demise also require preterm induction of labor.

As already explained, the prediction error using ultrasonographic measurement of fetal crownrump affects the time origin. Therefore it impacts both censoring times and the variable of interest which is the occurrence of a spontaneous onset of labor. This situation refers to the model studied in Section 2. In the following, we consider the measurement error as a Gaussian distribution with mean=0 and standard deviation of 0.35 weeks, as estimated by [START_REF] Stirnemann | Day-specific probabilities of conception in fertile cycles resulting in spontaneous pregnancies[END_REF].

The data we consider here is a sample of 8960 consecutive singleton pregnancies followed in the department of obstetrics, Necker teaching hospital in Paris, from the routine first trimester ultrasound around 12 weeks to delivery between 2011 and 2014. Dating of conception was performed by ultrasonographic measurement of crown-rump length in all cases. In this dataset, censoring occurred in 3410/8960 (38%) cases. The outcomes in the sample are presented in Table 3.

We estimated the hazard rate of the length of pregnancy for spontaneous delivery using estimator (11). The resulting hazard rate for spontaneous delivery is presented in Figure 1. This function increases rapidly from 37 weeks onwards reaching its maximum at 40 weeks and 6.5 days followed by a rapid decrease. In this population this result is markedly different from the usual estimate of 40 weeks that is considered in clinical practice. Therefore, our results would suggest that the true underlying length of pregnancy is longer than observed using noisy data. For comparison, in Figure 1, we also added the hazard rate function for misspecified models that do not model censoring and measurement error. In agreement with the simulation study (Table 2), the misspecified models show significant deviations from our estimator. Interestingly, the model that neglects dating error yields a maximum hazard at 40 weeks and 5 days, 1.5 days earlier than what we found using our estimator. Such a difference may impact on several aspects of clinical care and pregnancy management, when Hazard rate for spontaneous delivery estimated from the noisy and censored dataset of 8960 pregnancies. Solid line: our estimator (11) modeling both measurement error and censoring; broken line: censoring is modeled but not measurement error on time origin; dotted line: neither censoring nor measurement error is modeled seeking to optimize obstetric protocols at term by reducing unnecessary inductions and planning follow-up in prolonged pregnancies.

Discussion and perspectives

In this section, we discuss an alternative model that could be considered when both noise and censoring contaminate the data. Assume that we observe

(W j = (X j + ε j ) ∧ C j , ∆ j = 1 Z j ≤C j )
for j = 1, . . . , n and where Z j = X j + ε j . Despite its simple presentation, seemingly closely related to model (1), this model is more difficult to work with. In this case, the most natural quantity to estimate is the density f X of X, rather than the hazard function. However, we were not able to prove the upper bound of the L 2 risk of the estimator of f X .

More precisely, we consider the estimator by deconvolution of f X of X, as follows:

fX,m (x) = 1 2π πm -πm e -iux f * Z (u) f * ε (u) du where (12) f * Z (u) = 1 n n j=1 ∆ j S C (W j ) e iuW j .
The censoring correction ∆ j /S C (W j ) is standard for such data and sometimes called "Inverse Probability Censoring Weight" (IPCW) in the literature. Then setting f m such that f * m = f * 1 [-πm,πm] , we would easily get under the assumption

Assumption (A2): τ Z < τ C where τ L = sup{y, S L (y) > 0}, E( fX,m -f 2 ) ≤ f -f X 2 + 2E 1 S C (Z 1 ) J(m) n .
Taking the variance term as a penalty for model selection would give an adaptive estimator with good theoretical properties under (A2).

As S C is unknown, it can be estimated with the Kaplan-Meier estimator ŜC , with the modification suggested by [START_REF] Lo | Density and hazard rate estimation for censored data via strong representation of the Kaplan-Meier estimators[END_REF]:

ŜC (y) = W (i) ≤y n -i + 1 n -i + 2 1-∆ (i)
where (W (i) , ∆ (i) ) is ordered following the W j 's. Then this estimator can be plugged into (12) to obtain the estimator

(13) fX,m (x) = 1 2π πm -πm e -iux f * Z (u) f * ε (u) du, with f * Z (u) = 1 n n j=1 ∆ j ŜC (W j ) e iuW j .
However, Assumption (A2) is no longer adequate to handle the substitution of S C by its estimator, and contradictory conditions appear. To obtain an upper bound of the L 2 risk of fX,m , one would need to restrict the variables to belong a fixed compact set, which is rather standard in the the context of estimation with censored data (see the discussion in [START_REF] Gross | Nonparametric estimation and regression analysis with left-truncated and right-censored data[END_REF]), but this is not possible due to the additional deconvolution step: it leads to an assumption regarding the Z's that is contradictory. Thus, the estimation of f X remains an opened theoretical question in this model. Note that we implemented the method and provide in Table 4 the empirical results on simulation: from practical point of view, the procedure seems to work and gives coherent results.

Proofs

We recall the following version of Talagrand inequality.

Lemma 2. Let T 1 , . . . , T n be independent random variables and ν n (r) = (1/n) n j=1 [r(T j ) -E(r(T j ))], for r belonging to a countable class R of measurable functions. Then, for > 0, ( 14) Var(r(T j )) ≤ v.

E[sup r∈R |ν n (r)| 2 -(1 + 2 )H 2 ] + ≤ C v n e -K 1 nH 2 v + M 2 n 2 C 2 ( ) e -K 2 C( ) √ nH M with K 1 = 1/6, K 2 = 1/(21 √
Inequality ( 14) is a straightforward consequence of the Talagrand inequality given in [START_REF] Klein | Concentration around the mean for maxima of empirical processes[END_REF]. Moreover, standard density arguments allow us to apply it to the unit ball of spaces.

The following elementary inequalities will be also used:

(15) ∀u ∈ R, ∀a ∈ R, sin(u) u ≤ 1 and e iua -1 u ≤ |a|. 6.1. Proof of Lemma 1. Let us first remark that Ŝ * X∧C is well defined on R because lim u→0 e iuY j -f * ε (u) iu = Y j -E(ε 1 ). Moreover lim u→0 Ŝ * X∧C (u) = 1 n n i=1 Y j -E(ε 1 ) which tends a.s. when n grows to infinity to E(Y 1 -ε 1 ) = E(X 1 ∧ C 1 ) = S * X∧C (0).
Then we prove that Ŝ * X∧C is an unbiased estimate of S * X∧C . We have

E[ Ŝ * X∧C (u)] = 1 iu E[e iu(X∧C) -1] = 1 iu (e iuz -1)f X∧C (z)dz.
Then, noticing that (e iuz -1)/(iu) = z 0 e iuv dv and that

+∞ 0 +∞ 0 |e iuv f X∧C (z)1 v≤z |dvdz ≤ E(X ∧ C) < ∞, the Fubini Theorem implies that E[ Ŝ * X∧C (u)] = +∞ 0 z 0 e iuv dv f X∧C (z)dz = +∞ 0 e iuv +∞ v f X∧C (z)dz dv = +∞ 0 e iuv S X∧C (v)dv = S * X∧C (u). 6.2. Proof of Proposition 1. Let us set (S X∧C ) m = (S X∧C ) m -ψ m (x). Clearly, E( S X∧C -(S X∧C ) m 2 ) = S X∧C -(S X∧C ) m + ψ m 2 + E( (S X∧C ) m -(S X∧C ) m 2 ) ≤ 2 S X∧C -(S X∧C ) m 2 + 2 ψ m 2 +E( (S X∧C ) m -(S X∧C ) m 2 ), (16) where
(S X∧C ) m is such that (S X∧C ) * m = S * X∧C 1 [-πm,πm] . First we have, by Parseval formula, (17) 2 S X∧C -(S X∧C ) m 2 = 1 π |u|≥πm |S * X∧C (u)| 2 du.
Next, as ψ m is the Fourier transform of 1 |x|≥πm /(2iπx), we have ψ * m (u) = -1/(iu)1 |u|≥πm , and (18)

ψ m 2 = (1/2π) ψ * m 2 = 1/(π 2 m).
For the last term, we use Parseval equality again: 17), ( 18) and ( 19) into ( 16) gives the result of Proposition 1. and using that ∀x, y ≥ 0, 2xy ≤ 4x 2 + y 2 /4, we get Var(r t (Y j )) ≤ v.

E( (S X∧C ) m -(S X∧C ) m 2 ) = 1 2π πm -πm E(| Ŝ * X∧C (u) -S * X∧C (u)| 2 )du. Let us set Ŝ * X∧C (u) -S * X∧C (u) = 1 n 1 iu n j=1 Z j (u) -E(Z j (u)) f * ε (u) with Z j (u) = e iuY j -1. We notice that e iuY j -f * ε (u) -E(e iuY j -f * ε (u)) = Z j (u) -E(Z j (u)). Thus, E(| Ŝ * X∧C (u) -S * X∧C (u)| 2 ) = 1 nu 2 Var(Z 1 (u)) |f * ε (u)| 2 ≤ 1 nu 2 E(|e iuY 1 -1| 2 ) |f * ε (u)| 2 = 4 n E(sin 2 (uY 1 /2)) u 2 |f * ε (u)| 2 . Thanks to inequality (15), we bound this term by (1/n)E(Y 2 1 ) for |u| ∈ [0, 1] and, using | sin(z)| ≤ 1, ∀z ∈ R, by 4/(nu 2 ) for |u| > 1. We get (19) E( (S X∧C ) m -(S X∧C ) m 2 ) ≤ E(Y 2 1 ) πn 1 0 du |f * ε (u)| 2 + 4 πn πm 1 du u 2 |f * ε (u)| 2 . Plugging (
T ≤ 2 ( S X∧C ) m2 -(S X∧C ) m sup t∈S m∨ m2 , t =1 |ν n (t)| ≤ 1 4 (S X∧C ) m2 -(S X∧C ) m 2 + 4 sup t∈S m∨ m2 , t =1 |ν n (t)| 2 ≤ 1 2 (S X∧C ) m2 -S X∧C 2 + 1 2 S X∧C -(S X∧C ) m 2 + 4 sup t∈S m∨ m2 , t =1 |ν n (t)| 2 . ( 23
Similarly to previous computation, we get H 2 = J 2 (m ∨ m )/n, v = J 2 (m ∨ m ) and M = J 2 (m ∨ m ). Moreover we take = 6 log(n 2 ) ∨ 1 = 6 log(n 2 ) for n ≥ 2, and we get We notice that card(M n ) ≤ n and thus card(M n )e -K 2 √ n is bounded to get the result.

  Figure1. Hazard rate for spontaneous delivery estimated from the noisy and censored dataset of 8960 pregnancies. Solid line: our estimator (11) modeling both measurement error and censoring; broken line: censoring is modeled but not measurement error on time origin; dotted line: neither censoring nor measurement error is modeled

  1 and C a universal constant and where sup r∈R r ∞ ≤ M, E sup r∈R |ν n (r)| ≤ H, sup r∈R

  m ∨ m )e -log(n 2 ) + J 2 (m ∨ m ) n e -K 2 √ n using that ≥ 1. Now we have J 2 (m ∨ m ) ≤ n, by definition of M n,2 so that m ∈Mn E sup t∈S m∨m , t =1 |ν n,1 (t)| 2 -3p(m, m )

Table 1 .

 1 MISE×100 of the estimation of h X , compared with the MISE obtained when data are not censored, or not noisy, or neither censored nor noisy. MISE was averaged over 500 samples. Data are simulated with a Laplace noise, and an exponential censoring variable. Results. The values of the MISE are computed from 500 simulated data sets, for each density and simulation scenario and are given (multiplied by 100) in Table

	s2n = 10		0% censoring	20% censoring	40% censoring
			n = 400 n = 1000	n = 400 n = 1000	n = 400 n = 1000
	f X Mixed Gamma with noise	0.723	0.285	0.857	0.423	1.426	0.726
		without noise	0.793	0.287	1.120	0.381	2.071	0.690
	f X Beta	with noise	1.405	0.813	1.780	1.019	2.328	1.354
		without noise	1.288	0.641	1.767	0.793	2.454	1.050
	f X Gaussian	with noise	0.598	0.250	1.238	0.703	6.580	6.077
		without noise	0.656	0.201	2.009	0.628	8.934	5.808
	f X Gamma	with noise	0.805	0.361	0.843	0.364	0.988	0.408
		without noise	0.684	0.275	0.865	0.268	1.058	0.327
	s2n = 2.5		0% censoring	20% censoring	40% censoring
			n = 400 n = 1000	n = 400 n = 1000	n = 400 n = 1000
	f X Mixed Gamma with noise	1.129	0.497	1.235	0.645	1.727	0.915
		without noise	0.793	0.287	1.120	0.381	2.071	0.690
	f X Beta	with noise	2.075	1.154	2.878	1.455	3.437	1.847
		without noise	1.288	0.641	1.767	0.793	2.454	1.050
	f X Gaussian	with noise	0.950	0.437	1.978	1.031	5.526	4.940
		without noise	0.656	0.201	2.009	0.628	8.934	5.808
	f X Gamma	with noise	1.173	0.613	1.379	0.660	1.538	0.821
		without noise	0.684	0.275	0.865	0.268	1.058	0.327
	3.3.

Table 3 .

 3 Outcomes and main reasons for censoring in a sample of 8960 singleton pregnancies followed in Necker Enfants Malade teaching hospital

	Termination of pregnancy or intrauterine fetal demise 110
	Induction at term (≥ 41 weeks)	903
	Induction for maternal or fetal reasons before term	1490
	Planned cesarean section	907
	Spontaneous labor	5550

Table 4 .

 4 MISE×100 of the estimation of f X , compared with the MISE obtained when data are not censored, or not noisy, or neither censored nor noisy. MISE was averaged over 100 samples. Data are simulated with a Laplace noise, and an exponential censoring variable.

	s2n = 10		0% censoring	20% censoring	40% censoring
			n = 400 n = 1000	n = 400 n = 1000	n = 400 n = 1000
	f X Mixed Gamma with noise	0.203	0.097	0.381	0.171	0.445	0.185
		without noise	0.181	0.082	0.258	0.087	0.259	0.102
	f X Beta	with noise	0.271	0.193	0.353	0.258	0.432	0.255
		without noise	0.975	0.579	0.280	0.163	0.349	0.191
	f X Gaussian	with noise	0.139	0.054	0.527	0.255	1.719	0.973
		without noise	0.481	0.237	0.146	0.070	0.452	0.127
	f X Gamma	with noise	0.290	0.138	0.316	0.166	0.371	0.170
		without noise	0.549	0.235	0.196	0.083	0.211	0.114

  Now we split ν n (t) = ν n,1 (t) + R n (t) with where c is defined in Proposition 1. For the other term we use the following Proposition.Proposition 2. Let p(m, m ) = n -1 log(n 2 )J 2 (m ∨ m ), thenThe proof of Proposition 2 follows from Talagrand inequality and is proved below. Now we notice that 3κ 2 p(m, m ) ≤ 6pen 2 (m) + 6pen 2 (m ) so that 8E supt∈S m∨ m2 , t =1 |ν n,1 (t)| 2 -pen 2 ( m2 )/8 ≤ 8E sup t∈S m∨ m2 , t =1 |ν n,1 (t)| 2 -3p(m, m2 ) -1 ≤ 0 i.e. κ 2 ≥ 48. Plugging (25) and (26) in (24), we obtain, ∀m ∈ {1, . . . , m n,2 },E( (S X∧C ) m2 -S X∧C 2 + ψ m2 2 ) ≤ 3 S X∧C -(S X∧C ) m 2 + ψ m X∧C ) m2 -S X∧C 2 = (S X∧C ) m2 -S X∧C + ψ m2 2 ≤ 2( (S X∧C ) m2 -S X∧C Proof of Proposition 2.We set, for t ∈ S m , and we apply Inequality of Lemma 2 to R = S m∨m , by using standard arguments of continuity of t → ν n,1 (t) and density of a countable subset of S m∨m . Next we have to compute H 2 , M, v such that sup

								+
	) (26)					≤	+( c n + 48 κ 2 48 -1)E(pen 2 ( m2 )) + κ 2 pen 2 (m),	48 κ 2	pen 2 (m)
	Plugging (23) into (22) yields 1 2 (S X∧C ) m2 -S X∧C 2 + for 48/κ 2 2 1 2 ψ m2 2 ≤ 3 2 S X∧C -(S X∧C ) m 2 + 1 2 ψ m 2 + pen 2 (m) +4 sup |ν n (t)| 2 -pen 2 ( m2 ). (24) +2(1 + 24/κ 2 )pen 2 (m) + 8c . n t∈S m∨ m2 , t =1 To conclude, we notice that
	(S 2 + ψ m2	2 )
	R n (t) = which implies t We have 1 2π |u|≤1 E( (S X∧C ) m2 -S X∧C 2 ) ≤	inf m∈{1,...,m n,2 }	6 S X∧C -(S X∧C ) m	2 + 2 ψ m	2 + 6pen 2 (m) +	c n	,
	sup c = 16c, which is the announced result. |ν n (t)| 2 ≤ 2	sup		R 2 n (t) + 2	sup	|ν n,1 (t)| 2 .
	t∈S m∨ m2 , t =1 Moreover by Schwarz Inequality, we have t∈S m∨ m2 , t =1 E sup t∈S m∨ m2 , t =1 R 2 n (t) ≤ 6.4. r t (x) = 1 2π t * (u) 1 2π E , so that |u|≤1 e iux 1 |u|≥1 t∈S m∨ m2 , t =1 iuf * du ε (u) | Ŝ * X∧C (u) -S * X∧C (u)| 2 du , and from the proof of Proposition 1, we easily get (25) E sup t∈S m∨ m2 , t =1 R 2 n (t) ≤ 2 2π E(Y 2 1 ) n 1 0 du |f * ν n,1 (t) = n 1 n j=1 [r Classically we write c ε (u)| 2 = n E sup |ν n,1 (t)| 2 -3p(m, m2 ) ≤ E sup |ν n,1 (t)| 2 -3p(m, m )
	t∈S m∨ m2 , t =1			+	m ∈Mn	t∈S m∨m , t =1
	t∈S m∨m	sup	E	sup				t∈S m∨m	1 n	n j=1

* (-u)( Ŝ * X∧C (u) -S * X∧C (u))du, ν n,1 (t) = ν n (t) -R n (t). t∈S m∨ m2 , t =1 |ν n,1 (t)| 2 -3p(m, m2 ) + ≤ c n . t (Y j ) -E(r t (Y j ))]. + x∈R |r t (x)| ≤ M, E sup t∈S m∨m |ν n (r t )| ≤ H, sup

Acknowledgments

The authors thank Pr Jean-Christophe Thalabard for his advices, suggestions and support.

6.3. Proof of Theorem 1. Let S m = {t ∈ L 2 (R), Supp(t * ) ⊂ [-πm, πm]}. Note that the estimator (S X∧C ) m = (S X∧C ) m -ψ m (x) satisfies

with Ŝ * X∧C given by ( 8). One can see this by noticing that

On the other hand we have (S X∧C ) m 2 = (S X∧C ) m 2 + ψ m 2 , since the support of the Fourier transforms of the functions in the norms are disjoint. Therefore

Therefore, definition (10) of m2 can be written

The equality (20) for m2 implies that, ∀m ∈ {1, . . . , m n,2 } and for all t ∈ S m ,

Taking t = (S X∧C ) m and using (21), this can be rewritten

Then considering the term in the second line of ( 22


Jingtao Zhu

Yuchun Tu

Yanyan Yuan

Zhixiang Feng

Haochuan Li

Yizhi Zhang

Philippe Jonnard

Christian Mény

Karine Le Guen

Jean-Michel André

Zhanshan Wang
email: wangzs@tongji.edu.cn

Thermal stability of Co/C multilayers

Keywords: Multilayer, X-ray optics, interface, thin film structure, NMR

The structural and interface changes induced by thermal annealing in Co/C multilayers were investigated.

Co/C multilayers with period thickness of 4.1 nm and bi-layer number of 20 were deposited by direct currentmagnetron sputtering.We characterized all samples by using X-ray reflectivity, X-ray diffuse scattering, zero-fieldnuclear magnetic resonancespectroscopy and X-ray diffraction. The results indicate that Co and C atoms mixed during deposition and then after annealing both atoms separated from their mixed region. The annealing process also causes asharpincrease of roughness at interfaces, which can be attributed to the crystallization of Co layers.

Introduction

Multilayer mirrors are widely used as reflective optical elements for applications such as extreme ultraviolet (EUV) lithography [1], solar astrophysics imaging [2], synchrotron radiation [START_REF] Morawe | Proc.SPIE[END_REF] and x-ray spectrometry [START_REF] Hombourger | [END_REF].Near the "carbon window" region (λ~4.5 nm) where the absorption of carbon materials is low, Co/C multilayers have the rather hightheoretical reflectivity [5] and have been produced for Schwarzschild objective [6]. The peak reflectivity of Co/C multilayers was measured to be 14.8% at incident angle of 5°.In x-ray optics of synchrotron radiation beamlines and telescopes for astronomical observations, multilayers are used as monochromators and focusing elements. Co/C multilayers are promising near Co K-edges because of their fine optical properties [START_REF] Bellotti | Proc. SPIE[END_REF].

In order to improve the optical performance, it is necessary to deposit multilayer films that give rise to smooth and sharp interfaces, as the interface imperfections (i.e., roughness anddiffussion) will reduce the reflectance and increase non-specular scattering as well [START_REF] Bellotti | Proc. SPIE[END_REF].Furthermore, multilayer mirrors are usually exposed to high flux of incident photons or endure high heat loads in practical applications such as solar observation or synchrotron radiation, where they might reach temperatures and heat load in excess of 500°C [START_REF] Ziegler | [END_REF] and 100 W/mm 2 [9]. It is known that the structure and interfaces of multilayers may change, e.g. period expansion or contraction, interdiffusion, crystallization, which could reduce the optical performance and accelerate the aging of the multilayer [10]. Thus, it is important to investigate the thermal stability andstructure evolution of the multilayers.

In this paper, we combinedX-ray reflectivity (XRR), X-ray diffuse scattering(XDS), zero-field nuclear magnetic resonance (NMR) spectroscopy and X-ray diffraction (XRD)techniques to investigate the evolution of structure and interfacesof Co/C multilayers upon annealing.

Experiments

The Co/C multilayers were deposited onto polished Si wafer by ultra-high vacuum direct current (DC) magnetron sputtering system with targets of Co (purity 99.95%), C (purity 99.99%). The working gas is argon (99.999%). The base pressure was 1×10 -4 Pa and the sputtering gas pressure was 0.13 Pa. The deposited multilayer was designed with bilayer number N=20, period thickness Λ=4.1nm, and the thickness of Co layer, d Co =1.5nm.The layer thickness was calibrated by using grazing incident X-ray reflectivity measurement. After deposition, the samples were mounted on a plate heated by a wire-wound furnace in a vacuum chamber with a base pressure of 3×10 -4 Pa. The samples were heated from room temperature to 300 and 600°Ckeeping for 1 hour. A detailed description of the deposition and annealing setup can be found in Ref. [11].

The structure of the multilayers was analyzed by using grazing incident XRR on aX-ray diffractometer (D1 system, Bede Inc.) at Cu Kα line (0.154 nm)working in the θ-2θ mode. The angular resolution is 5/1000°. Bragg law corrected for refraction was used to determine the multilayer period. The fitting of the XRR curves performed with Bede Refs software was used to determine thickness, roughness and density of each layer.Rocking curves were measured by using the samediffractometerwiththe detector fixed at the angle of first order Bragg peak.XRD measurements in the symmetric (θ-2θ) geometry were performed on a Rigaku-Dmax-2550V powder diffractometer.

In order to probe the chemical state of the Co atoms within the multilayer, the samples are analyzed byNMR spectroscopy. The NMR spectra represent the distribution of the Co atoms as a function of their resonance frequency, i.e. the hyperfine field experienced by the Co nuclei [12][13]. The NMR resonance frequency is sensitive to the local environment of the probed atoms: nearest neighbor local structure and/or local chemical environment. To enhance the sensitivity, the testing temperature is 2 Kfor all the samples.All spectra have been recorded for different values of radio frequency field strengths allowing for correcting the NMR intensities for a frequency dependent enhancement factor.

Results and discussion

Grazing incident X-ray reflectivity

Thegrazing incident X-ray reflectivity curves of the Co/C multilayers are presented inFigure 1.All the XRR curves were fitted to obtain the roughness, thickness and density of individual layers. The fitted values are listed in Table .1. For the sample annealed at 300°C, the Bragg peaks slightlyshift to lower grazing incident angle with respect to those of as-deposited sample.This means a slightincrease of period thickness. When the sample was annealed at 600°C, we observe anexpansion of 17.5% in period thickness (from 4.09 nm to 4.81nm).

The fitting results reveal that the carbon layer plays a crucial role in the period expansion. After 600°C annealing, the thickness of C layer increases from 2.58nm to 3.37nm.Observation of period expansion in carbon-based multilayersafter annealingwasreported in Cr/C [14],Pt/C [15] and W/C [16] systems. This phenomenon is attributed to the graphitization of the amorphous carbon layers by annealing, which is confirmed by Raman scattering measurements [14,16]. There is a convention of sp3 bond to sp2 bond and amorphous carbon structure to nanocrystalline graphite in carbon layers during annealing [14]. The density of carbon changes linearly with sp 3 fraction [17]. That is to say, the density of carbon should decrease and the thickness of the carbon layers should increase upon annealing. The thickness of Co layers decrease from 1.51nm to 1.44nm. This may be due to the demixing of Co and C atoms in the Co layers after annealing which will be mentioned below. In Figure 1, three well-defined peaks are observed for the as-deposited and 300°C annealed samples up to 3.5°.In the case of the 600°C annealed sample, only first and second Bragg peaks are observed and the reflectance decreased.The intensity of Bragg peaks increases slightly after 300°C, meaning that there is aninterface smoothening effect in Co/C multilayers. The comparison of the reflectivity curves shows that after 600°C annealing, the roughness of interfaces increases dramatically, from 0.4 nm to about 1.0 nm. The roughness deduced from the fitting of XRR (σ)is an overall imperfection of interfaces, including the contributions from both geometrical roughness (σ r )and interdiffusion (σ d ). Possible reasons for thermally induced roughness change in Co/C multilayers are investigated in the following.

Zero-field nuclear magnetic resonance (NMR) spectroscopy

NMR spectra for the Co/C multilayers are presented in Figure 2. All spectra in Fig. 2 (a) are normalized to the surface area of each sample because the intensity is proportional to the size of samples. The spectrum of the Co/C as-deposited multilayer shows a weak signal without any defined feature in the frequency range of pure Co.This indicates that the Co layers are not pure but that alien carbon atoms are mixed with Co over the total Co thickness. The as-deposited multilayers are actually CoC x /C multilayer with probably some concentration gradient at the interface with the C layers. In addition some Co atoms are most probably also situated in non-ferromagnetic phases (and therefore give no NMR signal) with a large content of C.The metal-containing layer is a Co-C alloy with the density of 6.2-6.5 g/cm 3 that is less than density of massive Co and carbon enriched Co 2 C layers [START_REF] Artyukov | Proc. SPIE[END_REF]. As a result, Co/C multilayer mirrors consist of cobalt-carbon alloy and amorphous carbon layers.

The NMR spectrum of the Co/C 300°C annealed sample shows no defined peak either, meaning that Co and Catoms are still stronglymixed in CoC x layers. There is a broad structureat 110MHz showing that some atomic motion occurred during annealing but the exact origin of this signal is difficult to identify. However,compared to the result of as-deposited samples, Co and C has a tendency of demixing at a temperature of 300°C. Indeed the demixing processbetween the Co and C atoms increases thex-rayopticalcontrastof the Co/Cmultilayers. That is why the X-ray reflectivity is improved at temperature of 300°C.

The spectrum of the Co/C 600°C annealed sample shows a well-defined peak at about 220MHz which corresponds to bulk hcp /fcclikeCo [12,13]. This means that the pure Co regions appear and the diffusion decreasesat the interfaces.We infer from these results that Co and C layers in the Co/C stack are strongly mixed during deposition and then demixingoccursupon annealing. Based on the Miedema's macroscopic atom model [START_REF] Das | [END_REF], the mixing enthalpy of Co-C system is calculated and the results are shown in Fig. 3. It can be seen that the value of mixing enthalpy is positive at any carbon concentration. Thisindicates that the Co-C system is easy to separate into two phases from the view of thermodynamics. This is in agreement with the NMR results. On the other hand, the carbon atoms separating out from CoC x layers may be an additional reason for the expansion of carbon layers.

X-ray diffuse scattering

The NMR results confirm the decrease of intermixing after annealing. Thus we attribute the deterioration of the interfaces to the increase of geometrical roughness. Therefore, X-ray diffuse scattering measurement was performed. Figure 4 shows the rocking curves measuredaround the first Bragg peakfor the three samples. For the as-deposited and 300°C annealed sample there is no significant change of the intensity except a slight shift of peak position.After annealed at 600°C, a high intensity of diffuse scattered x-rays is observed, indicatingsignificant increase of the interfaces roughness.

X-ray diffraction

The x-ray diffraction patterns of Co/C as-deposited and annealed at 300°C and 600°C samples are presented in Fig. 5. A broad intensity can be observed around 44.5° around the hcp(0002)/fcc (111) [START_REF]The diffraction information for Co is obtained from ICDD[END_REF] peak position for the sample annealed at600°C. Compared with the amorphous layers of as-deposited and 300°C annealed sample, this peak indicatesa significant structural change:the amount of pure Co increases upon annealing. This is in good agreement with results of NMR spectra. Co and C atoms may not separate completely because the XRD peak is weak and its width is much broader than reported for Co/Cu multilayers with 1.5 nm thick Co layers [START_REF] Kubota | [END_REF].

That is to say, the Co layers still contain a certain amount of C upon600°C annealing.

The reduction in reflectivity induced by interface roughness becomes more significant as the multilayer period thickness decreases, because the ratio of the interface roughness to the multilayer period thickness directly affects the reflectivity according to the Debye-Waller factor [22]. Crystallization of metal layers in multilayers can cause an increase of geometrical roughness [23]. Crystallites can nucleate at some points during annealing. They have a size limited to the metal layer thicknessin the growth direction and arealso small in the in-plane direction. The spatial distribution of crystallites in the metal layer is replicated as a modulation of the surface height.Since the metal layer become rougher, the nonmetal layer on top of the metal layer also becomes rougher. For the Co/C multilayers annealed up to 600°C, the crystallization of Co layers increases the geometrical roughness and aggravates the overall imperfection of the interfaces.

Fig. 1 .Table. 1 .

 11 Fig. 1. Measured (black circle)and fitted (red solid line) XRR curves of as-deposited and annealed。 Co/C multilayers.

Fig. 2 .

 2 Fig. 2.Zero-field NMR spectra of as-deposited and annealed Co/C multilayers.

Fig. 3 .

 3 Fig.3.Theoretical calculated mixing enthalpy of the Co-C system.

Fig. 4 .

 4 Fig.4. Rocking curves around the first Bragg peak of as-deposited and annealed Co/C multilayers.

Fig. 5 . 4 Conclusion

 54 Fig.5.Diffraction patterns of as-deposited and annealed Co/C multilayers.The curveswere shifted vertically forthe sake of clarity.4 ConclusionSputter-deposited Co/C multilayers annealed up to 600°C (2.6 nm C /1.5 nm Co) 20 areinvestigatedby using XRR, XRD, XDS and NMR techniques. The multilayer is stable up to 300°C and there isa slightincrease in the period and reflectivity.We demonstrate through the NMR spectroscopy that the strong demixingbetween Co and C atoms occursafter annealing at 600°C.This phenomenon is also demonstrated by the Miedema'smodel which shows the phase separation of the Co-C system.The XDS and XRD results show that when the annealing temperature increases up to 600°C, although the demixingof Co and C atoms increases thex-rayopticalcontrastof the Co/Cmultilayers, the formation of small Co crystallinegraincauses the interfaces rough. The demixing and graphitization of the amorphous carbon layersinduceanexpansionof period thickness.

Acknowledgement:

This work was supported by CaiYuanpei Project (PHC N° 30248NF), National Basic Research Program of China (No.2011CB922203), and National Natural Science Foundation of China (Nos. 11375131 and 11305104).


Yunbo Li
email: yunbo.li@emn.fr

Anne-Cécile Orgerie
email: anne-cecile.orgerie@irisa.fr

Jean-Marc Menaud
email: menaud@emn.fr

Opportunistic Scheduling in Clouds Partially Powered by Green Energy

The fast growth of demand for computing and storage resources in data centers has considerably increased their energy consumption. Improving the utilization of data center resources and integrating renewable energy, such as solar and wind, has been proposed to reduce both the brown energy consumption and carbon footprint of the data centers. In this paper, we propose a novel framework oPportunistic schedulIng broKer infrAstructure (PIKA) to save energy in small mono-site data centers. In order to reduce the brown energy consumption, PIKA integrates resource overcommit techniques that help to minimize the number of powered-on Physical Machines (PMs). On the other hand, PIKA dynamically schedules the jobs and adjusts the number of powered-on PMs to match the variable renewable energy supply. Our simulations with a real-world job workload and solar power traces demonstrate that PIKA saves brown energy consumption by up to 44.9% compared to a typical scheduling algorithm.

I. INTRODUCTION

Data centers have been the key system infrastructure for cloud computing. With the emergence of the Future Internet and the dawning of new IT models such as cloud computing, the usage of data centers, and consequently their power consumption, increases dramatically [START_REF] Orgerie | A Survey on Techniques for Improving the Energy Efficiency of Large-scale Distributed Systems[END_REF]. A data center (DC) is a facility used to house tens to thousands of computers and their associated components. These servers are used to host applications available in the Internet, from simple web servers to multi-tier applications, but also some batch jobs [START_REF] Menache | On-demand, Spot, or Both: Dynamic Resource Allocation for Executing Batch Jobs in the Cloud[END_REF]. Besides the ecological impact, the energy consumption is a predominant criteria for DC providers since it determines the daily cost of their infrastructure. As a consequence, power management becomes one of the main challenges for DC infrastructures and more generally for large-scale distributed systems [START_REF] Orgerie | A Survey on Techniques for Improving the Energy Efficiency of Large-scale Distributed Systems[END_REF].

In parallel to the expansion of cloud computing, from this recent years, a new model emerges: decentralized cloud infrastructures [START_REF] Bertier | Beyond the Clouds: How Should Next Generation Utility Computing Infrastructures Be Designed?[END_REF]. To improve the performance of their cloud and to leverage their available infrastructure, telecom operators, like Orange, try to deploy micro data center (20 to 50 servers by micro-DC) at the network border, closer to customers. In this new model, by deploying data centers closer to the user, the response time and throughput would greatly improve.

From an energy point of view, these micro-data centers allow the study of new power supply solutions based on renewable energy, like wind or sun. Using these renewable energy sources can reduce the operating cost but, unfortunately, this kind of energy stays intermittent by nature. To address this problem, two solutions: investing in heavy expensive battery systems to smooth over the day the renewable energy production, or developing new applications management solutions adapted to the electricity production. In this paper, we propose to design a disruptive approach to Clouds resource management which takes advantage of renewable energy availability to perform opportunistic tasks.

The micro-DC receives a fixed amount of power from the regular electrical grid. This power allows it to run the usual tasks. In addition, the micro-DC is also connected to renewable energy sources (such as windmills or solar cells) and when these sources produce electricity, the micro-DC uses it to run more, less urgent, tasks. In order to achieve this energy-aware resource allocation, we distinguish two kinds of jobs to be scheduled on the data center: the web jobs which represent jobs requiring to run continuously (like web server), and the batch jobs which represent jobs that can be delayed and interrupted, but with a deadline constraint. The second type of jobs are the natural candidates of the opportunistic scheduling algorithm. This paper presents PIKA, a framework aiming at reducing the brown energy consumption (ie. from non-renewable energy sources), and improving the usage of renewable energy for mono-site data center. It exploits jobs with slack periods, and executes or suspends them depending on the renewable energy availability. By consolidating the virtual machines (VM) on the physical servers, PIKA adjusts the number of powered-on servers in order for the overall energy consumption to match with the renewable energy supply. Using simulations driven by real-world workloads and solar power traces, we demonstrate that PIKA consumes 44.9% less brown energy and increases by 110.1% the renewable energy integration ratio in comparison with the baseline algorithm from literature.

The remainder of the paper is organized as follows. Section II presents related work. Section III formalizes the problem. We give a brief overview of PIKA framework in Section IV and we presents our formalization of PIKA in Section V. Section VI gives different resource overcommit policies and the experimental setup is explained in Section VII. Section VIII evaluates the various policies based on simulations. Section IX concludes this work.

II. RELATED WORKS

We briefly describe related works in energy consumption management and renewable energy integration in data centers.

A. Energy consumption management in data center(DC)

Energy consumption management has become an important issue for data centers, but solutions mainly focus on the reduction of the brown energy consumption [START_REF] Orgerie | A Survey on Techniques for Improving the Energy Efficiency of Large-scale Distributed Systems[END_REF]. To save energy in a single data center, a common goal is to reduce the number of powered-on PMs (ON PMs). The performance of VM placement algorithm directly affects the number of ON PMs. The problem of VM placement is typically modeled as a n-dimensional bin-packing problem which is NP-Hard [START_REF] Chekuri | On Multi-Dimensional Packing Problems[END_REF]. In this paper, we only consider the PMs CPU and RAM as resource constraints such that the VM placement problem becomes 2-dimension. The sum of VMs resources requirement on a PM can not exceed this PM's capacity. In [START_REF] Wood | Black-box and Gray-box Strategies for Virtual Machine Migration[END_REF]

, Wood et al. propose a metric Volume with vol = 1 1-cpu × 1 1-mem × 1 1-net
where cpu, net and mem are the corresponding utilization levels of that resource for the VM on the PM. The higher utilization of a resource corresponds to a higher volume. However, this solution considers that all the PMs have identical hardware which is different from our context (heterogeneous environment). More important, if we swap two values among the three (e.g. cpu and mem), the volume is still holding the previous value, but may be ambiguous.

The above proposals of VM placement investigate the bin packing problem to use less PMs and to pack more VM. They consider the VM resources requirements as a bin size. But there is a a huge difference between the VM resources requirements and the real resources usage of a VM that may also lead a wastage of resources. In [START_REF] Zhang | Virtual machine migration in an over-committed cloud[END_REF], Zhang et al. design an algorithm (called Scattered) to reduce the number of ON PMs by resource over-commitment. Similarly, in [START_REF] Dabbagh | Efficient datacenter resource utilization through cloud resource overcommitment[END_REF] and [START_REF]Toward energy-efficient cloud computing: Prediction, consolidation, and overcommitment[END_REF], the authors use the resource over-commitment technique to increase the resources usage. Hence, the over-commit resources technique is required for our framework to increase the PM usage and reduce the number of PMs if the VMs are not fully loaded. In addition, a single PM energy consumption is related to its different components [START_REF] Quesnel | Estimating the Power Consumption of an Idle Virtual Machine[END_REF]. For example, the CPU and RAM highly impact the PM performance and the energy consumption, more than the other components. Furthermore, the energy consumption of a PM mainly depends on the number of active cores and threads in comparison with the others hardware such as RAM, network card [START_REF] Fan | Power provisioning for a warehouse-sized computer[END_REF]. Thus, over-committing the CPU resource is considered prior than the other resources in our framework. However, the resource over-commitment may lead to overloading of the PM. Then, we need to migrate the VMs from an overloaded PM to another one. In [START_REF] Beloglazov | Energy efficient resource management in virtualized cloud data centers[END_REF], an optimization power-aware migration has been developed. It allows to migrate the VM to the most energy-efficient PM. However, it does not consider the number of migrations. Moreover, the migration also has an energy overhead that may lead to extra energy consumption and performance degradation.

Since in Cloud environments, each job has different arrival time and lifetime, the static VM placement cannot satisfy the dynamicity of the system. It necessitates adequate optimization operations to adjust the number of ON PMs in the current system. The dynamic VM consolidation techniques have now been widely employed in modern data centers, which enable the data center to decrease the number of ON PMs [START_REF] Orgerie | A Survey on Techniques for Improving the Energy Efficiency of Large-scale Distributed Systems[END_REF].

B. Renewable energy in DC

The energy supply can be roughly divided into fossil fuels energy (brown energy) and renewable energy such as the wind energy and solar energy. Integrating the renewable energy into data center can further reduce the brown energy consumption. Due to the intermittent and variable nature of renewable energy, we aim at scheduling opportunistically the jobs according to the renewable energy availability. In [START_REF] Krioukov | Integrating Renewable Energy Using Data Analytics Systems: Challenges and Opportunities[END_REF] and [START_REF] Chen | Green-aware workload scheduling in geographically distributed data centers[END_REF], the authors point out the slack time is a key feature that enables the jobs to be delayed until the renewable energy becomes available. However, the authors only take into account batch jobs in their work. In real-world workloads, there are massive number of jobs which do not have slack, like web services for instance. Furthermore, they propose a simple scheduling algorithm without VM consolidation and migration while both techniques can improve the energy consumption.

III. PROBLEM FORMULATION

This paper targets small/medium-size mono-site data center (typically between 20 and 150 servers). The data center consists of several Physical Machines (PMs). We assume the computing environment in the data center to be heterogeneous, meaning that the PMs can have different hardware. Each PM has limited resources (CPU, RAM, Network) and has its own disk storage. We assume that the data center has no centralized storage system (such as a NAS for instance) as described in [START_REF] Beldiceanu | The EPOC project: Energy Proportional and Opportunistic Computing system[END_REF]. The PMs with different capacities and performances may potentially lead to different energy consumption for a given VM.

We assume the data center has dual brown and renewable energy power supplies. Each PM has a switch connected with both energy supplies and opts for renewable energy only if there is enough of it. Otherwise, the PMs consume the brown energy from the regular grid. We assume there is no battery or the batteries are only used for emergency cases.

A. Job

Our system does not only accommodate periodic jobs which means that jobs may have different lifetimes and can arrive at anytime. Once a job is submitted by a user, it is encapsulated into an individual Virtual Machine (VM). In our system, a VM is considered as the basic unit of resource allocation. It demands two types of resources from PM: CPU and RAM. The lifetime of a VM depends on the job it accommodates. When a job is finished, the VM s destroyed and its reserved resources on the PM are released.

B. Traces

We studied anonymized traces provided by the University of Nantes, France. Theses traces concern a VM hosting provider with 55 servers. The traces stretch from the 25th of March 2014 to the 6th of July 2014. They consist in the logs for real CPU, RAM, network and disk utilization of each server every 90 seconds. They also contain the client's requests for VMs with CPU and RAM sizes, and the submission dates. These traces present a realistic scenario in our context.

The Figure 1 illustrates the average CPU and RAM utilization of all the PMs during a normal week in the data center. IV. PIKA OVERVIEW Our proposed framework PIKA is designed as a centralized solution. It focuses on minimizing the brown energy consumption in a single small/medium-size data center. As the system is dynamic, PIKA performs the optimization operations periodically. The optimization cycle is defined as a slot, such that the time in our system is divided into a series of continuous slots. As shown in Figure 2, at the beginning of each slot, the broker executes the three main steps. First, the broker checks each PM's state and suspends some jobs from the overloaded PMs. Second, the renewable energy predictor predicts the amount of renewable energy for the current slot and informs the broker about it. Then, the broker determines the number of ON PMs that can be supported by the renewable energy supply. Finally, according to the available resources from these ON PMs, the broker schedules the jobs that can be executed during the current slot. Each component of PIKA is described in the remainder of this section.

Predictor

A. Renewable energy predictor

We make several simplifying assumptions. The renewable energy prediction performs at the beginning of each slot and the predicted renewable energy amount is used for only one slot. The short-term prediction's advantage is that it significantly reduces prediction errors caused by varying weather. Given an accurate prediction on renewable energy, the broker dynamically switches on and off PMs to adjust the energy consumption in order to maximize the renewable energy integration ratio.

B. Gap

PIKA is designed to be aware of the variable and intermittent nature of renewable energy supply, which makes the gap module in PIKA play an important role for the other operations. The renewable energy predictor provides the energy availability that can be consumed in the current slot. The key feature of gap enables the broker to dynamically adjust the number of ON PMs following by the renewable energy availability. We detail how the gap works in Section V-B.

C. Job pool

The key insight of PIKA is to align the energy consumption with the variable renewable energy supplies. Each job is submitted with the following parameters: (t b , T, t e , t d , v cpu i , v ram i ): beginning time, type, execution time, deadline, corresponding the i th VM CPU and RAM requirement. They are pushed into the job pool by the broker. Once the broker gets the energy availability through the renewable energy predictor, it estimates the maximum number of ON PMs via the gap function in the current slot. To increase the chance to exploit more renewable energy, we employ the slack time for jobs into PIKA.

The slack is the most crucial factor in affecting the renewable energy integration ratio. It is given in Equation 1:

j slack i = t d -t b -t e (1)
Through j slack i , the broker enables a job to be delayed and this increases the chance of exploiting the renewable energy. Consequently, we classify the jobs into two types according to their characteristics: web job and batch job.

The web-job is non-interruptible with a little slack (< 1slot). It has a higher priority than other jobs either on resource allocation or scheduling. For instance, when a web job is submitted, the broker pushes it into a specific job pool with higher priority than the waiting queue of the batch jobs. The web job is then placed on a PM which has sufficient resources to meet its VM resources requirements. Unlike the web job, the Batch job can be interrupted or delayed within a slack. Furthermore, because web jobs have higher priority, the batch jobs must wait for all the web jobs to be placed before being allocated only to a ON PM (PM already hosting some web jobs). The batch jobs are not allowed to switch on an OFF PM. Specially, when a batch job's slack is strict inferior to 1, it mutates as a web job that has the same priority as a regular web job. So, the job pool in PIKA is divided in two corresponding pools: the web pool and batch pool. Note that the VM placement process for web job and batch job is individual.

D. Flexible pooling of ON/OFF PMs

The PIKA offers a mechanism for dynamically adjusting the number of ON/OFF PMs that tries to follow the variable renewable energy supply if possible (according to the workload). Firstly, we divide the PMs into two categories: ON PM and OFF PM. Both web job and batch job are capable of being allocated on an already ON PM if it has sufficient resources to meet their VM resources requirements. An OFF PM can be switched on in the following two cases: 1) there is no more resources to meet the web job's VM resource requirements; 2) there is sufficient renewable energy for all the already ON PMs at that time and an extra amount from renewable energy supply to switch on new OFF PMs (called potential ON PM later in the paper). A batch job is allowed to switch on an OFF PM if and only if it is mutated into a web jobs. The general VM placement process as shown Figure 3. In this section, we have given an overview of PIKA's architecture. The next section is dedicated to the formalization of the various algorithms used by PIKA for the resource management and job scheduling.

V. RESOURCE MANAGEMENT AND JOB SCHEDULING

The general process of PIKA consists of four major steps detailed in Algorithm 1. At the beginning of each slot, the steps aim at finding the CPU/RAM utilization of a PM which exceeds its capacity. The gap value is calculated at step 2, the broker makes the decision of either proactive consolidating or switching-on OFF PMs in near future. Then, the broker updates the slack of each batch job and places all the web jobs and the selected batch jobs in this slot. The consolidation decision at step 4 is based on the renewable energy availability.

A. Overloaded PM detection

Since we introduce the resource over-commit policy in PIKA, the varying CPU/RAM load of the jobs at each slot may lead the PMs to overload situations. At the beginning of each slot, if a PM utilization exceeds its capacity, the broker first verifies the batch jobs utilization on this PM and suspend them if necessary.

Suspend the corresponding batch jobs, If U batch ≥ U exceed . Push the web job to the waiting queue, otherwise.

(2) We first consider suspending the current running batch job on overloaded PM until the PM returns to normal state. If all the batch jobs on overloaded PM has been suspended and the PM still exceeds its capacity, the broker migrate the web jobs to other ON PMs. Assume, for example, there is a PM which exceeds it CPU capacity. Firstly, the broker sorts the VMs in ascending order of CPU utilization. Then, the broker pushes the web job with largest CPU utilization to the web job waiting queue until the PM state becomes normal. If the PM backs to normal state, the broker finishes the overload detection process. The queued web jobs are scheduled onto other ON PMs before any new allocation. The suspended batch jobs are pushed to the batch job pool, then they are either migrated to other ON PMs if possible or put in the waiting list of the next slot to run.

B. Gap

At the beginning of each slot, the broker receives the total electricity generated from renewable energy for the next slot. Ideally, if there is sufficient renewable energy to cover all the electricity need of all already executed jobs, the broker keeps the current PM state. Otherwise, the broker proposes a plan for VM consolidation in order to decrease the number of ON PMs. The metric gap(t) is the key module of PIKA to aid the broker to dynamically adjust the number of ON PMs. The gap(t) (line 6) describes the difference between the renewable energy and current energy consumption in one slot. It is calculated as follows:

gap(t) = E r (t) - i E s on i (t) (3)
Where E r (t) denotes the predicted amount of renewable energy at the t th slot and i E s on i (t) denotes the current energy consumption of all the ON PMs (line 3-5). According to the gap(t), the broker decides whether to consolidate at the t th slot.

If gap(t) ≥ 0 (line 8-17), it means there is extra renewable energy that may be used to switch on n OFF PMs and to run more jobs (opportunistically). Otherwise, it is necessary to switch off m PMs in order to reduce the current energy consumption. To determine the value of n, the broker first sorts the powered-off PMs list by the metric as follow:

e efficient = 1 max(Es) max(smips) = max(s mips ) max(E s ) (4)
Where s mips is the CPU performance that is presented on megahertz (MHz), the max(E s ) is the maximum energy consumption of the PM s and the max(Es) max(smips) denotes energy consumption per unit of MHz. The lower the value of max(Es) max(smips) , the more efficient the PM is. The higher is the value of e efficient , the more efficient it is. Each time the broker sorts the OFF PMs list in decreasing order and selects the n first PMs to switch on where n is equal to count (line 27).

When gap(t) < 0 (line 18-26), there is no sufficient renewable energy to cover all the ON PMs energy consumption. It necessitates performing a consolidation to decrease the number of ON PMs. Due to VM migration also having an energy overhead, we need to minimize the number of migrations. Thus, the broker seeks a set of PMs which has fewer number of VMs. The size of this set is equal to the value of count (line 27).

C. VM placement

As described in Section II, by means of resource overcommit policy, the VM placement problem is transformed from 2-dimensional to 1-dimensional bin packing problem. In 1dimensional bin packing problem, FFD (First Fit Decreasing) is a classic greedy algorithm which is proved to use: maximum 11/9 ×n + 1 bins where n presents the number of bins in the optimal solution [START_REF] Yue | A simple proof of the inequality F F D(L) ≤ 11 9 OP T (L) + 1, ∀L for the FFD bin-packing algorithm[END_REF]. The FFD and over-commit policy are combined in PIKA to optimize the resource allocation, the broker first places all the web jobs and then places all or a part of the batch jobs depending on the remaining ON resources. If there are some batch jobs that cannot successfully be scheduled in this step, the broker updates their slack time and reschedules them in a later slot.

But if there is no sufficient resources to place all the web jobs at the first step, the broker suspends a part of or all the batch jobs which are executed on current ON PMs. Then the broker places the web jobs on the current ON PMs. If it still cannot meet all the web jobs resources requirement, the broker actives one or more OFF PMs directly.

D. Consolidation and migration

Through the above analysis presented in Section V-C, we detail our heuristic for the problem of dynamic consolidation in context of renewable energy. We subdivide the consolidation problem into the two following issues: when and how to consolidate.

1) When to consolidate:

There is not enough renewable energy to cover all the servers energy consumption for the current slot, i.e., the count is negative. In this way, the broker attempts to decrease the number of current ON PMs if possible.

2) How to consolidate: The broker gets the number of PMs that should be switched-off via count. The broker seeks the PM which has the least number of web jobs and tries to migrate all the web jobs from this PM to others. While a PM is selected to be switched-off, the broker tries to migrate all the web jobs from this PM to the others and interrupts all the batch jobs on this PM. The affected batch jobs are pushed into the batch pool and wait for the broker to complete the web job placement. If there is not enough free resources on other ON PMs to migrate all the VMs on this PM, the broker does not perform the migration and aborts the consolidation process. If each web job on this PM can find an another ON PM to host it, the broker migrates all of them and pushes all the batch jobs on this PM to the batch pool. Once all the VM migrations have been completed, the broker switches off this PM. The broker repeats this process until it reaches the expected number count if possible.

VI. DIFFERENT OVER-COMMIT RESOURCE POLICIES

As stated in Section II, the resource over-commitment is an efficient method to increase the CPU/RAM utilization in order to minimize the number of ON PMs. Figure 1 provides an analysis on the real world workload traces. The VMs average CPU utilization is 15% and 50% of RAM utilization in comparison with their original VM resource requirement. It leads to a significant wastage of the CPU/RAM resources. The over-commitment can be used to further optimize the utilization of PMs. We define four resource over-commitment policies, which cover all the possible cases.

A. Non over-commit policy

Basic case, the resources actually allocated to the VMs as requested initially.

 

 i v cpu i ≤ S cpu n , ∀v i ∈ S n . i v ram i ≤ S ram n , ∀v i ∈ S n . (5)

B. Over-commit RAM policy

In this case, the broker only guarantees the VM CPU resource requirement and over-allocate the RAM.

v cpu i ≤ S cpu n , ∀v i ∈ S n . (6)

C. Over-commit CPU policy

In contrast with the Over-commit RAM policy. The overcommit CPU policy over-commits CPU resource instead of RAM resource. The broker considers the RAM capacity of server and there is no upper-bound of CPU resource.

v ram i ≤ S ram n , ∀v i ∈ S n . (7)

D. Optimal Over-commit CPU/RAM policy

The optimal over-commitment policy that involved both CPU and RAM resources. This solution needs to analyze the history of each job and predict the job resource utilization for the near future. It is designed to keep the PM utilization always close to the upper bound. However, there is a risk that may lead the PM be overloaded if the prediction is not accurate. Consequently, a huge number of migrations may occur and lead to an additional energy consumption and performance degradation.

VII. EXPERIMENTAL SETUP

A. Trace-driven simulator

To evaluate the proposed algorithms in PIKA framework under different multiple-parameters configuration, we built a novel trace-driven simulator. The simulator is modeling a single energy-aware data center. It allows to simulate different resource allocation and scheduling policies.

B. Job workload

We use the trace as stated previously in Section III-B and use it for all the simulations in this paper. In this trace, each job consists of job id, time stamp, initial VM resource requirement, the Instantaneous CPU and RAM utilization. To eliminate the noise of CPU and RAM utilization for each job, the CPU and RAM utilization is averaged over the interval T (T = 1 slot) As shown in Table I, we extract a non-holiday week. It contains 750 jobs per day including 150 web jobs and 600 batch jobs. The slack time of a web job is defined as less than a slot. The slack time of batch job is defined as double time than its execution time. Half of jobs are submitted at anytime before noon (0h -12h), and the other half are submitted from noon to night (12h -24h).

C. Solar energy traces

To build the renewable energy workload in the simulator, we use the database provided by the University of Nantes, France. This database records the solar power data every hour (24 values per day). We choose the trace from the same week as the trace we used for the jobs.

VIII. EVALUATION

A. The performance of resource over-commitment policy

We implement the aforementioned resource overcommitment policies on two PMs. PM-1 consists of 12 cores and 48GB RAM (2933 MHz each core). PM-2 consists of 16 cores and 192 GB RAM (2933 MHz each core). We submit the same number of jobs for the two PMs. We combine these resource over-commitment policies with FFD algorithm as mentioned in Section V-C.

Table II and Table III show the number of hosted VMs for different resource over-commitment (OC) compared to the , this ratio is > 1). Besides for PMs, the ratio are 1/4 and 1/12 respectively (PM1 and PM2). When the broker maps the VMs to the PMs, the total CPU resource requirements of VMs will reach the PM's CPU upper bound before the RAM reaches its upper bound. Therefore, the result of CPU OC it better than the NO OC and RAM OC policies. However, the OPT outperforms up by 50% in comparison with CPU OC. Recall that, OPT dynamically adjusts the over-commitment threshold under assumption with high accurate VM utilization predictor. Figure 4 shows the CPU over-commit ratio and the mean value of real PM CPU load of CPU OC policy. The mean value of real CPU load of PM is between 45% and 70%. It over-commits near 4.1 times more the CPU resource than the PM original CPU capacity over 168 hours. Figure 5 shows the RAM over-commit ratio and the mean value of real PM RAM load of CPU OC policy. The mean value of real CPU load of PM is between 34% and 44%. It does not over allocate the RAM resources and the mean value of real PM RAM load is near 95%. The above Figures 4 and performance without a complete knowledge of future and rarely leads a PM to be overloaded. For this reason, we chose to implement CPU OC into PIKA.

B. Energy model

1) Physical machine model: The variation in energy consumption of a PM mainly depends on CPU utilization [START_REF] Orgerie | A Survey on Techniques for Improving the Energy Efficiency of Large-scale Distributed Systems[END_REF]. We experiment multiple tests on Taurus node at the Lyon site of Grid'5000, a large-scale and versatile test-bed for experimentdriven research. Each node has 12 cores and each core presents 8.3% overall CPU utilization. We use stress benchmark to vary CPU utilization in order to estimate the server's energy consumption according to its load. We activate one more core each time and keep the core at 100% utilization within 300 seconds. The test begins with 0 core (idle state). Figure 6 presents the energy consumption with different numbers of active cores. Note that an idle PM (0% CPU utilization) consumes about half of full charged PM energy consumption (100% CPU charge). Similar observations can be found in literature [START_REF] Quesnel | Estimating the Power Consumption of an Idle Virtual Machine[END_REF]. We model the PM energy consumption as successive steps. The total CPU resource is divided into 12 intervals (i.e. number of cores). The number of cores transforms on CPU utilization (i.e. 0 core = 0%, 1 core = 8.33%, 2 cores = 16.66% and so on till 12 cores = 100% CPU utilization). We create a sequence of number corresponding to the cores. The real PM utilization is easily falling into an interval of two consecutive numbers a, b (a < b) from this sequence. The PM energy consumption is the difference between a's energy consumption and b's energy consumption. Then, we construct a linear model between the two numbers a, b to calculate the PM energy consumption, this model is similar to the one used in [START_REF] Beloglazov | Energy efficient resource management in virtualized cloud data centers[END_REF].

E n = E a + (E b -E a ) × u n -a b -a , where a < u n < b (8)
where E a denotes the lower bound a energy consumption and E b denotes the upper bound b energy consumption. u n present the PM CPU utilization and un-a b-a denotes the percentage of difference between PM utilization and the lower bound. Moreover, there is a switching cost, as the broker needs to dynamically switch on PMs or switch them to sleep mode. We define a fix energy overhead in the later simulation when broker switches on a PM, this value has been measured experimentally.

2) VM energy consumption model: The VM energy consumption model consists of two modules: VM creation energy overhead and VM migration overhead. The VM creation energy overhead is defined as a fixed value in the simulator. The energy consumption of VM migration depends on the following parameters: migration duration time and the energy consumption per unit of time. The study of [START_REF] Huang | Power consumption of virtual machine live migration in clouds[END_REF] shows that the duration of live migration depends mostly on the memory and disk used by the migrated VM. The VM migration will sightly increase the CPU utilization on destination server. As abovementioned, the energy consumption of PM increases almost linearly with CPU utilization, we formalize the linear model for energy consumption of VM live migration as following:

E M igration = C RAM + C DISK B × E CP U (9)
where B denotes the bandwidth between the PMs. C RAM + C DISK denotes the migration overhead that is the sum of the memory size and the disk size of a VM , C RAM +C DISK B denotes the duration of migration and E CP U is the extra CPU energy consumption on destination PM per unit of time.

C. Energy consumption

The result of energy consumption for both baseline algorithm and PIKA are shown in Figure 7. The top curve presents the baseline result (with FFD algorithm to allocate the VMs) and the bottom corresponds to PIKA. The energy consumption of baseline is flat. The workload scheduling is not affected by the variable renewable energy supply (the green curve). The energy consumption of PIKA is following the renewable energy variations. PIKA significantly increases the renewable energy integration into data center. While the renewable energy becomes unavailable, the broker switches off some ON PMs and only launches some essential jobs (web job, mutated batch job and a few of batch jobs on remaining resources). Due to this behavior, PIKA finishes all the job 11 hours later than the baseline. Indeed, in PIKA, the broker opportunistically schedules the batch jobs. So, some batch jobs are delayed when there is no sufficient renewable energy.

D. Simulation results

Table IV shows the result of brown-, renewable-and totalenergy consumption for the baseline and PIKA. Compared to the baseline, PIKA reduces by 44.9% brown energy consumption and increases by 110.1% the renewable energy integration. The results show that PIKA significantly reduces the brown energy consumption in compare with the baseline, representing a typically energy-efficient algorithm (but not renewable-aware). The results also indicate that PIKA consumes 31% more energy in total. This is because PIKA performs dynamic VM consolidation to adjust the number of ON PMs and that leads to a large number of VM migrations compared with baseline (the migration in baseline is only in case of overloading PM). But all of this extra energy consumption comes from renewable energy supply. So, this extra energy is not used and thus wasted in the baseline case.

This work shows the opportunity created by medium-sized data centers partially powered by renewable energy in order to save energy for distributed Cloud infrastructures, such as the one promoted by the Future Internet.

IX. CONCLUSION AND FUTURE WORK

Data centers partially powered by renewables become attractive for the new generation cloud architectures. It significantly reduces the traditional energy consumption and CO 2 footprint. This paper is dealing with resource allocation and opportunistic job scheduling in a small/medium mono-site data center. Our proposal is called PIKA and the preliminary results demonstrate that it can outperform the classical energyefficient VM management algorithms.

Our next work involves studying the workload prediction, this method may still optimize the resource utilization to further reduce the brown energy consumption. Furthermore, we plan to integrate the management of the cooling system within PIKA in order to jointly and more efficiently tackle VM placement and job scheduling.

Fig. 1 .

 1 Fig. 1. CPU and RAM real utilization over one-week of real trace

Fig. 2 .

 2 Fig. 2. PIKA framework

Fig. 3 .

 3 Fig. 3. VM placement for batch and web

Algorithm 1

 1 General process of PIKA1: Step 1: Detect the overloaded PMs; 2: Step 2: Launch calculating gap process and make the decision of whether consolidate depending on the count; 3: Step 3: Update each batch job's slack;4: Step 4: Select the jobs that can be executed in the current slot and place them adequately on PMs (VM placement algorithm). The decision of VM consolidation relies on the result of 2nd step.

Fig. 4 .

 4 Fig. 4. (a) CPU real-time utilization; (b) total VM allocation ratio

Fig. 5 .

 5 Fig. 5. (a) RAM real-time utilization; (b) total RAM allocation ratio

Fig. 6 .

 6 Fig. 6. Energy consumption of node Taurus at different number of cores in Watts on grid'5000 Lyon site. It owns 12 cores x 2933 MHz

Fig. 7 .

 7 Fig. 7. Energy consumption : baseline vs PIKA

  The gap(t) function is shown in Algorithm 2.

	Algorithm 2 gap function
	2: pmList.splitTwoSubList()→powerOnPmList	and
		powerOffPmList
	3: for each pm ∈ powerOnPmList do
	4:	sumPmPower += CurrentPower(pm);
	5: end for
	6: gap = sumGreenPower -sumPmPower
	7: count=0;
	8: if gap >= 0 then
	9:	powerOffServerList.sortByServerScore();
	10:	for each server ∈ powerOffServerList do
	11:	while gap > 0 do
	12:	gap -= pm.getMaxPower();
	13:	pm.setPotentialOn(true);
	14:	count++;
	15:	end while
	16:	end for
	17: else
	18:	powerOffServerList.sortByNumVM();
	19:	for each server ∈ powerOnServerList do
	20:	while gap < 0 do
	21:	gap += pm.getCurrentPower();
	22:	pm.setPotentialOff(true);
	23:	count--;
	24:	end while
	25:	end for
	26: end if
	27: return count;

1: INPUT: pmList, renewable energy availability in the current slot: sumGreenPower OUTPUT: the number of PMs to be switched-on or -off

  The first inequation denotes the PM n can simultaneously host multiple VM i if and only if the both CPU/RAM resource request of VM can not exceed the PM n's capacity.

	Where VM v cpu i	i	v cpu i on PM S cpu denotes the CPU resource request of all the n .

TABLE I

 I

		.	JOB CHARACTERISTICS (HOUR)
	Type	Number of jobs	Execution time	Slack time
	Web		150	24	< 1
	Batch		600	6	12

TABLE II

 II

	.	NUMBER OF VMS	TABLE III.	NUMBER OF VMS
		PM-1			PM-2
	Policy	Mean	St.dev.	Policy	Mean	St.dev.
	NO.OC	6.356	2.154	NO.OC	8.141	5.337
	RAM.OC	6.724	2.651	RAM.OC	11.451	2.656
	CPU.OC	27.035	2.744	CPU.OC	44.893	4.515
	OPT	60.994	9.775	OPT	100.441	5.904
	solution solution combining both over-commitment techniques
	(OPT). On PM1, the first two policies (NO OC and RAM OC)
	have nearly identical performance. CPU OC hosts more than 3
	times more VMs than both NO OC and RAM OC policies . On
	PM-2, CPU OC also offers 400% better performance than both
	NO OC and RAM OC. This can be explained by an analysis
	of the trace: for the vast majority of VMs, the CPU resource
	requirement is higher than the RAM resource requirement (i.e.
	number of CPU cores
	RAM (GB)

TABLE IV .

 IV ENERGY SAVING RESULTS(KW)

	Algorithm	Total E. C.	Brown E. C.	Renewable E. C.
	Baseline	513.633	259.559	254.073
	PIKA	676.895	142.957	533.938
		31% ↑	44.9%↓	110.1%↑

ACKNOWLEDGMENT

This work has received a French state support granted to the CominLabs excellence laboratory and managed by the National Research Agency in the "Investing for the Future" program under reference Nb. ANR-10-LABX-07-01.