Bradley J Kavanagh

Re-examining the significance of the 750 GeV diphoton excess at ATLAS

The excess seen in the diphoton channel at around 750 GeV by both ATLAS and CMS has caused a great deal of excitement in the particle physics community. However, there has recently been much discussion about uncertainties in the significance of the peak seen by the ATLAS experiment. In this note, we aim to estimate this significance using a range of possible parametrisations for the smooth diphoton background. We obtain a local significance close to that reported by ATLAS and further demonstrate that the significance of the excess is not substantially reduced when more complicated background functions are considered. In particular, the background contribution is strongly constrained by the small numbers of events at large diphoton invariant mass. Future data releases will improve constraints on the diphoton background, as well as clarifying the true nature of the 750 GeV excess.

I. INTRODUCTION

In the first results from Run-II of the LHC, both the ATLAS and CMS collaborations reported an excess above the expected background in the diphoton channel [START_REF] Kavanagh | Search for resonances decaying to photon pairs in 3.2 fb -1 of pp collisions at √ s = 13[END_REF][START_REF]Search for new physics in high mass diphoton events in proton-proton collisions at 13TeV[END_REF] at a diphoton invariant mass of m γγ ∼ 750 GeV. The local significance of this excess was reported as 3.9σ (2.3σ global) by ATLAS and 2.6σ (1.2σ global) by CMS. The possibility that this excess could be due to a new bosonic particle has caused a great deal of excitement in the particle physics community and the list of papers exploring the implications for New Physics is of course too long to mention.

Since the announcement of these results [START_REF]CMS physics results from Run 2[END_REF], a number of authors have performed fits to the Run-II diphoton data from a theory perspective in order to estimate the mass, width and cross section times branching ratio of a possible new particle (see e.g. Ref. [START_REF] Falkowski | [END_REF]). However, these analyses typically assume a fixed form for the diphoton background whose parameters are not included in the fits. This procedure is not strictly correct, as parameter values extracted from the data should take into account uncertainties in the background. A full likelihood analysis of both the ATLAS and CMS binned diphoton data was performed in Ref. [START_REF] Buckley | [END_REF], assuming the same background models used by the respective collaborations. The local significances reported therein are in broad agreement with those reported by ATLAS and CMS.

However, a recent paper by Davis et al. [START_REF] Davis | [END_REF] explored the possibility of allowing more freedom in the modelling of the ATLAS continuum diphoton background. Their conclusion was that the background could be fit by a number of different functional forms and that some of these forms lead to a lower significance for the diphoton excess than that reported by ATLAS, perhaps as low as 2.0σ. Motivated by this discrepancy between the results of the AT-LAS collaboration and those of Davis et al., we perform an independent analysis of the binned ATLAS diphoton * bradley.kavanagh@lpthe.jussieu.fr invariant mass spectrum to estimate the significance of the excess under different background assumptions. In particular, we perform a Poissonian likelihood fit over 40 bins in the diphoton invariant mass m γγ , making this analysis distinct from that of Davis et al. The code used to perform this analysis is publicly available online [7].

The intention of this note is not to improve upon the ATLAS analysis or to report definitive values for the significance of the 750 GeV excess. Indeed, this is not possible with the information which is publicly available. Instead, we aim to explore the robustness of the excess to variations in background fitting, as well as to clarify certain aspects of the fits performed by ATLAS.

The different functional forms we consider for the diphoton background are given in Eq. 2, while the local significances we obtain for the 750 GeV diphoton excess (using each functional form) are given in Table III. We describe in detail the fitting procedure used in this work in Sec. II, followed by the results of background-only and signal + background fits in Sec. III and Sec. IV. Finally, in Sec. V we discuss the implications of these results for the significance of the 750 GeV excess.

II. FITTING PROCEDURE

Within a frequentist framework, the significance of a possible excess in the diphoton spectrum can be estimated by calculating the maximum likelihood obtained assuming that only backgrounds contribute to the observed events (background-only) and that obtained assuming that both signal and background have a contribution to the observed event rate (signal + background). By comparing the maximum likelihood obtained under these two assumptions, we can then calculate the significance of a possible signal in the data.

The likelihood L(θ) is simply the probability of obtaining the observed data assuming a particular set of parameters θ. For binned data, this has the Poisson form:

L(θ) = n bins i=1 (N i e ) N i o N i o ! exp(-N i e ) , (1)
where N i e = N i e (θ) is the expected number of events in the i th bin (for a given set of parameters under the background or signal + background hypothesis) and N i o is the observed number of events in the i th bin. In particular, we do not include any additional terms in the likelihood due to uncertainties on the diphoton invariant mass resolution (as in Ref. [START_REF] Kavanagh | Search for resonances decaying to photon pairs in 3.2 fb -1 of pp collisions at √ s = 13[END_REF]).

We use a total of n bins = 40 bins in the m γγ invariant mass, each with width 40 GeV, spanning the range m γγ ∈ [150, 1750] GeV. We obtain the number of events observed in each bin from Ref. [START_REF] Kavanagh | Search for resonances decaying to photon pairs in 3.2 fb -1 of pp collisions at √ s = 13[END_REF] (by digitising Fig. 1 therein using a publicly available plot digitiser [START_REF] Rohatgi | Webplotdigitizer v3.9[END_REF]). This procedure may introduce some digitisation error into the analysis, particularly for bins which contain a large number of events. However, we have explicitly checked the impact of digitisation errors on our results and, as described in Sec. V, these do not affect the conclusions we report.

The expected number of background events is obtained by integrating the background event distribution over each bin in m γγ , for a given choice of functional form and for a given set of background parameters θ b . In Ref. [START_REF] Kavanagh | Search for resonances decaying to photon pairs in 3.2 fb -1 of pp collisions at √ s = 13[END_REF], the ATLAS collaboration discuss a set of possible empirical functions for the continuum background. These have been adapted from functions used in multijet searches for New Physics [START_REF] Aad | [END_REF] and have been validated against both Monte Carlo and data samples. In Ref. [START_REF] Davis | [END_REF], Davis et al. introduced another possible parameterisation for the continuum diphoton background (also validated with a Monte Carlo study) which we discuss further in Sec. V. In this work, however, we focus on those background functions used by the ATLAS collaboration, given explicitly by

f k (x) = N (1 -x 1/3 ) b x k j=0 aj (log x) j (2)
where x = m γγ / √ s and f (x) is the differential distribution of expected background events (in units of events/40 GeV). 1 We will consider k = 0, 1, 2 in this work, for which the parameters b, a 0 , a 1 and a 2 determine the shape of the background and are allowed to vary in the analysis. The parameter N which controls the overall normalisation of the background is typically included when using this class of functions [START_REF] Aad | [END_REF]. While it is not explicitly described in Ref. [START_REF] Kavanagh | Search for resonances decaying to photon pairs in 3.2 fb -1 of pp collisions at √ s = 13[END_REF], the free normalisation parameter has been included in previous analyses of the diphoton channel [10] and its inclusion will allow us to better match the results reported by ATLAS. In light of this, we will assume that ATLAS allow for a free normalisation of the I. Prior ranges for the background and signal parameters used in this analysis. In each case, the prior is constant over the specified range. We note that because we use MCMC only to find the maximum likelihood point, the exact form of the prior does not affect the results presented.

background in their analysis, though for comparison we will consider cases where N is set to 1 ('fixed-N ') as well as where N is allowed to vary ('free-N ').

The expected number of signal events coming from a possible new resonance near m γγ ∼ 750 GeV is obtained using either the narrow width approximation (NWA) or by allowing the width of the resonance to vary freely (free-width). Under the NWA, the signal is modelled as a Gaussian centred on m γγ = m B with standard deviation σ given by the diphoton invariant mass resolution. This resolution is taken to increase linearly from σ = 2 GeV at m γγ = 200 GeV to σ = 13 GeV at m γγ = 2 TeV [START_REF] Kavanagh | Search for resonances decaying to photon pairs in 3.2 fb -1 of pp collisions at √ s = 13[END_REF], giving a value of around σ = 5.4 GeV near the putative resonance at 750 GeV. This is a rather crude approximation to the double-sided Crystal Ball function (DCSB) used by the ATLAS collaboration in their NWA analysis. However, given that the bin widths used in this analysis are much larger than the mass resolution, we do not expect the precise shape of the narrow resonance to strongly affect our results. Indeed, we are still able to recover a significance close to that reported in Ref. [START_REF] Kavanagh | Search for resonances decaying to photon pairs in 3.2 fb -1 of pp collisions at √ s = 13[END_REF]. The signal is normalised so as to contribute a total of N S events, meaning that for the NWA we have two free parameters describing the signal: θ s = (m B , N S ).

In the free-width analysis, we model the signal as a Breit-Wigner distribution with width Γ = αm B . We do not include the effects of finite diphoton invariant mass resolution in the free-width analysis. However, as in the NWA analysis, we can still recover a significance close to that reported in Ref. [START_REF] Kavanagh | Search for resonances decaying to photon pairs in 3.2 fb -1 of pp collisions at √ s = 13[END_REF]. In the free-width analysis, we have three free parameters:

θ s = (m B , N S , α).
In order to explore the signal and background parameter spaces, we use the publicly available Markov Chain Monte Carlo (MCMC) software emcee [11], which uses Affine Invariant MCMC Ensemble sampling [12]. We perform at least 20000 likelihood evaluations for each fit. The prior ranges for the parameters used in this work are listed in Table I, though we emphasise that emcee was only used to find the maximum likelihood points and therefore that the precise form of the priors should have little impact on the analysis.

We quantify the significance of a possible excess by constructing the test statistic q 0 [13]:

q 0 =        -2 log L(N S = 0, θb ) L( NS , m B , α, θb ) NS ≥ 0 0 NS < 0 , (3)
where θb is the set of background parameters which maximise the likelihood under the background-only hypothesis and θ is the set of parameters which maximise the full likelihood under the signal + background hypothesis for a given value of m B and (if we are considering the free-width analysis) α. Here, we restrict to positive parameter values for the number of signal events, so the best fit value NS will always be positive. The local p-value for the background-only hypothesis is obtained from the observed value of the test statistic q 0,obs using,

p 0 = ∞ q 0,obs f (q 0 ) dq 0 ( 4
)
where f (q 0 ) is the probability density function of q 0 under the background-only hypothesis. According to Wilks' Theorem [14], the log-likelihood ratio (appearing in the top line of Eq. 3) is asymptotically χ 2 m distributed if the background-only hypothesis is correct. The number of degrees of freedom m is given by the difference between the number of parameters in the background-only and the signal + background hypotheses. For a fixed value of m B (and α), there is only one free parameter (the number of signal events N S ) so we have m = 1. As detailed in Ref. [13], in this case q 0 follows a 'half chi-square' distribution and the local signal significance is simply given by Z = √ q 0,obs . The maximum value of the local significance can then be obtained by maximising over m B (or m B and α in the free-width analysis).

III. BACKGROUND-ONLY FITS

We begin by examining in Fig. 1 the fits to the ATLAS data (black points) which are obtained when no signal contribution is included. We show the best fit (maximum likelihood) background curves for k = 0 (red) and k = 1 (blue), both using a free normalisation N (solid lines) and when setting N = 1 (dashed lines). We note that the error bars on the data points are for illustrative purposes and denote the 1σ confidence intervals on the mean number of expected events in each bin given the number of observed events.

As pointed out by Davis et al. [START_REF] Davis | [END_REF], the k = 0 fixednormalisation background curve (dashed red) appears to underestimate the background above m γγ > 1000 GeV. It is perhaps not surprising that this parametrisation is not able to give a good fit to the background over the entire invariant mass range, given that it includes only two 200 400 600 800 1000 1200 1400 1600 Events / 40 GeV Background-only fit

m γγ [GeV]
n bins =40 k =0 k =1
Free norm. Fixed norm.

FIG. 1. Background-only fits to the ATLAS diphoton invariant mass spectrum. The observed numbers of events in each bin are shown as black circles, while the curves show the background distribution for the k = 0 (red) and k = 1 (blue) empirical functions defined in Eq. 2. We show background fits with both free normalisation (solid) and with fixed normalisation (dashed). The background parameters used are those which maximise the likelihood for the background-only hypothesis. Note that the k = 2 background fits are not shown as they lie close to the k = 1, free normalisation curve (solid blue). free parameter. When adding an additional parameter to the background function, either using k = 0 with free normalisation (solid red) or k = 1 with fixed normalisation (dashed blue), the fits tend to prefer a higher background at high mass in order to alleviate this possible tension. The background contribution is also increased in the region of the 750 GeV excess, suggesting that the significance of the peak will be reduced when including these additional parameters.

We note that both the k = 0, free-N background and the k = 1, fixed-N background appear to match closely the background-only fit reported in Fig. 1 of Ref. [START_REF] Kavanagh | Search for resonances decaying to photon pairs in 3.2 fb -1 of pp collisions at √ s = 13[END_REF]. The ATLAS collaboration state that they use the k = 0 background function in Ref. [START_REF] Kavanagh | Search for resonances decaying to photon pairs in 3.2 fb -1 of pp collisions at √ s = 13[END_REF] and we therefore assume that they also fit the normalisation N of the background, in order to match the results shown here in Fig. 1. It therefore appears then that the background curve used by ATLAS fits the data well, apart from in the region of the 750 GeV excess and above m γγ ∼ 1600 GeV (where ATLAS report an excess with 2.8σ local significance).

Following Ref. [START_REF] Davis | [END_REF], we have calculated the Bayesian Information Crierion (BIC) [15], a model selection criterion which be used to compare the fit to data obtained using different models, penalising models which have additional free parameters. The Akaike Information Criterion (AIC) [16] is a related model selection criterion, though in general the BIC penalises the addition of extra parameters more strongly. The BIC is defined as: The BIC is defined in Eq. 5. The background function used by the ATLAS collaboration in Ref. [START_REF] Kavanagh | Search for resonances decaying to photon pairs in 3.2 fb -1 of pp collisions at √ s = 13[END_REF] is marked with a dagger.

BIC = -2 log L + N p log n bins , ( 5
)
where L is the maximum likelihood and N p is the number of free parameters in the model. The lower the BIC, the stronger the evidence for the model in question. A difference of around 2 in the BIC between two different models is considered positive evidence that the model with the higher BIC should be rejected in favour of the one with the lower. A difference in BIC greater than 6 is considered strong evidence [17].

We show in Table II the maximum log-likelihood and the value of the BIC obtained in background-only fits using the three parametrisations k = 0, 1, 2 with both free and fixed normalisation. We find that the k = 0, fixed-N background gives the largest BIC, indicating that there is strong evidence to reject this background in favour of the alternative parametrisations. Note that the remaining BIC values are all rather close in value, indicating that there should be no strong preference amongst them. This matches the conclusion of the ATLAS collaboration [START_REF] Kavanagh | Search for resonances decaying to photon pairs in 3.2 fb -1 of pp collisions at √ s = 13[END_REF], which used both a Fisher test and a 'spurious signal' analysis2 to determine that no background models more complex than k = 0, free-N were necessary to fit the data.

We finally note that increasing the number of parameters further, with the k = 2 functional form, does not substantially improve the fits, as reflected by the slightly larger BIC value for k = 2, free-N . We do not show these functional forms in Fig. 1, because they lie very close to the solid blue k = 1, free-N curve. The smooth background at high m γγ cannot be increased further due to the tension with a number of empty bins. It therefore seems difficult to further increase the background 200 400 600 800 1000 1200 1400 1600 FIG. 2. Signal + background fits to the ATLAS diphoton invariant mass spectrum, allowing the width of the resonance near 750 GeV to vary freely. The colour scheme matches that of Fig. 1. The parameters used are those which maximise the likelihood for the signal + background hypothesis. As in Fig. 1, we do not show the fits for k = 2, as these lie close to the k = 1, free normalisation curve (solid blue).

m γγ [GeV]
contribution to the rate at either m γγ ∼ 750 GeV or m γγ ∼ 1600 GeV by going to more complex smooth background parametrisations.

IV. SIGNAL + BACKGROUND FITS

We now examine the signal + background fits to the ATLAS diphoton spectrum and determine the significance of the diphoton excess under different assumptions for the background. In Fig. 2, we show the maximum likelihood fits to the data for the signal + background hypothesis using the free-width analysis for the possible new resonance near 750 GeV. As before, we include fits for k = 0, 1 (see Eq. 2) using both free and fixed normalisation, with the colour scheme matching that of Fig. 1.

Adding a new resonance near 750 GeV clearly allows the data to be better fit with all background parametrisations. As in the background-only fits, the k = 0, free-N and k = 1, fixed-N backgrounds produce similar results. In all cases, the expected rate at high m γγ is reduced compared with the background-only fits. With the 750 GeV excess saturated by the signal contribution, the background above 1000 GeV can be reduced to better fit the large number of empty bins.

In Table III, we report the maximum local significance of the ATLAS 750 GeV excess obtained in this work, using each of the functional forms introduced in Eq. 2. For reference, the significances reported by ATLAS in Ref. [START_REF] Kavanagh | Search for resonances decaying to photon pairs in 3.2 fb -1 of pp collisions at √ s = 13[END_REF] are 3.6σ (NWA) and 3.9σ (free-width). The significances obtained in this work, using the same background function used by ATLAS, are 3.4σ (NWA) and 3.6σ (free- GeV diphoton excess obtained in this work using each of the background functions described in Eq. 2, assuming a freely varying resonance width (free-width) and under the narrow width approximation (NWA). The background function used by the ATLAS collaboration in Ref. [START_REF] Kavanagh | Search for resonances decaying to photon pairs in 3.2 fb -1 of pp collisions at √ s = 13[END_REF] is marked with a dagger. For comparison, we also give the local significance reported by the ATLAS collaboration.

width). We discuss in Sec. V the possible sources of this small discrepancy.

For the simplest background function (k = 0, fixed-N ) the significance of the excess is significantly larger. 3 However, in Sec. III, we found (in agreement with ATLAS) that this simple function fits the background significantly worse than the more complex ones. The remaining background parametrisations lead to very similar significances compared with the k = 0, free-N function. Thus, the significance of the excess appears to be robust against the choice of background.

V. DISCUSSION

As initially pointed out by Davis et al. [START_REF] Davis | [END_REF], the significance of the 750 GeV diphoton excess reported by AT-LAS is higher when the background is fit with the k = 0, fixed-normalisation parametrisation and is reduced when an extra parameter is added to the background fit. However, this appears to be consistent with the results reported by ATLAS. Assuming that the overall normalisation of the k = 0 background is also included in the fit, we recover the background-only fit presented in Ref. [START_REF] Kavanagh | Search for resonances decaying to photon pairs in 3.2 fb -1 of pp collisions at √ s = 13[END_REF] and obtain significances close to (but slightly smaller than) those reported by the ATLAS collaboration.

We also note that the data show no preference for an increase in complexity of the background function (as demonstrated by the Bayesian Information Criterion study in Table II). Furthermore, we find that adding additional parameters to the background fit does not have any impact on the significance of the excess. This is because of a number of bins above m γγ ∼ 1100 GeV which 3 We note that we find an even larger significance using this background than that reported by Davis et al. [START_REF] Davis | [END_REF], though given the different statistical approaches, this is perhaps not surprising.

see no events. Any smooth background is constrained not to overshoot these bins. Davis et al. introduce a different possible parametrisation for the background (which was also validated by a Monte Carlo study) and find that the significance of the excess is further reduced with respect to the k = 1, fixed-N case. However, the empty bins at high m γγ were not included in that analysis, leading to a background fit which overestimates the high m γγ event rate. Indeed, using the Davis et al. background parametrisation (with free normalisation) in this analysis gives a local significance of 3.8σ for a free-width resonance. This does not discount the possibility that exploring a wider range of possible background functions may impact the significance of the 750 GeV excess, but the correct constraints from the entire range of m γγ should be taken into account.

It is of course necessary to point out that the significances we report are only estimates and care must be taken when comparing with the official ATLAS analysis. In particular, the results reported by ATLAS use the full unbinned data set, while we consider here only binned data. This has the largest impact on our NWA results. The diphoton invariant mass resolution (∼5.4 GeV) is substantially smaller than the bin width (40 GeV), meaning that a narrow resonance cannot be resolved in the binned data. The fitted signal resonance is effectively widened by integration over a single bin, improving the fit to the relatively wide observed excess (∼3 bins). This means that the significance we obtain for the NWA is larger than it would be using the unbinned data. However, the ATLAS analysis includes an additional nuisance parameter for the invariant mass resolution, which we have omitted. This mass resolution parameter can be increased to improve the fit to the data, leading to a significance similar to that obtained here.

We also emphasise that the binned data was obtained by digitising the results released in Ref. [START_REF] Kavanagh | Search for resonances decaying to photon pairs in 3.2 fb -1 of pp collisions at √ s = 13[END_REF]. However, we have investigated the possible impact of digitisation error on our analysis. In order to do this, we added random noise to the first 10 bins in m γγ (distributed uniformly between -3% and +3% of the number of events in each bin) in order to simulate digitisation errors. 4 Ten such 'randomised' datasets were generated and the peak significance for each was calculated assuming the k = 0, free-N background. The resulting local significances were in the range: NWA: 3.3σ-3.5σ Free-width: 3.5σ-3.8σ .

These tests indicate that digitisation could have induced an error of order 0.2σ in the analysis, and may also explain some of the discrepancy between our results and those reported by ATLAS. We further note that in this work we have used only approximate functional forms for the signal distributions and while we approximately recover the significances reported by the ATLAS collaboration, our analysis does not capture many of the important details involved in fitting the signal and background (for example, uncertainties in the diphoton mass reconstruction). Furthermore, we have taken an empirical approach and allowed a wide range of values for the background parameters N , b, a 0 , a 1 , a 2 . This may not accurately reflect the background distributions seen in MCMC and data samples. However, restricting the possible ranges of the background parameters is likely only to increase the significance of a possible excess.

In spite of these simplifications, the broad message of this note still holds. The significance of the 750 GeV excess does not appear to be strongly affected by different choices of smooth background model. While more complicated background distributions could be explored to ensure the robustness of the significance reported by ATLAS, these are likely to be highly constrained by the lack of observed events at high diphoton invariant mass.

Note that we use the notation log for the natural logarithm log e .

The 'spurious signal' analysis requires that (for a given functional form for the background) the bias on the fitted signal yield is significantly smaller than the statistical uncertainty on the signal yield, as determined from Monte Carlo samples. See Refs.[START_REF] Kavanagh | Search for resonances decaying to photon pairs in 3.2 fb -1 of pp collisions at √ s = 13[END_REF]10] for further details.

At larger values of mγγ the small number of events means that integer numbers of events can accurately be read off the figure. However at small values of mγγ digitisation error could induce variations of order 10-100 events.

ACKNOWLEDGMENTS

The author would like to thank the authors of Ref. [6] for stimulating and helpful discussions. The help of Adam Falkowski was invaluable in the completion of this work and the preparation of this manuscript. The organisers and participants of the Magic Journal Club are also thanked for interesting (if long) discussions on the 750 GeV excess. The author is supported by the John Templeton Foundation Grant 48222 and by the European Research Council (Erc) under the EU Seventh Framework Programme (FP7/2007-2013)/Erc Starting Grant (agreement n. 278234 -'NewDark' project). The author also acknowledges the hospitality of the Institut d'Astrophysique de Paris, where part of this work was completed.


N Maculan
email: maculan@cos.ufrj.br

H Ouzia

An efficient cutting plane algorithm for the minimum weighted elementary directed cycle problem in planar digraphs

Keywords: digraph, cycle, linear relaxation, valid inequality, polytope

In this paper, we study the efficiency (both theoretically and computationally) of a class of valid inequalities for the minimum weighted elementary directed cycle problem (MWEDCP) in planar digraphs with negative weight elementary directed cycles. These valid inequalities are called cycle valid inequalities and are parametrized by an integer called inequality's order. From a theoretical point of view, we prove that separating cycle valid inequalities of order 1 in planar digraph can be done in polynomial time. From a computational point of view, we present a cutting plane algorithm featuring the efficiency of a lifted form of the cycle valid inequalities of order 1. In addition to these lifted valid inequalities, our algorithm is also based on a mixed integer linear formulation of the MWEDCP. The computational results are carried out on randomly generated planar digraph instances of the MWEDCP. For all 29 instances considered, we obtain in average 26.47% gap improvement. Moreover, for some of our instances the strengthening process directly displays the optimal integer elementary directed cycle.

Introduction

Let D = (V, A) be a connected planar digraph with V as vertex set, A as arc set and an arc weight function w : A → R. A directed cycle is a sequence (v 0 , a 1 , . . . , a k , v k ),

where k is an integer (assumed to be ≥ 1), v 0 , v 1 , . . . v k are vertices such that v 0 = v k , for every index i belonging to {1, . . . k}, a i is an arc connecting the vertices v i-1 and v i and, finally, all arcs a i have the same orientation. An elementary directed cycle is a directed cycle (v 0 , a 1 , . . . , a k , v k ) in which each vertex v i , for every index belonging to {0 . . . k}, appears once. We denote by P the polytope of all elementary directed cycles in D. That is, the convex hull of the set of incidence vectors of elementary directed cycles of the digraph D.

The minimum weighted elementary directed cycle problem consists in finding of an elementary directed cycle γ * such that:

l(γ * ) = min {l(γ) : γ ∈ C} ( 1
)
Where C is the set of all elementary directed cycles in D and l(γ) = a∈A(γ) w(a) where A(γ) denotes the arc set of the elementary directed cycle γ.

This problem is known to be NP-hard, since it includes the asymmetric traveling salesman problem as a special case. As a related work, Balas and Oosten (see [START_REF] Balas | On the cycle polytope of a directed graph[END_REF]) present a linear description of the elementary directed cycle polytope. They also introduce some facets of the polytope P. In [START_REF] Balas | On the cycle polytope of a directed graph and its relaxations[END_REF], the authors consider the dominant of the polytope P and derive others facets of P. The polytope P is also studied in undirected graphs case (see [START_REF] Coullard | On cycle cones and polyhedra[END_REF]). The integer linear formulation considered in [START_REF] Balas | On the cycle polytope of a directed graph[END_REF] and [START_REF] Balas | On the cycle polytope of a directed graph and its relaxations[END_REF] is a non compact formulation, since it contains an exponential number of constraints as the so-called multiple cycle exclusion constraints and linear ordering constraints. Other cycle polytopes with cardinality restrictions defined on graphs and digraphs are already well studied. For an exhaustive list of those publications, one can refer to [START_REF] Kaibel | On cardinality constrained cycle and path polytopes[END_REF]. Notice that all these researches investigate, in common, a partial description of the polytope P in general undirected and directed graphs. But they do not treat the resolution question.

From an algorithmic point of view, although the problem is very fundamental, the state of the art for it dates back to a seminal paper by Itai and Rodeh (see [START_REF] Itai | Finding a minimum circuit in a graph[END_REF]) from the 70s, that deals only with the unweighted variant of the problem. Itai and Rodeh present an O(n w )-time algorithm for an n-node unweighted undirected graph and an O(n w logn)-time algorithm for an n-node unweighted directed graph. In the same paper, Itai and Rodeh pose the question whether similar results exist for weighted graphs. Yuster and Zwick [START_REF] Yuster | Finding even cycles faster[END_REF] present an algorithm to find the shortest even cycle in unweighted undirected graphs in O(n 2 ). After, Lingas and Lundell [START_REF] Lingas | Efficient approximation algorithms for shortest cycles in undirected graphs[END_REF] address efficient approximation algorithms for the shortest cycles in unweighted and weighted undirected graphs. In the same way, Roditty and Williams (see [START_REF] Roditty | Minimum weight cycles and triangles: equivalences and algorithms[END_REF]) present a O(M n w )-time algorithms for directed graphs with integral edge weights in [-M, M ] (and no negative weight cycles) and for undirected graphs with integral edge weights in [1, M ], where w (less than 2.376) is the exponent of square matrix multiplication (see [START_REF] Coppersmith | Matrix multiplication via arithmetic progressions[END_REF]). On the other hand, Letchford and Pearson [START_REF] Letchford | A fast algorithm for minimum weight odd circuits and cuts in planar graphs[END_REF] give a simple O(n 3/2 logn) algorithm for finding a minimum weight odd cycle in positive weighted undirected planar graphs.

In this paper, we propose an efficient cutting plane algorithm based, in part, on the linear relaxation of a compact mixed integer linear formulation of the MWEDCP in planar digraphs. Such an algorithm possibly gives the optimal integer solution. The linear formulation of the MWEDCP considered is due to Maculan et al. (see [START_REF] Maculan | Integer linear models with a polynomial number of variables and constraints for some classical combinatorial optimization problems[END_REF]). Even in presence of negative weight cycles, at optimality, such a formulation provides the optimal integer elementary directed cycle with an additional vertex s / ∈ V and arc (s, i), i ∈ V . Generally, the structure of any optimal fractional solution of the linear relaxation of the considered MWEDCP linear formulation is such that by deleting the vertex s from the supporting sub-digraph, one disconnects it. But, each connected component of this supporting sub-digraph is actually strongly connected. From the dual digraph of every connected component of the supporting sub-digraph, we show that to any u-v elementary directed path corresponds to what we call a cycle valid inequality of order 1 (see [START_REF] Ibrahim | A new class of valid inequalities for the minimum weighted elementary directed cycle problem in digraph[END_REF]). Recall that, in strongly connected digraphs, the separation problem of cycle valid inequalities of order β is NP-hard (even in the case β = 1, (see [START_REF] Ibrahim | A new class of valid inequalities for the minimum weighted elementary directed cycle problem in digraph[END_REF]). We will see that in the case of planar digraphs, valid inequalities that correspond to u-v elementary directed paths in the dual digraph (of the supporting digraph of the optimal fractional solution) can be separated in polynomial time.

As the cycle valid inequalities of order 1 are generated w.r.t. the connected components of the sub-digraph supporting the optimal fractional solution, we therefore resort to a lifting technique to derive valid inequalities for the polytope P. The lifting technique used is close in spirit to the one presented in [START_REF] Ibrahim | Valid inequalities and lifting procedures for the shortest path problem in digraphs with negative cycles[END_REF] for the shortest path problem and applied to derive the so-called simple-lifted-valid-inequalities. However, in this work (contrary to the simple-lifted-valid-inequalities) we do not consider all the arcs of the complementary sub-digraph of the supporting digraph of the optimal fractional solution. After, we present computational results featuring the efficiency of these lifted valid inequalities. Within a cutting plane algorithm, we show that such lifted valid inequalities are capable to strengthen the linear relaxation of the evoked formulation.

The paper is organized as follows. In section 2, we present a compact linear formulation of the minimum weighted elementary directed cycle problem and we show that the connected components of the subgraph, that supports the optimal fractional solution of its linear relaxation without the additional vertex s / ∈ V and all the arcs (s, i), i ∈ V , are strongly connected. In section 3, we present a class of valid inequalities for the polytope P and we give a polynomial time algorithm that separates these valid inequalities in case of strongly connected planar digraphs and for β = 1. We also address an efficient lifting technique of the cycle valid inequalities of order 1 based on the incident arcs of some given vertices. In Section 4, we discuss the computational results.

A compact linear formulation for MWEDCP

The compact linear formulation of the minimum weighted elementary directed cycle problem given below is introduced in [START_REF] Maculan | Integer linear models with a polynomial number of variables and constraints for some classical combinatorial optimization problems[END_REF].

Given a planar connected digraph D = (V, A), we build an auxiliary digraph (denoted by D s ) obtained from D by adding an additional vertex s that does not belong to V and zero weighted arcs (s, k) for every vertex k belonging to V and: (i) for every vertex i belonging to V ∪ {s}, we associate a binary variable x i such that x i equals 1 if and only if the minimum weighted elementary directed cycle visits the vertex i; (ii) for every arc (i, j) belonging to A ∪ {(s, k) : k ∈ V }, we associate a binary variable y ij such that y ij equals 1 if and only if the minimum weighted elementary directed cycle passes through the arc (i, j); (iii) we define variables f k ij ≥ 0 as the flow passing through the arc (i, j), from the source vertex s to the terminal vertex k.

The mixed integer linear formulation then we obtain is as follow:

M in (i,j)∈A w ij y ij s.t. j∈Γ + (s) f k sj = x k , k ∈ V (2)
j∈Γ + (i)

f k ij - j∈Γ -(i) f k ji = 0, k ∈ V, i ∈ V \ {k} (3)
j∈Γ + (k)

f k kj - j∈Γ -(k) f k jk = -x k , k ∈ V (4)
f k ij ≤ y ij , (i, j) ∈ A, k ∈ V (5)
j∈Γ + (k)

y kj = x k , k ∈ V (6)
j∈Γ -(k)

y jk = x k , k ∈ V (7)
j∈Γ + (s)

y sj = 1, (8)
k∈V x k ≥ 2, (9)
f k ij ≥ 0, k ∈ V, (i, j) ∈ A ∪ {(s, j) : j ∈ V } (10) x i , y ij ∈ {0, 1}, i ∈ V ∪ {s}, (i, j) ∈ A ∪ {(s, j) : j ∈ V } (11)
Where, for every vertex i, Γ + (i) and Γ -(i) denote the sets of outgoing arcs from the vertex i and ingoing arcs to the vertex i, respectively. Constraints (2)-( 4) and ( 6)- [START_REF] Kaibel | On cardinality constrained cycle and path polytopes[END_REF] are flow conservation and degree constraints, respectively. The constraints (5) establish the link between the flow passing through an arc and the arc. Constraints (9) say that every elementary directed cycle may include at most 2 vertices. At optimality, the above model produces an optimal elementary directed cycle with an additional arc of type (s, j), j ∈ V . So, to find the optimal solution of the MWEDCP, we have to delete the additional arc (s, j).

Consider an optimal fractional solution (x, ŷ, f ) of the MWEDCP, let D s [ŷ] be its supporting directed sub-digraph. That is, the arc set of D s [ŷ] contains only arcs (i, j) such that ŷij > 0. On the other hand, let D[ŷ], with V as  as arc set, be the directed subgraph obtained from the digraph D s [ŷ] by deleting the vertex s and all the arcs of type (s, k) such that ŷsk > 0. Notice that the digraph D[ŷ] is not necessarily connected. But, as shown in the next lemma, each of its connected component is strongly connected.

Lemma 1 Each component of the directed sub-digraph D[ŷ] is strongly connected. Proof. Let C be a connected component of D[ŷ].
Let u and v be two vertices from C. Since the optimal fractional solution (x, ŷ, f ) satisfies the flow constraints ( 6) and ( 7) then there exists a directed path from u to v and a directed path from v to u. This completes the proof.

In the next section, we recall what we call cycle valid inequality of order β and discuss its separation problem in planar digraphs and for β = 1.

3 Separating cycle valid inequalities of order β = 1 in planar digraphs

Concept of digraph duality and cycle valid inequalities of order β

In this part, we define the concept of duality of a given planar digraph and we recall what we call cycle valid inequalities introduced in [8].

Concept of planar digraph duality

Let D be a connected planar digraph. Since the digraph D can be drawn on the plane, a face of D is defined as the area surrounded by an elementary cycle. An Unbounded face is the area surrounding the entire digraph D. Assume that the digraph D has r faces denoted g k , where k belongs to {0, . . . , r}, where g 0 is its unbounded face. The dual digraph D * associated with D is defined as follows: (i) For every face g k , k = 0, . . . , r in D, we associate a vertex v k , k = 0, . . . , r in D * ; (ii) For every arc a = (u, v) adjacent to both faces g k and g j , k = j, we associate an arc (v k , v j ) in D * . The orientation of the arc (v k , v j ) is chosen such that it crosses the original arc a = (u, v) from left to right, here left means the left side when we traverse the arc a from its tail u to its head v. As an example, in Figure 1: the vertex v 1 in the dual digraph D * corresponds to the elementary directed cycle (1, 2, 3, 6, 5, 4, 1) in D; the vertex v 2 corresponds to the elementary directed cycle [START_REF] Fortune | The directed subgraph homeomorphism problem[END_REF][START_REF] Coullard | On cycle cones and polyhedra[END_REF][START_REF] Kaibel | On cardinality constrained cycle and path polytopes[END_REF][START_REF] Ibrahim | A new class of valid inequalities for the minimum weighted elementary directed cycle problem in digraph[END_REF][START_REF] Fortune | The directed subgraph homeomorphism problem[END_REF] in D; the vertex v 3 corresponds to the elementary directed cycle (7, 8, 10, 9, 7) in D and v 0 represents the unbounded face of D. One can observe that all elementary directed cycles in D are represented by directed stars (or distars) in the dual digraph D * . We mean by a distar a set of arcs oriented in the same direction and incident to the same vertex. In Figure 1, with respect to the vertex v 2 , the arc set

{(v 2 , v 1 ), (v 2 , v 3 ), (v 2 , v 0 )} constitutes a distar.

Cycle valid inequalities of order β

Let us consider an arc set F from the planar digraph D such that:

|F ∩ A(γ)| ≤ β, ∀γ ∈ C, (12)
where A(γ) denotes the set of arcs of the elementary directed cycle γ. Such an arc set F induces valid inequalities for the polytope P.

Theorem 1 see [START_REF] Ibrahim | A new class of valid inequalities for the minimum weighted elementary directed cycle problem in digraph[END_REF] Let β be a nonnegative integer. If F is an arc set satisfying the condition [START_REF] Itai | Finding a minimum circuit in a graph[END_REF] then

a∈F y a ≤ β (13)
is a valid inequality for the polytope P.

The valid inequality [START_REF] Letchford | A fast algorithm for minimum weight odd circuits and cuts in planar graphs[END_REF], called cycle valid inequalities of order β, is introduced and studied in [START_REF] Ibrahim | A new class of valid inequalities for the minimum weighted elementary directed cycle problem in digraph[END_REF]. If the digraph D is strongly connected, then the separation problem of cycle valid inequalities of order β is NP-hard (even in the case β = 1, see [START_REF] Ibrahim | A new class of valid inequalities for the minimum weighted elementary directed cycle problem in digraph[END_REF]). Indeed, the separation problem of cycle valid inequalities of order 1 can be reduced to the existence of τ -vertex-disjoint simple directed paths in the digraph D (see [START_REF] Fortune | The directed subgraph homeomorphism problem[END_REF]), with τ ≥ 2 and τ > β. However, if the digraph D is planar, we show that the separation of the cycle valid inequalities of order 1 can be done in polynomial time.

Separation algorithm

In the sequel, we prove that separating cycle valid inequalities of order 1 in planar digraphs can be done in polynomial time.

Given an optimal fractional solution (x, ŷ, f ) of the linear relaxation ( 2)-( 11) of the MWEDCP and its corresponding sub-digraph D[ŷ], the separation problem of the cycle valid inequalities of order β consists in finding in the digraph D[ŷ] an arc subset F ⊂ Â such that:

a∈ F ŷa > β.
We recall that ŷa is the optimal (fractional) value of the arc a. As the sub-digraph D[ŷ] is possibly disconnected, let C i [ŷ], i = 1, . . . , c be its different connected components and consider the dual components

C * i [ŷ], i = 1, . . . , c drawn from the different components C i [ŷ], i = 1, . . . , c. Proposition 1 For every vertex pair (u, v) of the dual component C * i [ŷ], i = 1, . . . , c, a u-v elementary directed path in C * i [ŷ], i = 1, . . . , c, defines a cycle valid inequality of order 1 w.r.t the component C i [ŷ], i = 1, . . . , c.
Proof. Let u and v be two vertices belonging to the dual component

C * i [ŷ]. Let P u,v be an elementary directed path from u to v in C * i [ŷ].
Let Fi be the arc set of the component C i [ŷ] containing all arcs that are crossed by an arc of A(P u,v ), the arc set of the u-v directed path P u,v . Since every arc of P u,v crosses an arc from C i [ŷ] that belongs to two distinct elementary directed cycles of C i [ŷ], then we have

| Fi ∩ A (γ) | ≤ 1, ∀ γ ( 14
)
where γ is an elementary directed cycle of the component In Figure 1, in the dual digraph D * , the elementary directed paths (v 2 , v 0 , v 1 ) and (v 2 , v 0 , v 3 ) correspond to the cycle valid inequalities of order 1 y 3,6 + y 8,6 ≤ 1 and y 8,6 + y 8,10 ≤ 1 in D, respectively. One can easily verify in D that the arcs [START_REF] Bauer | A Polyhedral Approach to the Weighted Girth Problem[END_REF][START_REF] Fortune | The directed subgraph homeomorphism problem[END_REF] and [START_REF] Ibrahim | A new class of valid inequalities for the minimum weighted elementary directed cycle problem in digraph[END_REF][START_REF] Fortune | The directed subgraph homeomorphism problem[END_REF] belong to distincts elementary directed cycles in D.

C i [ŷ]. Therefore, a∈ Fi y a ≤ 1, ∀ i,
Generally, since the valid inequalities derived on a given sub-digraph, (as D[ŷ] the sub-digraph that supports the optimal (fractional) solution (x, ŷ, f ) of the linear relaxation of the model ( 2)-( 11)), are not valid considering the entire digraph D = (V, A), we may resort to a lifting technique to extend such valid inequalities for the whole polytope P. This fact is illustrated in what follows Let the following planar digraph D = (V, A) be an instance of the MWEDCP. Consider C = {γ 1 , γ 2 , γ 3 , γ 4 , γ 5 , γ 6 , γ 7 } the set of all elementary directed cycles in D. We represent an elementary directed cycle by an ordered sequence of its vertices. Thus, γ 1 = (5, 10, 11, 6, 5), γ 2 = (5, 10, 11, 12, 7, 6, 5), γ 3 = (5, 10, 11, 12, 13, 8, 7, 6, 5), γ 4 = (8, 9, 14, 13, 8), γ 5 = ( [START_REF] Ibrahim | Strong flow-based formulation for the shortest path problem in digraphs with negative cycles[END_REF][START_REF] Ibrahim | Valid inequalities and lifting procedures for the shortest path problem in digraphs with negative cycles[END_REF], [START_REF] Ibrahim | Valid inequalities and lifting procedures for the shortest path problem in digraphs with negative cycles[END_REF][START_REF] Fortune | The directed subgraph homeomorphism problem[END_REF], (6, 5)}, {(7, 6), [START_REF] Ibrahim | Valid inequalities and lifting procedures for the shortest path problem in digraphs with negative cycles[END_REF][START_REF] Fortune | The directed subgraph homeomorphism problem[END_REF], [START_REF] Ibrahim | Valid inequalities and lifting procedures for the shortest path problem in digraphs with negative cycles[END_REF][START_REF] Itai | Finding a minimum circuit in a graph[END_REF], (12, 7)}, {(12, 13), [START_REF] Letchford | A fast algorithm for minimum weight odd circuits and cuts in planar graphs[END_REF][START_REF] Ibrahim | A new class of valid inequalities for the minimum weighted elementary directed cycle problem in digraph[END_REF] In the following lines, we address a lifting technique that can be polynomially implemented and enables to transform the valid inequality (15) into the following inequality y 11,6y 15,11 + y 12,7 + y 12,13 + y 9,14y 4,9 ≤ 1 [START_REF] Minoux | An Hierarchy of Strong Block-Decomposable Linear Relaxations for 0-1 MIPs[END_REF] that is satisfied by all incidence vectors of the elementary directed cycles of D. Moreover, the lifted valid inequality ( 16) is violated by the fractional solution supported by the sub-digraph represented on the left side of Figure 3. Thus, ŷ11,6 -ŷ15,11 + ŷ12,7 + ŷ12,13 + ŷ9,14 -ŷ4,9 = 0.5 -0 + 0.25 + 0.25 + 0.75 -0 = 1.75 > 1.

Given an instance D = (V, A) the MWEDCP and its corresponding sub-digraph D[ŷ]. Consider the arc subset F of A such that F verifies the condition [START_REF] Lingas | Efficient approximation algorithms for shortest cycles in undirected graphs[END_REF]. For any arc a = (u, v) ∈ F , we denote by Fthe set of all arcs in A \ Â such that these arcs have the node u as tail. Similarly, we denote by F + the set of all arcs in A \ Â such that these arcs have v as head. That is is valid for the polytope P.

F -= {(x, u) ∈ (A \ Â) : (u, v) ∈ F } and F + = {(v, x) ∈ (A \ Â) : (u, v) ∈ F }.
We call the valid inequality ( 17) the F -lifted cycle valid inequality of order 1. Such a valid inequality is close in spirit to the so-called simple lifted valid inequalities introduced in [START_REF] Ibrahim | Valid inequalities and lifting procedures for the shortest path problem in digraphs with negative cycles[END_REF] for the s-t elementary directed path problem. However, here contrary to the simple lifted valid inequalities, we do not consider all the arcs of A c = A \ Â, the arcs of the complementary digraph of D[ŷ]. We recall that  is the arc set of the sub-digraph

D[ŷ].
Proof. According to the condition ( 14), every elementary directed cycle in the subdigraph D[ŷ] intersects with at most one arc of the arc set F . Consider an elementary directed cycle γ ∈ C, we can distinguish cases:

1. if γ is an elementary directed cycle of the sub-digraph D[ŷ], it is that the inequality ( 17) is valid for the polytope P.

2. if γ is an elementary directed cycle such that A(γ) ∩ Â = ∅, the inequality (15) remains valid. Thus, as F = ∅, the sum of terms on the left hand side of the inequality are non-positive. 17) is satisfied by the incidence vector of such an elementary directed cycle γ as generally every directed cycle γ of C which passes through an arc a = (u, v) ∈ F has already passed through an arc ingoing to the head of the arc a and/or must pass through an arc outgoing from the tail of a.

if γ is an elementary directed cycle such that

This completes the proof.

Notice that from a given cycle valid inequality of order 1 generated on the sub-digraph D[ŷ], the computational time required to construct its corresponding F -lifted cycle valid inequality of order

1 is O(|V |) in grid digraphs and is O(|V | 2 ) in complete digraphs.
The steps of our cutting plane procedure are summarized in algorithm 1 where: the procedure solveLP returns an optimal solution of the linear optimization problem over the polyhedron Q given as its argument; the procedure isInteger returns true if its argument is a binary vector and false otherwise; the procedure SupportDigraph generates the supporting digraph of the fractional solution given as its argument; the procedure ConnectedComponents returns the connected components of the graph given as its argument; the procedure DualDigraph builds and returns the dual digraph of the digraph given as its argument; the procedure stViolatedPath returns true if there is a violated directed path in the dual sub-digraph given as its first argument (the path is returned as its second argument) and false otherwise. We say that a directed path is violated if its corresponding cycle valid inequalty is violated by the optimal fractional solution (x, ŷ); The procedure addLiftedIneq appends to the description of the polyhedron Q the lifted form of the cycle valid inequality of order 1 derived from the violated directed path π given as its second argument. Finally, the procedure Weight gives the value of the optimal solution after the cutting plane phase.

Computational results

In this section, we present computational results showing the efficiency of the cutting plane algorithm devised for the MWEDCP.

We deal with planar grid digraphs that contain negative weighted elementary directed cycles. We recall that the realistic instances from some libraries like TSPLIB are not planar in general. and does not contain negative weighted elementary cycles, as the of its arcs are positive. instances are randomly generated directed grids D of different number of vertices: 50 (with 85 arcs), 100 (with 180 arcs) and 200 (with 370 arcs). As the neighboor set of every vertex in a grid D is known in advance, we have just to choose randomly the orientation of every arc (i, j). For this, we generate uniformly two numbers x and y, if x < y then the arc is oriented from i to j otherwise the arc is oriented from j to i. In regard to the arc weights, with each arc (i, j), we associate a weight z 1 -10 such that 0 ≤ z 1 ≤ 20. In this way, the weights have a probability equals to 1 2 to be whether negative or not. So, all our generated instances contain negative weight elementary directed cycles (see the values of the third column of Table 2 are all different to 0). Therefore, our instances cannot be solved using the algorithms proposed by Roditty and Williams ([17]). Notice that among all randomly generated instances, we select only instances having a significant integrality gap (see the fourth and fifth columns of Table 2). Our algorithm is implemented in C and all computations have been carried out on a computer equipped with a 1.50 GHz Intel (R) core (TM) 2 CPU. All our instances are solved using the open-source software glpk. The computational results are presented in Table 2 (the entries of Table 2 are given in Table 1). We can make the following comments. Consider the first line of Table 2 represented by the instance named 50-1, such an instance has 50 vertices and 85 arcs. After, the resolution of the linear relaxation of model ( 2)-( 11) applied to 50-1, as shown by the second column of Table 2, the sub-digraph D s [ŷ] that supports the optimal (fractional) solution has 29 vertices, 38 arcs and contains 21 negative weighted directed cycles (see column 3). Assuming that the sub-digraph D[ŷ], obtained from the subdigraph D s [ŷ] by deleting the source vertex s and all arcs {(s, j) : ŷsj > 0}, is connected and according to algorithm 1 described above, we build its associated dual sub-digraph D * [ŷ]. After, we generated all directed paths in the dual sub-digraph corresponding to cycle valid inequalities of order 1. Among all generated cycle valid inequalities, as shown on the 8 th column we only keep 43 are violated by the optimal (fractional) solution (x, ŷ, f ). As these valid inequalities are derived w.r.t the sub-digraph D[ŷ], we use the lifting technique described in the previous section to transform the generated cycle valid inequalities into valid inequalities for the whole polytope P. The lifted valid inequalities are called the F -lifted cycle valid inequalities of order 1. Therefore, in a cutting plane process, we resort to the F -lifted cycle valid inequalities of order 1 to strengthen the linear relaxation of the MWEDCP model ( 2)- [START_REF] Ibrahim | Valid inequalities and lifting procedures for the shortest path problem in digraphs with negative cycles[END_REF]. Considering all tested instances, strengthening with corresponding F -lifted cycle valid inequalities of order 1 produces in average about 26.47% gap improvement showing the efficiency of these lifted inequalities. However, for some of our generated instances, we do not obtain any violated lifted cycle valid inequality. Thus, for such instances (50-7, 100-10, 200-1 and 200-7), we have #cvi = 0 as shown on the 8 th column of Table 2. This explains in part the fact that these instances do not present any gap improvement (see on corresponding 9 th column, gap = 0). Nevertheless, we observe that for the instances 50-2 and 100-3 the gap is totally closed (thus, gap = 100%).

Table 1: Meaning of the entries of Table 2.

Column Description

Instances Name of the instance. The name follows the format x-y where x indicates the number of vertices and y the instance number.

(|V S |, |A S |) Number of vertices and number of arcs in the supporting digraph of the first optimal fractional solution.

#NDC Number of negative weight elementary directed cycles. The time limit is set to 24 hours. The symbol na indicates that the value is not available.

#DC Number of elementary directed cycles. The computational time limit is set to 24 hours. The symbol na indicates that the value is not available.

z of a optimal fractional solution of the linear relaxation (2)-( 11)

z * The weight of an optimal integer elementary directed cycle zcvi Value of an optimal solution obtained after strengthening with all F -lifted cycle valid inequalities of order 1.

#cvi Number of violated F -lifted cycle valid inequalities of order β = 1 generated gap(%) The gap closed (in percent), that is gap = zcvi -z z * -z × 100.

Conclusion

We address an efficient cutting plane algorithm for the minimum weighted elementary directed cycle problem (MWEDCP) in planar digraphs containing negative weight elementary directed cycles. The algorithm is based on the linear relaxation of a mixed integer linear formulation of the MWEDCP and the use of a lifting technique applied to a class of valid inequalities called cycle valid inequalities that are parametrized by an integer called the order (the inequality right-hand side). We prove that separating cycle inequalities of order 1 in planar digraphs can be done in polynomial time. In addition, as the inequalities are derived over the sub-digraph that supports the optimal (fractional) solution, with respect to the entire digraph, we deal with a lifted form of cycle valid inequalities of order 1 named the F -lifted cycle valid inequalities of order 1. The lifting technique performed is also efficient, as it is based on some of the incident arcs of every vertex in the considered sub-digraph.

The computational results carried out randomly generated planar digraphs (grids) of size set between 50 and 200 feature that a strengthening with F -lifted cycle valid inequalities of order 1 can significantly improve the integrality gap. Indeed, in average, for all 29 instances, we obtain 26.47% gap improvement. Moreover, there exists instances for which the strengthening process directly displays the optimal integer elementary directed cycle.

Figure 1 :

 1 Figure 1: A digraph D (left) and its dual D * (right).

Remark 1 y a ≤ 1 . 2

 112 is a valid inequality w.r.t the component C i [ŷ], this completes the proof. The above inequality is also valid for all elementary directed cycles in D[ŷ] = c i=1 C i [ŷ]. Moreover, as the arcs of each Fi are in a different component, namely C i [ŷ], we also have the following valid inequality a∈( c i=1 Fi ) Theorem The separation of cycle valid inequalities of order 1 induced by elementary directed paths of the dual D * [ŷ] of a planar sub-digraph D[ŷ] is solvable in polynomial time. Proof. W.l.o.g. we assume that the sub-digraph D[ŷ] is connected. By proposition 1, any u-v elementary directed path in the dual digraph D * [ŷ] can be associated with a cycle valid inequality of order 1 in D[ŷ]. Using a depth first search (DFS) algorithm from a vertex u ∈ D * [ŷ] in view to reach a vertex v ∈ D * [ŷ], we generate in polynomial time a u-v elementary directed path that induces a corresponding cycle valid inequality of order 1, as the complexity time of a DFS algorithm is O(|V | + |A|). To really solve the separation problem, we consider only the u-v elementary directed paths of the dual digraph such that the sum of the values of the of the arcs of the corresponding cycle valid inequality is up to 1.

Figure 2 :

 2 Figure 2: An instance D of the MWDECP.

Figure 3 :

 3 Figure 3: A sub-digraph of D (left) and its dual (right).

Proposition 2 F y a -

 2 Let F be an arc subset from  verifying the condition (14) w.r.t. the sub-digraph D[ŷ]. Then the following lifted inequality a∈

  A(γ) ∩ Â = ∅ and A(γ) ∩ A c = ∅, it means that the elementary directed cycle γ has some of its arcs in D[ŷ] and the others in the complementary subgraph of D[ŷ]. The inequality (

Algorithm 1 : 2 Separate 3 Q ← P 4 while 7 NotIntegerSolution ← False 8 else 9 D 10 L

 123478910 Strengthening of the linear relaxationData: P // The linear relaxation of (2)-(11)Result: LB // Bestlower bound 1 begin ← True; NotIntegerSolution ← True Separate and NotIntegerSolution do 5(x, ŷ, f ) ← solveLP(Q) 6 if isInteger((x, ŷ)) then ← SupportDigraph((x, ŷ))

  According to the proposition 1, from the v 1v 4 directed pathp 1 = (v 1 , v 2 , v 3 , v 0 , v 4 ) of D * [ŷ]represented by an ordered sequence of its vertices, we deduce that the inequalityy11,6 + y 12,7 + y 12,13 + y 9,14 ≤ 1 (15) is valid over the sub-digraph D[ŷ]. Indeed, it is satisfied by the incidence vectors of all elementary directed paths γ 1 , γ 2 , γ 3 and γ 4 of the sub-digraph D[ŷ]. However, it is violated by the incidence vectors of the elementary directed cycles γ 5 , γ 6 and γ 7 of the digraph D. Inequality (15) is then not valid w.r.t the digraph D.

	, (8, 7),
	(12, 7)} and {(8, 9), (9, 14), (14, 13), (13, 8)}, respectively. With v 0 associated to the un-
	bounded face. The elementary directed cycles of D[ŷ] are γ 1 = (5, 10, 11, 6, 5), γ 2 =
	(5, 10, 11, 12, 7, 6, 5), γ 3 = (5, 10, 11, 12, 13, 8, 7, 6, 5) and γ 4 = (8, 9, 14, 13, 8).

Table 2 :

 2 Efficiency of F -lifted cycle valid inequalities in strengthening of the linear relaxation of the MWEDCP model (2)-[START_REF] Ibrahim | Valid inequalities and lifting procedures for the shortest path problem in digraphs with negative cycles[END_REF].

	Instances (|V S |, |A S |) #NDC #DC z	z *	zcvi	#cvi gap(%)
	50-1	(29,38)	21	41	-55.5	-38	-46.54	43	51.20
	50-2	(8,10)	8	11	-27	-22	-22	8	100
	50-3	(32,37)	33	44	-72.5	-41	-70		7.94
	50-4	(17,21)	9	13	-61	-36	-59	5	8
	50-5	(11,14)			-25.5	-16	-23.5	3	21.05
		(26,33)	26	27	-69.66 -44	-59.5	27	39.59
	50-7	(22,24)	70	90	-75.5	-68	-75.5	0	0
	50-8	(16,20)	12	17	-44.5	-22	-36.5	5	35.56
	50-9		14	27	-51.5	-40	-51.5	3	0
	50-10	(32,38)	34	42	-66.75 -51	-60.5	27	39.68
	100-1	(33,43)	24	44	-52.25 -26	-50.5	19	6.67
	100-2	(30,34)	66	100 -63.66 -57	-61.5	9	32.43
	100-3	(17,20)	16	24	-49	-44	-44	11	100
	100-4	(26,29)	39	44	-56	-43	-50.25	8	44.23
	100-5	(28,31)	212	309 -65	-50	-54.75	8	68.33
	100-6	(35,40)	34	53	-56.5	-43	-54.5	6	14.81
	100-7	(43,46)	408	902 -77	-61	-74.25	2	17.19
	100-8	(25,31)	29	41	-57.25 -51	-56.75	10	8
	100-9	(39,43)	57	70	-85.5	-70	-85	5	3.23
	100-10	(21,24)	70	87	-65	-57	-65	0	0
	200-1	(30,34)	67	84	-104.5 -68	-104.5	0	0
	200-2	(15,18)	68	106 -53.5	-46	-51.16	3	31.20
	200-3	(44,50)	na	na	-137.5 -106 126.25 14	35.71
	200-4	(30,34)	na	na	-48.5	-39	-44	10	47.37
	200-5	(39,47)	78	102 -73	-55	-71.62	2	7.67
	200-6	(33,40)	na	na		-33	-45.66	15	6.22
	200-7	(20,22)	na	na	-58	-37	-58	0	0
	200-8	(35,45)	na	na	-72.5	-62	-70	46	23.81
	200-9	(54,60)	na	na		-98	-110.75 9	17.74

Acknowledgments We gratefully acknowledge the referees for their careful reading and constructive comments.


Bruno Galerne

Arthur Leclaire

Lionel Moisan

MICROTEXTURE INPAINTING THROUGH GAUSSIAN CONDITIONAL SIMULATION

Keywords: Inpainting, Gaussian texture, texture synthesis, conditional simulation, kriging

Image inpainting consists in filling missing regions of an image by inferring from the surrounding content. In the case of texture images, inpainting can be formulated in terms of conditional simulation of a stochastic texture model. Many texture synthesis methods thus have been adapted to texture inpainting, but these methods do not offer theoretical guarantees since the conditional sampling is in general only approximate. Here we show that in the case of Gaussian textures, inpainting can be addressed with perfect conditional simulation relying on kriging estimation. We thus obtain a microtexture inpainting algorithm that is able to fill holes of any shape and size in an efficient manner while respecting exactly a stochastic model.

INTRODUCTION

Since the seminal papers [START_REF] Lewis | Texture Synthesis for Digital Painting[END_REF][START_REF] Lewis | Methods for Stochastic Spectral Synthesis[END_REF][START_REF] Van Wijk | Spot noise texture synthesis for data visualization[END_REF], random phase fields have proven to be a powerful tool for microtexture modeling. In particular, given a microtexture, a Gaussian model can be easily derived and sampled [START_REF] Galerne | Random Phase Textures: Theory and Synthesis[END_REF] and latter approximated by a spot noise model to allow for fast by-example synthesis [START_REF] Galerne | A Texton for Fast and Flexible Gaussian Texture Synthesis[END_REF]. In addition, the theoretical benefits of the Gaussian model have been exploited to address texture analysis [START_REF] Desolneux | A compact representation of random phase and Gaussian textures[END_REF] and texture mixing [START_REF] Xia | Synthesizing and Mixing Stationary Gaussian Texture Models[END_REF]. In this paper, we take advantage of the Gaussian model in order to address microtexture inpainting.

The inpainting problem consists in recovering hidden parts of an image based on the surrounding content. Many methods have been proposed for inpainting, falling into three categories. The first are deterministic methods which use a variational or PDE-based framework to fill the missing regions by connecting the edges or level lines in a way that respects the Gestaltist's good continuation principle [START_REF] Masnou | Level lines based disocclusion[END_REF][START_REF] Bertalmio | Image Inpainting[END_REF][START_REF] Ballester | Filling-in by joint interpolation of vector fields and graylevels[END_REF][START_REF] Chan | Nontexture inpainting by curvaturedriven diffusions[END_REF][START_REF] Chan | Mathematical models for local nontexture inpaintings[END_REF][START_REF] Esedoglu | Digital inpainting based on the Mumford-Shah-Euler image model[END_REF][START_REF] Tschumperlé | Fast anisotropic smoothing of multi-valued images using curvature-preserving PDE's[END_REF][START_REF] Bornemann | Fast image inpainting based on coherence transport[END_REF]. The common drawback of these methods is their inability to handle texture precisely.

The second are stochastic methods. Inpainting is by essence an ill-posed problem. However, if one has a pertinent stochastic model p(u) for the whole image u, a perfect inpainting solution can be obtained by conditional simulation, that is simulating an image u with distribution p(u) given the known values of u. This stochastic point of view is not really helpful in general since one cannot sample the distribution of realistic images. However, this approach is of particular interest when dealing with texture images since they are often well represented by tractable stochastic models. The first works in the direction of stochastic inpainting suggest to fill the masked zone by matching specific statistics, like co-occurrence matrices [START_REF] Malzbender | A Context Sensitive Texture Nib[END_REF] or histograms of filter responses [START_REF] Igehy | Image replacement through texture synthesis[END_REF]. Truly conditional methods later appeared based on Markov random fields, which are inherently adapted to conditional sampling [START_REF] Winkler | Image Analysis, Random Fields and Markov Chain Monte Carlo Methods: A Mathematical Introduction[END_REF]. Using patch comparison to approximate local conditional sampling, the authors of [START_REF] Efros | Texture synthesis by nonparametric sampling[END_REF][START_REF] Wei | Fast texture synthesis using treestructured vector quantization[END_REF] were led to a celebrated non-parametric texture synthesis algorithm which can also be used for hole filling. However, these synthesis algorithms may degenerate after a spatial boundary (growing garbage effect), and besides are quite dependent on the filling order, thus not ensuring the prolongation of large-scale features.

Taking profit of the first attempts in deterministic and stochastic inpainting, hybrid methods have emerged, which consist in filling both geometric and textural contents in a coherent manner. After the proposition of [START_REF] Bornard | Missing data correction in still images and image sequences[END_REF] to extend non-parametric sampling to general image synthesis, several authors improved this work by varying the pixel filling order [START_REF] Drori | Fragment-based image completion[END_REF][START_REF] Criminisi | Region filling and object removal by exemplar-based image inpainting[END_REF][START_REF] Cao | Geometrically guided exemplar-based inpainting[END_REF]. Other methods suggest to separate geometry and texture and to inpaint them either with two different techniques [START_REF] Bertalmio | Simultaneous structure and texture image inpainting[END_REF] or to minimize a unified functional [START_REF] Elad | Simultaneous cartoon and texture image inpainting using morphological component analysis (MCA)[END_REF]. The last work paves the way for inpainting with adapted dictionaries [START_REF] Guleryuz | Nonlinear approximation based image recovery using adaptive sparse reconstructions and iterated denoising-part II: adaptive algorithms[END_REF][START_REF] Mairal | Sparse Representation for Color Image Restoration[END_REF][START_REF] Cai | A framelet-based image inpainting algorithm[END_REF][START_REF] Peyré | Texture Synthesis with Grouplets[END_REF]. More recently, several authors gave a unified view of example-based image inpainting in a variational framework [START_REF] Aujol | Exemplar-based inpainting from a variational point of view[END_REF][START_REF] Arias | A variational framework for exemplar-based image inpainting[END_REF].

These methods are able to produce a satisfying result on a large class of images with a reasonable computational time, but in general do not provide any theoretical guarantee in the synthesized content. In contrast, we demonstrate in this paper that if one restricts to a Gaussian random model, microtexture inpainting can be achieved by perfect conditional sampling. Our algorithm relies on kriging estimation [START_REF] Lantuéjoul | Geostatistical Simulation: Models and Algorithms[END_REF] and fills holes with two independent components, namely the kriging component which extends long scale structures from the hole boundary, and the innovation component which adds an independently generated textural content containing small scale details. The proposed algorithm can inpaint holes of any shape and size, provided that a faithful Gaussian texture model can be estimated within the unmasked image area. For large holes, this algorithm does not suffer from the growing garbage effect encountered in [START_REF] Efros | Texture synthesis by nonparametric sampling[END_REF], and besides, edge-like structures in the textures are naturally restored in a way that respects the covariance.

Texture inpainting was addressed in [START_REF] Chandra | Texture Interpolation Using Ordinary Kriging[END_REF][START_REF] Jassim | Image Inpainting by Kriging Interpolation Technique[END_REF] with kriging estimation, but this does not amount to sampling a global Gaussian model since the innovation component is missing (see Section 4). Also, texture synthesis was addressed in [START_REF] Raad | Conditional Gaussian Models for Texture Synthesis[END_REF] by progressive filling with Gaussian conditional patch models, but in this case too the global model is not Gaussian.

GAUSSIAN CONDITIONAL SIMULATION

In this section, following [START_REF] Lantuéjoul | Geostatistical Simulation: Models and Algorithms[END_REF], we rely on kriging estimation to explain the perfect conditional sampling of Gaussian random vectors that is at the heart of our microtexture inpainting algorithm.

Let Ω be a finite set. For a set A ⊂ Ω and a function f : Ω → R we denote by |A| the cardinality of the finite set A, and f |A the restriction to A of the function f . Let (F (x)) x∈Ω be a Gaussian vector with mean zero (for the sake of simplicity) and covariance function

Γ(x, y) = Cov(F (x), F (y)) = E(F (x)F (y)), x, y ∈ Ω.
Given prescribed values ϕ : C → R on a set C ⊂ Ω of conditioning points, conditional Gaussian simulation consists in sampling the conditional distribution of F given that F |C = ϕ.

One defines the kriging estimator at x ∈ Ω as the conditional expectation

F * (x) = E( F (x) | F (c) , c ∈ C ).
This means that F * (x) is the best least-square estimation of F (x) that can be obtained as a measurable function of (F (c)) c∈C .

In particular for all x ∈ C, F * (x) = F (x). A standard result of probability theory [START_REF] Doob | Stochastic processes[END_REF] ensures that in the Gaussian case F * (x) is the orthogonal projection of F (x) on the subspace of linear combinations of (F (c)) c∈C (for the L 2 -distance between square-integrable random variables). Hence, there exist coefficients (λ c (x)) c∈C such that

F * (x) = c∈C λ c (x)F (c). (1)
These deterministic numbers (λ c (x)) c∈C are called the kriging coefficients. The following theorem is the core result for conditional Gaussian simulation. Proof. Due to the orthogonal projection, (F * , F -F * ) is a Gaussian vector whose components F * and F -F * are uncorrelated and thus independent.

A conditional sample of F given F |C = ϕ can thus be obtained with ϕ * + G -G * where the kriging component ϕ * is defined by ϕ * (x) = c∈C λ c (x)ϕ(c), x ∈ Ω and where the independent component G -G * is called the innovation component. To do so one needs to compute the kriging coefficients (λ c (x)) c∈C introduced in (1). Since F * (x) -F (x) is orthogonal to F (c), we obtain that E(F * (x)F (c)) = E(F (x)F (c)) = Γ(x, c), and using (1) to substitute F * (x) shows that (λ c (x)) c∈C is a solution of the following |C| × |C| linear system

∀c ∈ C, d∈C λ d (x)Γ(d, c) = Γ(x, c). (2)
In the following, we will assume that the matrix Γ |C×C is invertible (which turned out to be true in all our inpainting experiments) so that this system admits a unique solution.

If we denote by Λ = (λ c (x)) x∈Ω, c∈C the |Ω| × |C| matrix of all kriging coefficients, the previous system rewrites Λ = Γ |Ω×C Γ -1 |C×C , where Γ |Ω×C and Γ |C×C are the restrictions of the covariance matrix Γ. Finally, writing the Gaussian vectors in column, sampling F given F |C = ϕ amounts to sampling

Λϕ + G -ΛG |C , where Λ = Γ |Ω×C Γ -1 |C×C (3)
where G has the same distribution as F and is independent of F . The required number of operations is O(|C| 3 + |Ω||C|) in addition to the simulation cost for G.

MICROTEXTURE INPAINTING

In this section, we use the Gaussian conditional sampling scheme described in the previous section to perform microtexture inpainting. The input is a masked texture u : Ω → R the values of which are unknown on the mask domain M ⊂ Ω, and we want to sample these values based on the unmasked values in M c . This amounts to adopting a Gaussian texture model U on Ω, and to sample it on M given the surrounding values. This raises two main difficulties: 1) One must estimate the Gaussian model for U from the masked image u; 2) The computation of the kriging components must be practical, without a cubic complexity in |M c | 3 .

Gaussian Texture Model Estimation

The Gaussian model that we use for conditional simulation is given by an asymptotic discrete spot noise (ADSN), which is the convolution of a finitely-supported function h : Z 2 → R with a normalized Gaussian white noise on the discrete plane Z 2 . This ADSN has zero mean and covariance Γ(x, y) = (h * h)(x -y), x, y ∈ Z 2 , where h(x) = h(-x). Note that Γ(x, y) only depends on x -y since the ADSN is stationary. In order to estimate an ADSN from the masked exemplar, we consider a subimage v of u defined on a subdomain ω ⊂ M c . Following [START_REF] Galerne | Random Phase Textures: Theory and Synthesis[END_REF][START_REF] Galerne | A Texton for Fast and Flexible Gaussian Texture Synthesis[END_REF], we compute the empirical mean v and the normalized spot As one can see on the result, the conditional inpainting is able to fill both small and large holes, whatever the regularity of the boundary.

t v = 1 √ |ω| (v -v) that we
extend by zero-padding to Z 2 . The ADSN is then simply given by the convolution t v * W where W is a normalized Gaussian white noise on Z 2 . We must choose ω large enough to capture a representative piece of the desired texture.

Textural Inpainting Algorithm

Ideally one would consider the conditioning set C = M c , but because of the O(|C| 3 ) complexity, we restrict C to a band of width w pixels along the outside boundary of M (Let us mention that in the case where U has a Markov property with a neighborhood of width w, this restriction of the conditioning set would still give a perfect conditional sampling given all the values in M c ). We then sample U on M given U |C = u |C with the algorithm explained in Section 2. However, we exploit the stationarity of our Gaussian model to compute efficiently the kriging component. Indeed, thanks to stationarity, a matrix-vector multiplication with the covariance matrix Γ is a convolution by the image c v = t v * tv . Hence to compute the kriging component

ϕ * = Λϕ = Γ |Ω×C (Γ -1 |C×C ϕ) we first compute ψ = Γ -1
|C×C ϕ that we extend by zero-padding to get a function Ψ on Ω, and then compute ϕ * = c v * Ψ. The whole process is summarized in Algorithm 1.

Algorithm 1: Microtexture inpainting

Inputs: Mask M ⊂ Ω, values of the texture u on M c -From a restriction v of u to a subdomain ω ⊂ M c , compute the mean v, t v = 1 √ |ω| (v -v) ∈ R Ω
(extended by 0 outside ω), and

c v = t v * tv -Store the matrix Γ |C×C (c, d) = c v (c -d), c, d ∈ C -Draw an ADSN sample U = t v * W -Compute ψ 1 = Γ -1 |C×C (u |C -v) and ψ 2 = Γ -1 |C×C U |C -Extend ψ 1 and ψ 2 by zero-padding to get Ψ 1 and Ψ 2 -Compute the kriging components (u -v) * = c v * Ψ 1 and U * = c v * Ψ 2 Output: v + (u -v) * + U -U * Concerning the complexity, Γ -1
|C×C ϕ is computed with standard numerical techniques in O(|C| 3 ). Then, the kriging components is obtained with the discrete Fourier transform in O(|Ω| log |Ω|), which is the complexity to generate the ADSN sample U . Thus, the number of conditioning points must stay low for the conditional simulation to be feasible.

Extension to Color

Images Algorithm 1 can be straightforwardly adapted to color textures. Indeed, the simulation of the color ADSN model is identical to the grayscale case, provided that we convolve the three channels by the same Gaussian white noise [START_REF] Galerne | Random Phase Textures: Theory and Synthesis[END_REF]. Besides, the kriging framework can be adapted by taking the conditional expectation with respect to the three channels of conditional values. In the end, it amounts to computing a larger matrix Γ containing all the cross-correlations between the three channels. Then, the matrix Γ -1

|C×C is of size 3|C| × 3|C|, and the convolution involved in Γ Ω×C is now a convolution by a matrix-valued function.

RESULTS

First, in Fig. 1 we validate the methodology based on conditional simulation by inpainting a masked texture with a previously estimated Gaussian model. One can see that, in contrast to PDE-based methods, this algorithm is able to fill holes in textures of any size and shape. In particular, the kriging component extends edge-like structures and the innovation component adds local details, while globally preserving the texture covariance. This reveals the crucial difference with methods solely using the kriging component [START_REF] Chandra | Texture Interpolation Using Ordinary Kriging[END_REF][START_REF] Jassim | Image Inpainting by Kriging Interpolation Technique[END_REF].

Let us also discuss on this example the influence of the parameter w which is the width of the conditioning points. Intuitively, this value must not be too small in order to get a sufficient summary of the surrounding context based on the conditioning set C, but also not too large to keep a reasonable computational time (see Subsection 3.2). In our experiments we chose the value w = 3. We observed in other experiments (not shown here) that this is the minimal value ensuring a good line preservation. As already said, with a Gaussian Markov model, restricting the conditioning points to the mask border would not affect the conditional distributions. Let us mention however that the ADSN models that we use here do not satisfy the Markov property. But the visual impact of this conditioning restriction seems visually harmless. Now, let us present microtexture inpainting with a Gaussian model estimated on the masked exemplar. In Fig. 2, we inpaint a circular hole in the left part of the texture while learning the Gaussian model on the right part. One can observe that the algorithm produces a satisfying result on the two first examples because the hole cannot be perceived in the inpainted texture and it may be difficult to distinguish the original and the output. The result is however slightly less satisfying on the third example because the corresponding Gaussian model exhibits interference patterns which are not visible in the original.

In Fig. 3, we compare our method to the one of [START_REF] Chan | Mathematical models for local nontexture inpaintings[END_REF] based on the total variation (TV), and to the example-based method of [START_REF] Criminisi | Region filling and object removal by exemplar-based image inpainting[END_REF]. As expected, textures are completely lost with the TV-based model. In contrast, the result of [START_REF] Criminisi | Region filling and object removal by exemplar-based image inpainting[END_REF] is more precise, but the border of the inpainted domain is still visible. It is also the case with our method, which produces a more textured pattern.

Our Gaussian inpainting algorithm is also compared to Efros-Leung algorithm in Fig. 4. As one can see, the frequency content of the estimated Gaussian model is precisely reproduced in the filled zone. In contrast, Efros-Leung algorithm fails to reproduce the long-range correlations which exist in this microtexture. In Fig. 4 we also show the kriging component v+(u-v) * involved in the conditional simulation; one can see that this component contains long-range diagonal correlations that are characteristic of this wood texture.

To conclude, this paper shows that Gaussian conditional sampling is an interesting algorithm to fill large holes in microtextures provided that the Gaussian model has been properly estimated. For small to medium holes, the proposed algorithm is very efficient (the results of Figures 2,3 were obtained in about 1 sec.). Its main limitation is that the Gaussian model is only valid for a limited class of texture. Still, extending Gaussian conditional method to non-stationary Gaussian random fields could be of potential interest to inpaint nontexture images.

Original

TV inpainting [START_REF] Chan | Mathematical models for local nontexture inpaintings[END_REF] Our result Criminisi et al. [START_REF] Criminisi | Region filling and object removal by exemplar-based image inpainting[END_REF] Fig. 3. Comparison with [START_REF] Chan | Mathematical models for local nontexture inpaintings[END_REF][START_REF] Criminisi | Region filling and object removal by exemplar-based image inpainting[END_REF]. On this example taken from [START_REF] Criminisi | Region filling and object removal by exemplar-based image inpainting[END_REF], we compare our result (bottom left) to [START_REF] Chan | Mathematical models for local nontexture inpaintings[END_REF] (top right, obtained with the implementation available at [START_REF] Getreuer | Total Variation Inpainting using Split Bregman[END_REF]) and to [START_REF] Criminisi | Region filling and object removal by exemplar-based image inpainting[END_REF] (bottom right). The Gaussian model that we used was estimated in the red rectangle shown on the original (top left). As one can see, the TV inpainting is not able to preserve texture. In contrast, the result of [START_REF] Criminisi | Region filling and object removal by exemplar-based image inpainting[END_REF] is comparable to ours.

Original Kriging component

Efros-Leung Our result

Theorem 2 . 1 (

 21 [START_REF] Lantuéjoul | Geostatistical Simulation: Models and Algorithms[END_REF]). F * and F -F * are independent. Consequently, if G is independent of F and has the same distribution, then H = F * + (G -G * ) has the same distribution as F and satisfies H |C = F |C .

Fig. 1 .

 1 Fig. 1. Inpainting with the oracle Gaussian model. The masked texture is inpainted with a Gaussian model estimated on the non masked input using a conditioning set of 3-pixel-width along the mask boundary. The kriging component and the innovation component are displayed.As one can see on the result, the conditional inpainting is able to fill both small and large holes, whatever the regularity of the boundary.

Fig. 2 .

 2 Fig. 2. Examples of microtexture inpainting. Here, a Gaussian texture model was estimated on the masked exemplar (middle column) in the red rectangle, and conditionally sampled to get the inpainted texture (right column). The inpainting is satisfactory because the output has similar aspect than the original unmasked texture (left column) while being different on the mask.

Fig. 4 .

 4 Fig. 4. Comparison with Efros-Leung [19]. The original image (top left) was inpainted by learning a Gaussian model in the red rectangle. The kriging component is shown in the top right. In the second row, we show the result of [19] (left) and our result (right). As one can see, our algorithm better preserves the frequency content of the estimated texture, and the filled region naturally blends with the surrounding content.

Acknowledgments. This work has been partially funded by the French Research Agency (ANR) under grant nro ANR-14-CE27-0019 (MIRIAM).

This work has been partially funded by the French Research Agency (ANR) under grant nr ANR-14-CE27-001 (MIRIAM)


Bruno Combourieu

Pascale Pascale Besse

P Besse-Hoggan Besse-Hoggan

Martine Sancelme

Edmund Maser

A M Delort

Evidence of metyrapone reduction by two Mycobacterium strains shown by 1H NMR

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.