Jiaxin Xu

Jérôme Bobin
email: jerome.bobin@cea.fr

Anne De

Vismes Ott

Christophe Bobin

Paul Malfrait

Analysis of gamma-ray spectra with spectral unmixing -Part I: Determination of the characteristic

Keywords: gamma-ray spectrometry, decision threshold, uncertainty, full spectrum analysis, spectral unmixing

come

Analysis of gamma-ray spectra with spectral unmixing -Part I: Determination of the characteristic limits (decision threshold and statistical uncertainty)

Introduction

Gamma-ray spectrometry is one of the main techniques used for the identification and quantification of gamma-emitting radionuclides. The measurements can be performed with different types of detectors, including semi-conductor detectors such as High Purity Germanium detectors (HPGe) and the scintillation detectors, e.g., Sodium Iodide (thallium doped), NaI(Tl). A gamma-ray spectrum is the histogram of the number of detected events as a function of the energy that is deposited by the gamma-ray or the X-ray in the detector. In this context, radionuclides are classically identified and quantified with peak-based or deconvolution-based spectrum analysis methods (e.g.,Mirion Technologies (Canberra)'s Genie 2000 gamma analysis software1 ), which however has limited performances in terms of low-level activity measurements. In recent studies, with the purpose of providing more accurate and sensitive activity measurements, spectral unmixing methods have been investigated for gamma-ray spectrum analysis in [START_REF] Xu | Sparse spectral unmixing for activity estimation in γ-ray spectrometry applied to environmental measurements[END_REF] applied to HPGe measurements and in [START_REF] Paradis | Spectral unmixing applied to fast identification of γ-emitting radionuclides using NaI(Tl) detectors[END_REF], [START_REF] André | Metrological approach of γ-emitting radionuclides identification at low statistics: application of sparse spectral unmixing to scintillation detectors[END_REF] applied to NaI measurements.

In contrast to peak-based analysis only from the full energy peaks, spectral

unmixing is an analysis tool that identifies and quantifies radionuclides from the full energy spectrum (i.e., peaks and their associated continua) of a gammaray measurement. Furthermore, it allows accounting for the Poisson statistics that underlie the detection process. In this framework, gamma-ray spectrum analysis is recast as a regularized inverse problem, with the aim of decomposing a measured spectrum into the individual spectral contribution of each radionuclide present in the measured sample. A measured spectrum x is composed of M channels: x = [x 1 , ...x M ]. For ∀i ∈ [1, ...M ], the Poisson process of radioactive decay leads to model the problem as:

x i ∼ Poisson ([Φa] i + b i ) (1)
where we note the radionuclides' spectral signatures as Φ = [φ 1 , ...φ N ] (containing N radionuclides' spectral signatures) and the background spectrum as b. The spectral unmixing addresses the estimation of radionuclides' mixing weights vector a that are proportional to radionuclides' activities. According to the Poisson statistics, the likelihood of measuring a spectrum x knowing that the spectra linear combination is Φa + b can be written as follows:

P (X = x|Φa + b) = M i [Φa + b] i xi e -[Φa+b]i x i ! (2)
Minimizing the neg-log-likelihood of this distribution leads to the following estimator of the mixing weights [START_REF] Xu | Sparse spectral unmixing for activity estimation in γ-ray spectrometry applied to environmental measurements[END_REF]:

â = argmin a≥0 Φa + b -x log (Φa + b) (3)
where is the Hadamard product and â (i.e. entry-wise product) represents the estimate of the mixing weights vector a. The above optimization problem has been tackled as a regularized inverse problem, where the non-negativity is imposed on the solution to be estimated.

The Poisson-statistics based spectral unmixing has been shown to be a more accurate activity estimation method than peak-based analysis. To use such estimation algorithms for experimental data analysis, quantifying uncertainties is a critical step for assessing the statistical confidence of the estimation procedure as well as for decision making purposes. In radioactivity metrology, it is conventionally called the characteristic limits of a measurement. These limits comprise the decision threshold, the detection limit and the limits of the coverage interval, they are essential to achieve a metrological use of the spectral unmixing procedure. The determination of these limits are defined in the ISO 11929 (2010).

Their assessment has been studied for commonly used peak-based gamma-ray spectrum analysis method. [START_REF] Korun | Calculation of the decision thresholds in gamma-ray spectrometry[END_REF] studies the calculation of the decision thresholds. [START_REF] Boshkova | Decision threshold and relative uncertainty of the measured net signal in radiation measurement[END_REF] investigates the statistical significance of the measured net signal through the determination of characteristic limits. So far, and to the best of our knowledge, no study has considered the full spectrum analysis case. The aim of this work is to explore how the evaluation of the characteristic limits can be generalized to the Poisson-statistics based spectral unmixing procedure.

The paper is organized as follows: Section 2 presents the concept of characteristic limits for radioactivity measurements. In Section 3, we describe the data used in this work. Section 4 explores the calculation of the decision threshold.

Next, we discuss the assessment of the coverage intervals in Section 5. Section 7 provides the conclusions of this work.

Characteristic limits in radioactivity measurements

Definitions

Before going further, we first introduce some useful notations to describe the characteristic limits:

• A: Measurand, the quantity of interest.

• â: Determined value of the measurand Y (i.e., the estimate of A).

• ã: True value of the measurand.

[JB: Pourquoi ne pas mettre les mixing weights en gras, ils le sont plus loin ?] JX: dans ce contexte, il s'agit l'estimation de l'activity d'un radionucleide, i.e., mixing weigt, donc c'est un scalar Quantifying the characteristic limits in gamma-ray spectrum analysis is carried out within the classical statistical hypothesis framework. It considers testing hypotheses with the two following alternatives (associated with type I error and type II error described in Table 1):

• H 0 : the null hypothesis, where a given radionuclide is not "active", i.e.,

namely not present in the sample.

• H 1 : the alternative hypothesis, where the radionuclide is present in the mixture.

The standardization document (ISO 11929, 2010) defines the characteristic limits for ionizing radiation measurements and provides a framework for the computation of these limits. Referring to [START_REF] Weise | Determination of the detection limit and decision threshold for ionizing-radiation measurements: fundamentals and particular applications[END_REF], [START_REF] Michel | Measuring, estimating, and deciding under uncertainty[END_REF], the definition and interpretation of these limits for some estimate â of a measurand

A are as follows:

• Decision threshold (DT) allows a decision to be made on whether or not the physical effect quantified by the measurand is present.

The determination of DT is related to the Type I error described in Table 1. When the quantity y exceeds the critical value (DT), the null hypothesis H 0 should be rejected with respect to a given false positive rate (FPR).

It can be described with:

α = P (â ≥ DT |ã = 0) ( 4
)
where ã is the true value of the measurand and α is the desired critical FPR.

• Detection limit (DL) indicates the smallest true quantity value of the measurand, which can still be detected with the applied measurement procedure.

The determination of DL is related to the Type II error described in Table 1. It is selected with respect to a desired false negative rate (FNR) based on the decision threshold level.

More precisely, the detection limit (DL) is the smallest value that provides a desired Type II error probability β:

β = P (â ≤ DT |ã = DL) (5)
where the DT is given and ã is the true value of the measurand.

• The coverage interval for the estimate â is an interval that has a probability 1 -γ of containing the true value ã:

1 -γ = P (â l < ã < âr ) (6)
where [â l , âr ] stands for the coverage interval.

The regular radioactivity monitoring throughout France is ensured by the operators, regulated by the ASN (French Nuclear Safety Authority) and supplemented by IRSN (French Radiation Protection and Nuclear Safety Institute).

The environmental samples measurements are performed in this framework and the measurements data are fed into the RNM network1 . For that purpose ASN imposes that data contain the activity, the expanded uncertainty and the decision threshold. Indeed the decision threshold serves to conclude if the measurement is statistically significant or not, given a fixed false positive rate α.

The detection limit is not required because it is supposed to be compared with some regulatory levels, while the radioactivity monitoring aims to detect rapidly any increase in environmental radioactivity.

In this paper, we thus focus on the determination of the decision threshold and the coverage interval for Poisson-statistics based spectral unmixing algorithm. Firstly, one can determine whether the resulting activity of a radionuclide is significant by comparing it to the decision threshold. Secondly, the coverage interval provides the statistical uncertainty of the estimation.

Quantification of characteristic limits with Monte Carlo simulations

A traditional approach to quantify the characteristic limits for some estimation method is to make use of Monte Carlo simulations. The main drawback is their massive computational cost since Monte-Carlo simulations are needed for each new spectrum to be analyzed. Furthermore, the true value of the measurand in a measured spectrum is not known in advance. Therefore, we focus on alternatives which are less computationally demanding and yet precise to derive the characteristic limits without resorting to Monte-Carlo simulations. To evaluate and validate these methods, we propose to use synthetic spectra in which the radionuclides' mixtures are perfectly known. This approach was employed because the characteristic limits of these synthetic mixtures can be precisely calculated by Monte Carlo simulations, which allows assessing the statistical uncertainty of the proposed approach.

Decision threshold from Monte Carlo simulations. Let the probability distribution of the mixing weight under null hypothesis of the radionuclide indexed by j be defined as:

P (a j |ã j = 0) (7)
According to Eq.( 4), the decision threshold can be derived from 1 -α percentiles of this probability distribution, for α being the false positive rate with respect to the DT level. In practice, α = 2.5% is usually used in environmental measurements. Assuming that the linear combination model of radionuclides is known, one can make use of Monte Carlo simulations that mimic the mixture under the null hypothesis of each radionuclide to evaluate their DT levels.

In the experiments of this paper, 1000 simulations are performed to calculate an accurate DT value according to α = 2.5% for each radionuclide. In each random simulation process, we set the radionuclide to zero and other radionuclides to their actual levels. The spectral unmixing algorithm is subsequently applied to estimate radionuclides' mixing weights from these simulated spectra. The results obtained for the radionuclide set to zero (i.e., null hypothesis) allow us to compute its decision threshold value.

Coverage interval from Monte Carlo simulations. Since each estimate of the mixing weights is assumed to have a limited estimation bias or error, a coverage interval with respect to a given probability γ for individual radionuclides can also be derived by performing Monte-Carlo simulations. In practice, γ = 5% is usually taken into account. More precisely, we first perform 1000 Monte Carlo simulations that mimic the actual mixing scenario of radionuclides. Next, the spectral unmixing algorithm is applied to estimate mixing weights from these simulated spectra. Finally, the coverage intervals of radionuclides can be derived from the distribution of their estimated mixing weights.

Data description

The main contribution in this paper is the characteristic limits calculation with spectral unmixing analysis, applied to aerosol measurements performed with HPGe detectors. The validated methods are further studied in Xu et al.

(2021) for quantitative analysis of experimental data. In this context, we make use of simulations of aerosol filter measurements of a Mirion Broad Energy Germanium (BEGe s and normalized to the simulation counting time. The activity of the blank filters is disregarded since it has been checked experimentally to be negligible.

Table 2 shows the activities of these radionuclides that are customary in real aerosol measurements (i.e., measured during 320000 s). Since each simulated spectral signature corresponds to the energy response with unit particle (i.e.

one disintegration), the relationship of the radionuclides' mixing weights vector a and their activities can be described as follows (see details in [START_REF] Xu | Analysis of gamma-ray spectra with spectral unmixing -Part II: Calibrations for the quantitative analysis of HPGe measurements[END_REF]):

g = a t ( 8
)
where g is the radionuclides' activities in becquerel (Bq) (i.e., number of disintegrations per second) and t is the counting time.

We focus on the assessment of characteristic limits for 4 radionuclides: 7 Be,

22 Na, 137 Cs, 212 Pb, since these radionuclides cover the whole energy range and different statistical regimes. The simulation model of 10 radionuclides and the contributions of these 4 radionuclides are displayed in Fig. 1.

Next, we also generate simulations of measurements of a 3"x3" NaI(Tl) detector without shielding using point sources placed at a distance of 1 m [START_REF] Paradis | Spectral unmixing applied to fast identification of γ-emitting radionuclides using NaI(Tl) detectors[END_REF], which is made of 1024 channels. The spectral signatures corre-spond to the detector response of 4 gamma-emitting radionuclides with photon emissions covering a range of energies between 40 keV and 2 MeV: 60 Co, 134 Cs, 137 Cs, 152 Eu (see Fig. 2). We make use of simulated spectra with respect to the numbers of counts shown in Table 3.

Quantifying the decision threshold

In peak-based gamma-ray spectrum analysis, the decision threshold is usually derived from some statistical test based on the measured spectrum. This means to evaluating how much the estimated quantity associated with a radionuclide's activity departs from the background (i.e., other contributions composed in the measured spectrum) and is therefore statistically consistent or not with this background.

In the full spectrum analysis context, the spectral unmixing decomposes a gamma-ray spectrum into individual spectra of radionuclides. To determine the decision threshold of a single radionuclide indexed by j in the unmixing model, we reformulate the true linear mixing model with this radionuclide and an equivalent background:

Φa + b → φ j a j + m (9)
where φ j a j represents the contribution of the j th radionuclide in the spectrum.

The equivalent background is then composed of the background spectrum b and the contribution of all but the tested radionuclide:

m = l =j φ l a l + b (10)
Take an example of 137 Cs in the measurement presented in Fig. 1, the spectral contributions are illustrated in Fig. 3.

Recall the definition in Section 2, the decision threshold of the measurand a j (i.e., the mixing weight of the j th radionuclide) is derived from:

α = P (a j ≥ DT |ã j = 0) (11)
The null hypothesis of ãj (i.e., true value of a j ) implies that m is the mean value of the distribution:

x ∼ Poisson (m) (12)
We assume that m is accurately estimated by the spectral unmixing method.

This estimate is denoted by m. The decision threshold of the j th radionuclide can then be derived by the following statistical test:

α = P (t ≥ T (λ)|λ i = mi , ∀i ∈ [1, ...M ]) ( 13
)
where M is the number of channels in the measured spectrum. The quantity T (λ) is some test statistics that is defined as a scalar function of some background model λ. The decision threshold can be derived from the statistical distribution of such a test statistics under the null hypothesis: the tested radionuclide is not present and the observed spectrum is described by the equivalent background only. Next, we will consider different choices for the λ and T (λ). More precisely, we will considered test statistics T (λ) that will be defined of weighted sums of the total number of counts:

T (λ) = M i=1 w i x i (14)
Thanks to the statistical independence of each channel of the measured spectrum, the distribution of such general test statistics can be written by the following generic expression:

P M i=1 w i x i M i=1 w i mi ( 15
)
where w = [w 1 , ...w M ] stands for a vector of weights. The goal of this general approach is that it allows to select weights that are better adapted to distinguish between the radionuclides' contribution in the measured spectrum and its equivalent background. The DT (noted a * j ) of the j th radionuclide with respect to a given false positive rate α, are studied with 4 different choices of w.

1. A simple statistical test to consider is based on the total number of counts, as measured by the sum of observed counts under the null hypothesis ãj = 0, which leads to:

w sum = {w i = 1, ∀i ∈ [1, ...M ]} (16)
Eq.( 15) can be formulated with the following Poisson distribution:

M i=1 x i ∼ Poisson M i=1 mi (17)
The DT, a * j can be derived from its cumulative distribution function (CDF):

α = P M i=1 x i ≥ M i=1 [φ j a * j ] i + M i=1 mi (18)
2. By using w sum , we take into account the information carried out by the full spectrum, it however poorly distinguishes the radionuclide to be tested from the background. For instance, in Fig. 3, the Compton continuum of the 137 Cs is under its equivalent background, while it is better distinguished from the background in the peak region. Therefore, we rather use the pre-specified channels in a region of interest (ROI -e.g, peak region of the radionuclide), which can be written as:

w ROI =      w i = 1, if i ∈ ROI w i = 0, otherwise (19)
Indeed, this approach is carried out in the same manner as peak-based analysis, in which the DT level is derived based on the number of background counts (ISO 11929, 2010) for more details).

We calculate the DT from the CDF of the following distribution:

i∈ROI x i ∼ Poisson i∈ROI mi (20)
3. In order to better distinguish between the radionuclide to be tested and its equivalent background in the full energy range, we further investigate the DT calculation based on the measured spectrum x weighted by the projection of x onto φ j parallelly to m, which is graphically illustrated in Fig. 4. In constrast to a standard orthogonal projection, a projection onto φ j parallely to m allows to measure the contribution of x explained by the signature φ j while guaranteeing its contribution along the direction m to be zero. This eventually enhances the contrast between the contribution of the j-th radionuclide and the equivalent background m in

x, thus providing a more sensitive evaluation of the detection threshold.

Such a parallel projection is denoted by w pb .

In pratice it can be evaluate by first defining the 2-columns matrix Ψ = φ j m : the spectral signature φ j and the estimated equivalent background m. A simple weights vector to be considered is:

w pb = {(Ψ T Ψ) -1 Ψ T } φ j (21)
where {} φ j defines the entry related to φ j .

However, the resulting statistical test does not follow a Poisson distribution. Thankfully, since it is defined as a linear combination of a large number of observed channels, one can approximate the distribution to be Gaussian. The resulting statistical test with observed counts x and estimated equivalent background m will be calculated as follows:

M i=1 [w pb ] i x i ∼ N M i=1 [w pb ] i mi , M i=1 [w pb ] 2 i mi (22)
while in each channel, the mean value µ = [w pb ] i mi and variance σ 2 = [w pb ] 2 i mi are considered due to the Poisson statistics of the model. The DT, a * j can be derived from the CDF of the above distribution:

α 2 = P M i=1 [w pb ] i x i ≥ M i=1 [w pb ] i [φ j a * j ] i + M i=1 [w pb ] i mi (23) 4.
The above statistical test assumes the underlying noise is additive, white and Gaussian. To consider the covariance matrix underlying the Poisson noise, we further make use of an approach based on the a re-weighted projection parallel to the background, noted w rpb :

w rpb = Ψ T diag 1 φ j âj + m Ψ -1 Ψ T diag 1 φ j âj + m φ j (24)
and calculate the DT from the CDF of:

M i=1 [w rpb ] i x i ∼ N M i=1 [w rpb ] i mi , M i=1 [w rpb ] 2 i mi (25)
Experimental results

These four proposed statistical tests are applied to assess the decision thresh- For HPGe spectra, the different statistical tests are compared in Fig. 5, 6, 7 and 8 for 7 Be, 22 Na, 137 Cs and 212 Pb respectively. For NaI spectra, the statistical tests are compared in Fig. 9, 10, and 11, and 12 for 60 Co, 134 Cs,

137 Cs and 152 Eu respectively. The results show the distributions of:

• the accurate DT level calculated from Monte Carlo (i.e., 1 -α percentile of the distribution under the null hypothesis) is displayed with a dotted line (red).

• Distribution of calculated DT based on the statistical test of equivalent background weighted by w sum (red), w ROI (green), w pb (blue) and w rpb (orange).

From these results, we can draw the following conclusions:

• Firstly, the statistical test based on the equivalent background weighted by w sum is less precise to derive the DT than other approaches, since it poorly distinguishes the radionuclide's spectrum from the background.

• Compared to statistical tests based on the projected equivalent background (i.e., weighted by w pb and w rpb ), the ROI-based statistical test(i.e., weighted by w ROI ) is less precise since it only uses information at the vicinity of the peak regions.

• For HPGe measurement, for 22 Na, 137 Cs and 212 Pb, statistical tests based on the equivalent background weighted by w rpb derives similar DT levels compared to the actual 1 -α percentile of the distribution (i.e., accurate DT level), which is better than the choice of w pb . This is not valid for 7 Be due to the fact that 7 Be is predominant in the measured spectrum, while its equivalent background is poorly estimated. However, the determination of an accurate DT level for low-level radionuclides is more important for decision making. Indeed, in such a case, the dominant spectral contribution, such as 7 Be, generally yield a more precise estimation of the equivalent background, which entail a more accurate DT determination for low-level radionuclides.

• When we further focus on the results of NaI measurements, the choice of using w rpb is shown to be more consistent than w pb . It allows to better distinguish a radionuclide from its equivalent background when the spectra have a significant overlap in the whole energy range.

The DT assessment is evaluated with Monte Carlo simulations of realistic measurements of aerosol samples performed with HPGe detectors. The proposed analysis involving a statistical test based on re-weighted projections of the equivalent background is shown to be an efficient and rapid tool to determine DT. We apply this procedure to real measurements in the companion article [START_REF] Xu | Analysis of gamma-ray spectra with spectral unmixing -Part II: Calibrations for the quantitative analysis of HPGe measurements[END_REF].

Coverage interval

Evaluating the statistical reliability of a parameter estimation procedure is generally composed of two elements: the estimation bias and an interval of values that contains the true value of the parameter up to some probability.

In the spectral unmixing approach, the statistical uncertainty assessment can be carried out from Monte Carlo simulations that mimic the actual mixture of radionuclides as presented in Section 2.2. The distribution of estimated values provides both the estimation bias and the coverage interval. However, this requires a massive amount of simulations, making this approach hard to implement in practice. In this study, we focus on the metrological use of the Poissonbased spectral unmixing for analyzing the described HPGe detectors performed aerosol measurements. In [START_REF] Xu | Sparse spectral unmixing for activity estimation in γ-ray spectrometry applied to environmental measurements[END_REF], we showed that estimation biases with spectral unmixing analysis are low for routine aerosol measurements.

In this context, we propose to make used of the Fisher information matrix of the estimated mixing weights to compute the coverage interval. This allows for an analytic expression of the coverage interval, which does not require computationally expensive Monte Carlo simulations. Such method has been also investigated in [START_REF] André | Metrological approach of γ-emitting radionuclides identification at low statistics: application of sparse spectral unmixing to scintillation detectors[END_REF] for the coverage intervals assessment of spectral unmixing approach applied to NaI measurements.

For the observed variable x distributed as f (x|θ), the asymptotic distribution of the maximum likelihood estimator (MLE) of the parameter θ is a Gaussian distribution:

N θ, I(θ) -1 ( 26
)
where I(θ) is the Fisher information [START_REF] Fisher | Statistical methods and scientific inference[END_REF]) defined as:

I(θ) = E θ ∂ 2 log f (x|θ) ∂θ 2 (27)
where E θ stands for the expectation taken with respect to f (x|θ). The true value of the parameter is generally not known. When the estimation bias | θ -θ| is assumed to be low, the estimation of standard error of θ M LE can be obtained by replacing the unknown true value θ by the estimated value θ. To validate the metrological use of such an approach, we will compare the coverage interval obtained from the information matrix associated with θ to the one we would obtain from Monte Carlo simulations.

Recall that the Poisson-based spectral unmixing provides a maximum likelihood estimate of the mixing weights, noted â. We propose a first approximation by considering that the estimator of the mixing weights vector is non-biased and follows a Gaussian distribution centered at â. According to the Poisson likelihood Eq.( 3), the covariance matrix is approximated with the inverse of the Fisher information matrix:

I(â) -1 = Φ T diag x (Φâ + b) 2 Φ -1 (28)
where is the element-wise division operator. The coverage interval of the estimated mixing weights â can be derived from the diagonal elements of the inverse of the Fisher matrix; the resulting standard deviation for âi is then defined as [I(â) -1 ] i,i , which approximates the standard deviation of the maximum likelihood estimator distribution. This approximation will adapted when estimated activity is high enough so that the coverage interval does not contain 0. In Section 6, we discuss how the proposed approach can be extended to this case.

Experimental results

The coverage intervals assessment with Fisher information matrix is evaluated with the simulations of HPGe and NaI measurements described in Section 3. To validate this approach, the Gaussian distribution approximated with 4 and Table 5 for HPGe and NaI measurements respectively, which confirms that for this measurement scenario, the coverage intervals computed from the Fisher information matrix are quite accurate.

The results in Table 4 show that the quantiles of the estimation distribution of Monte Carlo simulations within: a 0 ±σ f , a 0 ±2σ f and a 0 ±3σ f are consistent with their expected values for 7 Be, 137 Cs and 212 Pb, whereas the results for

22 Na show a slight deviation. As shown in Fig. 1, the spectrum of 22 Na has a significant continuum, which is impacted by the background. It should be noted that the coverage interval assessment with Fisher information matrix has limited performance when the number of counts is low. This is further confirmed with results of NaI simulations shown in Table 5, in which the number of counts of radionuclides is low. In particular, with a number of counts of only 500, the coverage interval evaluation of 134 Cs from Fisher information matrix has limited performance. The limited performances of using the Fisher information matrix is due to: i), the larger estimation bias at low statistics which is shown in recent studies [START_REF] Xu | Sparse spectral unmixing for activity estimation in γ-ray spectrometry applied to environmental measurements[END_REF]. ii), the underlying optimization problem of the spectral unmixing is a maximum likelihood estimation with non-negativity constraint which makes the Gaussian distribution no longer true at low statistics.

Discussion about the limitations of the proposed methods

In this paragraph, we discuss some potential limitations of the proposed approach.

Statistical limitation

The derivations of the detection threshold and the coverage interval are both based on assumptions that we discuss hereafter:

• Computation of the decision threshold : In the proposed approach, the detection threshold for a given radionuclide is defined from the statistical distribution of some test statistics T that depends upon an equivalent background. This equivalent background is function of the contribution of all the radionuclides except the tested one. Since the such an equivalent background depends upon the result of the unmixing process, it is implicitly a function of the estimated activity of all the radionuclides. However, it is important to highlight that this dependency is weak. Indeed, for each radionuclide to be tested, the equivalent background derives from the (weighted) sum of the contributions of all radionuclides except the tested one. These linear combinations are dominated by the radionuclides with the highest activities, which are generally accurately estimated thus providing a reliable equivalent background estimation.

Traditional peak-based methods also depend on the estimation of an equivalent background at the vicinity of the energy peaks of the tested radionuclides. In the proposed method, the equivalent bakground estimates bear more physical information as the method allows to account for the compton continuum of each radionuclide in the full observed energy range.

• Determination of the coverage interval In the context of spectral unmixing, coverage intervals can be determined in several ways. One standard way would be to perform Monte-Carlo simulations, assuming some prior model about the mixture (e.g. the maximum likelihood estimate of the mixture weights). Similarly to the previous point dedicated to the detection threshold, the resulting coverage intervals will weakly depend on the estimates of the mixing weights. While Monte-Carlo simulations allow to determine potentially complex, non-symetric, coverage intervals, they also require a significantly large computational cost to derive reliable coverage intervals, which might limit its applicability for fast analysis purposes. To that end, we rather proposed to opt for an approximation of the distribution of the mixing weights distribution based on the Fisher information matrix. The inverse of the Fisher information matrix provably provide the asymptotic covariance matrix of the mixing weights distribution about the maximum likelihood estimate. The quality of this approximation will strongly depend on the actual statistics of the data, at very low statistics, this approximation is unlikely to be valid. Similarly, for radionuclides with very low activities, the corresponding coverage intervals are likely to be non-symetric.

Determining coverage intervals only make sense when for radionuclides that have activities larger than the detection thresholds. In the aforementioned limit cases, we made the assumption that the measured activities are high enough so that the coverage interval does not include 0. This is generally not the case for radionuclides whose activities are close to the detection threshold. In this case, the non-negativity constraint used in the spectral unmixing method prevents the estimated mixing weights to be negative. Consequently, the distribution of the estimated mixing weights can be described by a truncated Gaussian distribution defined as follows:

f (θ) = pδ(0) + Π + (θ) 1 2πσ 2 i exp - (θ -a M L i ) 2 2σ 2 i (29)
where δ(0) is the Dirac function centered about 0 and the probability weight at 0 is defined as:

p = 0 -∞ 1 2πσ 2 i exp - (θ -a M L i ) 2 2σ 2 i dθ
The function Π + (θ) is a step function that takes the value 0 when θ < 0 and 1 when θ > 0. The coverage interval for a given probability α can then be defined as the smallest interval [σ 1 , σ 2 ] so that:

σ2 σ1 f (θ)dθ = α (30)
As an illustration, we perform a comparison of the above analytic model with Monte-Carlo simulations for 137 Cs with an activity taken at the detection threshold. Figure 14 shows the distribution obtained from 15000

Monte-Carlo simulations in blue and the analytic Fisher-based approximate distribution defined in 29, which show a very good agreement. This experiments highlights that the Fisher-based approximate distribution can be generalized to low-level cases, close to the detection threshold.

Modeling limitations

[AVO: A RELIRE. TROP LONG ???]

In standard gamma-ray spectrometry analysis the detection efficiency calibration is required for each radionuclide and each detection configuration (specific detector and counting geometry). Alike the determination of the spectral signatures is required for each radionuclide and each detection configuration in the spectral unmixing approach. Due to the unavailability of standard sources containing all the radionuclides of interest in single emitter standard source, the spectral signatures are calculated using Monte Carlo type particle transport codes. MCNP-CP simulates the interaction between the photons emitted by the radionuclide of the source (or the sample) and the germanium crystal, and the surrounding materials. In this paper the simulations were performed considering the detection system used for the aerosol filters measurements : the same HPGe detector (61% Mirion-Canberra Broad Energy Germanium detector), the same counting geometry (10 mL) and the same sample matrix (polypropylene fibers). The sample height or position variability will be taken into account in the uncertainty budget described in the associated paper [START_REF] Xu | Analysis of gamma-ray spectra with spectral unmixing -Part II: Calibrations for the quantitative analysis of HPGe measurements[END_REF]) Jiaxin et al. (2021). Another component is also added in the uncertainty budget to account for the simulation process : indeed the modeling of the detection configuration is adjusted to reproduce the full energy peak (FEP) efficiency which is experimentally determined by the measurement of a standard source whose activities are given with uncertainties. Neither the total detection efficiency nor the full spectrum shape can be compared to experimental data, except with single radionuclide source measurements. This lack of knowledge is certainly a limitation of the proposed approach in the case when the spectral signatures must be determined by simulation. Forthcoming study will be dedicated to the potential differences between simulated and experimental spectral signatures, using single gamma-ray emitter standard sources.

Conclusion

In this paper, we investigate how the characteristic limits, i.e. the decision thresholds and coverage intervals, can be determined for spectral unmixing algorithms. These two quantities, which are related to the statistical limits for decision making purposes in gamma ray spectra analysis, are first investigated for such full spectrum analysis methods.

  old for HPGe and NaI measurements described in Section 3. The results are compared to those carried out with Monte Carlo simulations, noted as MC (described in Section 2.2). The false positive rate is fixed to α = 2.5%. More precisely, the proposed statistical tests are used to assess the decision threshold of radionuclides in the simulations of the HPGe and NaI measurements for each Monte Carlo simulation. The statistical test with w ROI is only used for HPGe measurements. Indeed, NaI measurements have much lower spectral resolution, which makes it much harder to define an accurate ROI. Experiments on the described selection of gamma-ray emitters (see Fig. 1) of the HPGe measurement considers the following peak regions: 7 Be at 477 keV, 22 Na at 1274 keV, 137 Cs at 662 keV and 212 Pb at 238 keV.

Fisher

  information matrix associated with the estimated mixing weights need to agree with the expected distribution (i.e., estimated distribution obtained from Monte Carlo simulations). Firstly, we make use of the Q-Q (quantile-quantile) plots of the Normal distribution generated from the standard deviation calculated with the Fisher information and the distribution of estimations derived from the Monte Carlo simulations, by plotting their quantiles one against another. The aim is to i), test the normality of the estimator, which allows val-idating the coverage interval with a standard uncertainty in terms of Normal distribution. ii), compare the standard deviation of each estimation calculated with the Fisher information to those obtained with Monte Carlo simulations.More precisely, the Q-Q plots are shown in Fig.13and Fig.15for the HPGe measurement and the NaI measurement, respectively.The Q-Q plots allows comparing the distribution of the scales and centered estimated parameters with a centered normal distributed, which is its expected asymptotic distribution. The straight lines show that taling the asymptotic Normal approximation provides a very accurate evaluation of the distribution of the estimated mixing weights, at least in the interval [0, 3]. [JB: taling dans ce context ?] Additionally, we can asses how accurate are the coveral intervals derived form the Fisher information matrix. For that purpose, we perform Monte-Carlo simulations to compute the probability associated with the intervals : a 0 ±σ f , a 0 ±2σ f and a 0 ±3σ f , where a 0 , σ f stand for the expected mixing weight (i.e., simulated value), which we comapre with the error determined by the Fisher matrix. The corresponding percentages are expected to be 68.27 %, 95.45 % and 99.73 % for Normal distribution. The results are shown in Table

Firstly, we propose

  to determine the decision threshold of each radionuclide based on different statistical tests of its equivalent background. This background spectrum, consisting of other spectral contributions provides a mean value of the distribution under null hypothesis, allowing quantifying the decision threshold based on a given false negative rate. Next, we propose to assess the limits of coverage intervals with the Fisher information matrix of estimated mixing weights, which approximates the standard deviation of the estimation distribution. This statistical uncertainty must be added to other uncertainty components which are introduced and described in[START_REF] Xu | Analysis of gamma-ray spectra with spectral unmixing -Part II: Calibrations for the quantitative analysis of HPGe measurements[END_REF] to obtain the final uncertainty.

Fig. 1 :

 1 Fig. 1: Spectral unmixing model used to illustrate the evaluation of the characteristic limits.The simulated measurements (gray) are composed of 10 radionuclides: 7 Be, 22 Na, 40 K, 137 Cs,

Fig. 4 :

 4 Fig. 4: [JB: Figure à refaire, la projection sur m ne devrait pas être // à gauche] ist of parallel projection with respect to the equivalent background. The vector φ j is the spectral signature of the radionuclide. The vectors x, m are the measured spectrum and the equivalent background spectrum. w pb is the weight vector.

Fig. 5 :

 5 Fig. 5: Decision threshold assessment for HPGe measurements: 7 Be. In Fig. 5, 6, 7 and 8: the right plot shows the distribution of calculated DT based on the statistical test of equivalent background weighted by wsum (red), w ROI (green), w pb (blue) and w rpb (orange). The accurate DT level calculated from Monte Carlo (i.e., 1 -α percentile of the distribution under the null hypothesis) is displayed with a dotted line (red). The left panels show the zoom in of the distributions in the right plot.

Fig. 6 :

 6 Fig. 6: Decision threshold assessment for HPGe measurements: 22 Na.

Fig. 7 :

 7 Fig. 7: Decision threshold assessment for HPGe measurements: 137 Cs.

Fig. 8 :Fig. 9 :

 89 Fig. 8: Decision threshold assessment for HPGe measurements: 212 Pb.

Fig. 11 :

 11 Fig. 11: Decision threshold assessment for NaI measurements: 137 Cs.

Fig. 12 :Fig. 13 :

 1213 Fig. 12: Decision threshold assessment for NaI measurements: 152 Eu.

Fig. 14 :

 14 Fig. 14: Evaluation of coverage intervals assessment with Fisher information matrix: HPGe measurements with 137 Cs at decision threshold level. The distribution of estimated mixing weight of 137 Cs of 15000 Monte Carlo simulations (blue) is compared to the probability density function of the truncated normal distribution N a 0 , σ f .

  T M ) (61 % relative efficiency) investigated in Xu et al. (2020), covering the energy range from 20 keV to 1640 keV composed in 16384 channels. The choice of gamma-ray emitters consists of 10 radionuclides: 7 Be, 22 Na, 40 K, 137 Cs, 210 Pb, 208 Tl, 212 Bi, 212 Pb, 214 Bi, 214 Pb. The background radiation spectrum is according to a measurement without sample during 560000

Table 1 :

 1 The proposed approaches have been evaluated with Monte Carlo simulations of both HPGe and NaI detectors. They have been shown as efficient and reliable methods to both of these simulated data. As the simulated spectra are calculated for radionuclides in samples measured on the HPGe detector of the routine aerosols measurements in the laboratory, the methods proposed in this paper can be used to provide metrological results of experimental data, which is presented in[START_REF] Xu | Analysis of gamma-ray spectra with spectral unmixing -Part II: Calibrations for the quantitative analysis of HPGe measurements[END_REF]. Two types of errors of hypotheses test. Cs 210 Pb 208 Tl 212 Bi 212 Pb 214 Bi 214 Pb

		H 0 is true	H 1 is true
	rejecting	Type I error: the error of
	H 0	rejecting H 0 when it is true,
		the probability of committing
		a type I error is denoted by α,
		called false positive rate.
	accepting		Type II error: the error of
	H 0		accepting H 0 when H 1 is true,
			the probability of committing
			a type II error is denoted by β,
			called false negative rate.
	7 Be 137 391 22 Na 40 K 0.078 1.24 0.015 27	0.098 0.425 0.239 0.054 0.056

Table 2 :

 2 Activity (Bq) of radionuclides of HPGe simulations.

	Radionuclide	60 Co 134 Cs 137 Cs 152 Eu background
	Number of counts 1500	500	1500	2500	4000

Table 3 :

 3 Number of counts of NaI simulations.

		7 Be	22 Na 137 Cs 212 Pb expected value
	percentage within a 0 ± σ f	68.0	67.3 68.87 68.25	68.27
	percentage within a 0 ± 2σ f 96.49 94.28 95.85 95.02	95.45
	percentage within a 0 ± 3σ f	99.9	99.8 99.76	99.9	99.73

Table 4 :

 4 Standard deviation from Fisher information matrix comparing to Monte Carlo simulations (HPGe measurements). 24 60 Co 134 Cs 137 Cs 152 Eu expected value percentage within a 0 ± σ f 68.38 65.81 68.7 65.58 68.27 percentage within a 0 ± 2σ f 95.34 94.78 94.62 95.23 95.45 percentage within a 0 ± 3σ f 99.73 99.54 99.95 99.8 99.73

Table 5 :

 5 Standard deviation from Fisher information matrix comparing to Monte Carlo simulations (NaI measurements).

https://www.mirion.com/products/genie-2000-gamma-analysis-software

https://www.mesure-radioactivite.fr/en#/expert

Pb, 208 Tl, 212 Bi, 212 Pb, 214 Bi, 214 Pb, while individual spectra of 7 Be, 22 Na, 137 Cs, 212 Pb will be evaluated.

Fig. 10: Decision threshold assessment for NaI measurements: 134 Cs.


Brigitte Jaumard

Yan Wang

David Coudert

Dantzig-Wolfe Decomposition for the Design of Filterless Optical Networks

Filterless optical networks use passive splitters and combiners with coherent optics, providing wavelength selection in the digital domain, while forming a passive fiber-tree topology between nodes. In this study, we investigate the optimal design of filterless optical networks while minimizing the number of required wavelengths.

We propose a Dantzig-Wolfe decomposition model in which each subproblem aims to generate a potential filterless optical subnetwork, with a directed tree topology. The master problem then selects the best combination of subnetworks.

Numerical experiments demonstrate significant performance improvement over previous work, reducing previous computational results by a factor of 2 to 10 depending on the size of the data instances.

I. INTRODUCTION

Optical networks are widely deployed today in all kinds of telecommunication networks, with a volume of user traffic on carrier networks growing at a rate between 20% and 40% per year during the last decade [START_REF]Cisco Annual Internet Report[END_REF]. As exponential traffic growth was set to slow, the rise of 5G introduced new use cases, with capacity growth set to continue. In this context, from a techno-economic point of view, filterless optical networks offer opportunities for cost and energy reduction, at least for certain types of networks, for example access networks [START_REF] Pavon-Marino | Techno-economic impact of filterless data plane and agile control plane in the 5G optical metro[END_REF].

There is a compromise to be found between the usage of more spectrum efficient technologies thanks to, e.g., pure optical layer switching using expensive Reconfigurable Optical Add-Drop Multiplexers (ROADMs) (see, e.g., Simons [START_REF] Simmons | Optical Network Design and Planning[END_REF]), and filterless optical networks which rely on cheap couplers/splitters at the expense of a reduction in the spectral efficiency performance [START_REF] Shariati | Real-time optical spectrum monitoring in filterless optical metro networks[END_REF].

Tzanakaki et al. [START_REF] Tzanakaki | Analysis of filterless wavelength converters employing cross-gain modulation in semiconductor optical amplifiers[END_REF] was one of the pioneers to investigate the filterless network components. The basic design principle of filterless optical networks was described by Savoie et al. [START_REF] Savoie | Physical layer validation of filterless optical networks[END_REF] in the early work on filterless optical networks with a motivation for cost reduction and environment friendly networks. In spite of a decrease in the spectrum efficiency, more authors are studying filterless optical networks today, and we can see studies investigating the wasted spectrum on the static planning case (filterless) vs. dynamic planning case (filtered node), see, e.g., Pavon-Marino et al. [START_REF] Pavon-Marino | Techno-economic impact of filterless data plane and agile control plane in the 5G optical metro[END_REF]. The addition of blockers offers a mean to improve the efficiency of spectrum utilization, see, e.g., Dochhan et al. [START_REF] Dochhan | Flexible metro network architecture based on wavelength blockers and coherent transmission[END_REF].

Several researchers study the design of optimized filterless optical networks. Tremblay et al. [START_REF] Tremblay | Passive filterless core networks based on advanced modulation and electrical compensation technologies[END_REF], Ayoub et al. [START_REF] Ayoub | Filterless and semi-filterless solutions in a metro-HAUL network architecture[END_REF] and Jaumard et al. [START_REF] Jaumard | Optimal design of filterless optical networks[END_REF] all proposed a two-step solution process with, for the first two references, a first step aiming at generating potential/promising fiber trees, and then a second step to provision the resulting trees. The quality of the resulting filterless optical networks is however difficult to assess, both for the two-step process, but also because one or both steps are solved using heuristics. For instance, Tremblay et al. [START_REF] Tremblay | Passive filterless core networks based on advanced modulation and electrical compensation technologies[END_REF] use heuristics for both steps: a genetic algorithm for the fiber link interconnexion (definition of filterless trees and subnetworks), and a Tabu Search meta-heuristic for the wavelength assignment.

In [START_REF] Jaumard | A two sub-problem decomposition for the optimal design of filterless optical networks[END_REF], there is a first attempt for designing an exact method for the optimal design of filterless optical networks, using a decomposition mathematical program. However, computational times are in the order of hours. In this paper, we improve on that latter study, with an improved mathematical programming model, as well as an improved exact solution process. This translates into significant reduced computational times and further improved accuracy for the output designs.

More recently, other studies have addressed additional features of filterless networks, e.g., survivability [START_REF] Ayoub | Survivable virtual network mapping in filterless optical networks[END_REF], bidirectional single-fiber filterless optical networks [START_REF] Uzunidis | Bidirectional single-fiber filterless optical networks: modeling and experimental assessment[END_REF], programmable filterless networks [START_REF] Ayoub | Trafficadaptive re-configuration of programmable filterless optical networks[END_REF].

The paper is organized as follows. In Section II, we present the problem statement and the concept of filterless optical network design. In Section III, we propose a new mathematical model DW FOP for the optimal design of filterless optical networks. The solution process, which relies on Dantzig-Wolfe decomposition (or column generation techniques), is described In Section IV. Computational results are summarized in Section V, including a comparison with the results of previous works. Conclusions are drawn in the last section.

II. DESIGN OF FILTERLESS NETWORKS

Consider an optical network represented by its physical network G = (V, L), where V is the set of nodes (indexed by v), and L is the set of fiber links (indexed by ). We denote the set of incoming and outgoing links of v as IN(v) and OUT(v), respectively.

The traffic is described by a set of unit requests K, where each request k is characterized by its source and destination

V 1 V 2 V 4 V 3 (a) Physical network V 1 V 2 V 3 Splitter Combiner V 4 FSN 2 FSN 1 (b) Two FSNs V 1 V 2 V 4 V 3 (c) Wavelength provisioning of FSN 2 Fig. 1: Construction of a filterless network solution V 1 V 5 V 4 V 3 V 2 V 8 V 7 V 6 k 1 k 3 k 4 k 2 k 1 k 3 k 2 k 4 Wavelength conflict graph k 5 k 5
Broadcast part of k 3 Fig. 2: Requests with wavelength conflicts nodes, v s and v d respectively. Indeed, non unit requests are considered as a set of unit requests.

The architecture provides at least one optical path between each node pair [START_REF] Pavon-Marino | Techno-economic impact of filterless data plane and agile control plane in the 5G optical metro[END_REF] and is expected to offer the following significant advantages [START_REF] Simmons | Optical Network Design and Planning[END_REF]: 1. Eliminating or reducing the usage of active photonic switching elements. 2. Reducing initial deployment costs. 3. Ease of maintenance and reconfigurability, good resilience, and multicast capability.

Filterless networks corresponds to passive networks in which passive optical splitters and combiners are used for interconnecting fiber links. Their advantages lie in: (i) the elimination/reduction of the usage of active photonic switching elements, (ii) the reduction of the initial deployment costs and (iii) the ease of of maintenance and reconfigurability, good resilience, and multicast capability.

Due to signal broadcast and the absence of filtering, some of the signals appear on unintended links which can lead to an overhead in spectrum usage and an increase of the blocking probability. Indeed, although blockers can be used, the unfiltered channels propagate after reaching their destination, the signal from each transmitter is broadcasted to all branches of the corresponding fiber tree (see definition below) downstream of the source node.

A filterless network solution consists of a set of filterless optical sub-network (FSN for short) solutions, where each FSN relies on an undirected tree, called the backbone of an FSN, deduced from a subgraph of the physical network G. FSNs need to be pairwise fiber link disjoint in order to avoid closed loops, i.e., undesirable laser effects in optically amplified links [START_REF] Tremblay | Passive filterless core networks based on advanced modulation and electrical compensation technologies[END_REF].

Wavelength assignment follows the usual rule of avoiding wavelength conflicts, see the discussion below. On the undirected tree, we use the term "edge" to represent an undirected link, which can host multiple directed links or links for short, possibly in both directions. We denote the set of adjacent links of v as ω(v).

Figure 1 illustrates the construction of a filterless optical network made of two FSNs: the red one (FSN 1 ) is supported by the undirected tree T 1 = {{v 1 , v 3 }, {v 3 , v 4 }}, and the green one (FSN 2 ), by the undirected tree

T 2 = {{v 1 , v 2 }, {v 2 , v 3 }, {v 2 , v 4 }}.
In Figure 1(b), fiber links are represented by directed arrows, where a bidirected arrow indicates two fiber links, one in each direction. The first FSN, i.e., FSN 1 , is made of one bi-directional arrow, and a single fiber link on the second tree edge supporting the FSN. The second FSN, i.e., FSN 2 has two opposite fiber links on each edge of its supporting undirected tree.

On the FSN intermediate nodes, passive splitters and combiners are required to interconnect the fiber links supporting the request provisioning, see, e.g., node v 2 requires three splitters/combiners with two ports in order to accommodate the provisioning of the requests using the green filterless network, see Figure 1(b).

On our example, each request can be provisioned on any one of the sub-networks. While this is not necessarily the case for all requests in general, this may happen in practice for some requests. In the example of Figure 1(b), the request from node v 1 to node v 4 can be provisioned either on FSN 1 or FSN 2 . The decision is usually driven by the objective, e.g., number of wavelengths, number of fiber trees, number of passive optical splitters. In this study, we are interested in establishing a set of FSNs such that the provisioning of all the requests minimizes the number of wavelengths.

Figure 1(c) depicts a possible wavelength assignment, with the provisioning of 4 requests on FSN 2 in Figure 1(b): k 1 →3 from v 1 to v 3 with the red wavelength (plain lines) and its broadcast effect (dashed lines), and with the same notations, k 1 →2 with the blue wavelegnth, k 4 →1 with the blue wavelength, and k 3 →2 with the red wavelength. Note that in the broadcast part, we can allow wavelength re-use, e.g., the red wavelength on link (v 2 , v 4 ) (a first broadcast effect from k 1 →2 , and a second one from k 3 →2 ). We now consider the broadcast effect. When the destination is reached, the lightpath will continue downstream. For instance, when the optical signal supporting request k 1 →2 reaches v 2 through the fiber link (v 1 , v 2 ), the signal is broadcast through fiber links (v 2 , v 3 ) and (v 2 , v 4 ).

Wavelength assignment should be done in order to avoid wavelength conflicts. Such conflicts occur for two requests if there is a shared fiber link on one of their routing path, i.e., between source and destination nodes. However, requests that share links for broadcasting beyond the destination, are not considered conflicting and may share the same wavelength. For instance, in Figure 2, there are conflicts for the requests k 1 and k 4 because their routes between the source and destination nodes share a fiber link. Requests k 1 and k 2 are also wavelength conflicting as one of the fiber links of k 2 lies on the broadcast route of k 1 with unfiltered channels. Similarly, requests k 3 and k 4 are wavelength conflicting, as k 4 is on the broadcast of k 3 . Non conflicting requests are for instance k 1 and k 5 , or k 4 and k 5 . Note that once an FSN is defined, wavelength assignment can be done with a graph coloring algorithm (see, e.g., [START_REF] Tremblay | Passive filterless core networks based on advanced modulation and electrical compensation technologies[END_REF]), using a wavelength graph conflict, as illustrated in Figure 2.

In Figure 1(c), requests k 1 →3 and k 3 →2 are not conflicting and so can be assigned the same wavelength. The same holds between requests k 1 →2 and k 4 →1 . For this example, 2 wavelengths are sufficient to satisfy the 4 requests.

III. MATHEMATICAL MODEL A. Definitions and Notations

We introduce the concept of FSN configuration as a provisioned directed tree characterized by:

• a F = 1 if link is in configuration F, 0 otherwise. • x F k = 1 if request k is routed on configuration F, 0 otherwise. Indeed, if it
happens that a request is provisioned in two different FSNs due to constraints (3), then one of them is arbitrarily chosen (using a post-processing). Note that the reason of the inequalities in constraints ( 3) is motivated by the greater easiness of solving the mathematical program with inequalities, than with equalities. Expressing constraints (3) with equality would eliminate this issue, at the expense of longer computational times.

• β F λ = 1 if wavelength λ is used in configuration F, 0 otherwise.
Denote by F the overall set of FSN configurations, and let n F be an upper bound on the number of FSNs in the network.

B. Mathematical Model

We propose an integer linear program that has two sets of decision variables:

• z F = 1 if FSN configuration is selected as a filterless subnetwork, 0 otherwise. Note that the conflicts between the requests routed on F have been taken into account in the definition of F and the wavelength assignment is encoded in F.

• x λ = 1 if wavelength λ is used by any of the selected FSN configurations, 0 otherwise. The objective corresponds to the minimization of the total number of wavelengths.

min λ∈Λ x λ (1)
subject to:

F∈F β F λ z F ≤ n F x λ λ ∈ Λ (2) F∈F x F k z F ≥ 1 k ∈ K (3) F∈F a F z F ≤ 1 ∈ L (4) F∈F z F ≤ n F (5) z F ∈ {0, 1} F ∈ F (6) x λ ∈ {0, 1} λ ∈ Λ. (7)
Constraints [START_REF] Pavon-Marino | Techno-economic impact of filterless data plane and agile control plane in the 5G optical metro[END_REF] indicate that a wavelength is used in the solution if it is used by any of the selected FSN. Since the FSNs are link disjoint, the same wavelength λ can be used by several FSNs. The constraint is presented as a "big-M" constraint, that is a constraint in which the right-hand side is a large enough constant (i.e., that is always larger than the value of the left-hand side) times a binary variable. Here, the value n F is sufficient since Constraint (5) limits the number of FSNs to be selected to n F . Constraints (3) ensure that each request is routed on at least one selected FSN. Observe that, although a request should be routed on "a unique" FSN, it is faster in practice for LP/ILP solvers to find a solution when the constraint is relax to "at least one" FSN, and then to use a simple post-processing procedure to extract a unique route per request from the solution. Constraints (4) ensure that selected FSNs are link disjoint (i.e., that a link can be used by at most one selected FSN). Constraint (5) limits the number of FSNs to be selected. Finally, Constraints ( 6) and ( 7) define the domains of the variables.

As can be observed, the above model ( 1) -( 7) has an exponential number of variables, and corresponds to a Dantzig-Wolfe decomposition model, which can be solved by a column generation algorithm, and, e.g., another algorithm to derive an ILP solution, without the requirement of explicitly enumerating all the variables.

IV. SOLUTION PROCESS

We propose to use a solution process, in which the column generation technique is used to solve the linear relaxation of problem ( 1)- [START_REF] Dochhan | Flexible metro network architecture based on wavelength blockers and coherent transmission[END_REF]. Section III presented the mathematical model of the so-called master problem in the context of Dantzig-Wolfe decomposition associated with column generation algorithm. We proposed two pricing models to compute the new improved FSN configuration. In section IV-A, we propose a one step full pricing model, which includes routing and wavelength assignment inside the same mathematical formulation. In Section IV-B, we present a simplified version of pricing model, which leads to a routing only configuration as a first step. Then, we apply graph coloring on the obtained configuration to assign wavelengths to requests, and finally we adjust the reduce cost calculation accordingly. In Section IV-C, we provide a detailed description on how to coordinate these two models.

A. FSN Pricing Problem

This pricing problem aims to find new "improving" FSN configurations, i.e., a configuration such that if added to the current Restricted Master Problem (RMP) will improve the value of its linear relaxation, where an FSN configuration is defined by the FSN backbone (undirected tree supporting the FSN), its set of links supporting the routing of the requests, and the wavelength assignment. Variables. Note that several variables of the pricing problem correspond to parameters of the Master Problem. In order to alleviate the notations, and with a slight abuse of notations, we denote them the same way.

• a = 1 if link is in the FSN under construction, 0 otherwise, for all ∈ L • a v = 1 if node v belongs to the FSN under construction, 0 otherwise, for all v ∈ V • α e = 1 if edge e belongs to the backbone of the selected tree under construction, 0 otherwise, for all e ∈ E. • x k = 1 if request k is routed on the FSN under construction, 0 otherwise, for all k ∈ K • ϕ k = 1 if the routing of request k goes through link , or if the channel used in the routing of request k propagates on link because it is not filtered, 0 otherwise, for all ∈ L, k ∈ K • ψ k = 1 if the routing of request k goes through link between its source and destination, 0 otherwise, for all ∈ L, k ∈ K • θ kk = 1 if request k and request k are in conflict on the FSN under construction, 0 otherwise.

• w kk = ψ k ϕ k , for all ∈ L, k, k ∈ K.
In other words, with w kk = 1 identifying a link on which k and k are conflicting either between their source and destination nodes, or with one of them being routed to the broadcast of the other request, and 0 otherwise (no conflict). • β λ = 1 if the wavelength λ is used on the FSN under construction, 0 otherwise, for all λ ∈ Λ • β λk = 1 if request k is assigned the wavelength λ on the FSN under construction, 0 otherwise, for all k ∈ K, λ ∈ Λ The objective of the pricing problem with respect to the z F variables (or columns made up of their coefficients in the constraint matrix of the master problem) corresponds to the reduced cost of these variables, following the Dantzig-Wolfe decomposition algorithm which equals the simplex algorithm with an implicit enumeration of variables. Therefore, variables of the pricing problem are weighted with the values of dual variables of the master problem, and we refer the reader who is not familiar with (delayed) linear programming (aka Dantzig-Wolfe algorithm) to the seminal book of Chvátal [START_REF]Linear Programming[END_REF]. Objective:

min 0 + λ∈Λ u (2) λ β λ - k∈K u (3) k x k + ∈L u (4) a + u (5) (8)
subject to:

Construct an undirected tree

e={v,v }∈E: v,v ∈V α e ≤ |V | -1 V ⊆ V, |V | ≥ 3 (9)
v∈V

a v = e∈E α e + 1 ( 10
)
2α e ≤ a v + a v v, v ∈ V, e = {v, v } (11) e∈ω(v) α e ≥ a v v ∈ V (12) a ≤ α e , a ≤ α e = (v, v ), = (v , v) v, v ∈ V, e = {v, v } (13)
Routing of the requests

ϕ k ≤ a k ∈ K, ∈ L (14)
a ≤ k∈K ϕ k ∈ L (15)
ϕ k + ϕ k ≤ 1 = (v, v ), = (v , v) : v ∈ V, v ∈ V (16)
Flow constraints for broadcast stream

∈IN(d k ) ϕ k = ∈OUT(s k ) ϕ k = x k k ∈ K (17) ∈IN(s k ) ϕ k = 0 k ∈ K (18)
ϕ k ≤ ∈IN(v) ϕ k + 1 -a k ∈ K, v ∈ V \ {s k }, ∈ OUT(v) ( 19
)
ϕ k ≤ ϕ k + 2 -a -a k ∈ K, v ∈ V, ∈ IN(v), ∈ OUT(v) \ { } (20)
Flow constraints between source and destination

∈IN(d k ) ψ k = ∈OUT(s k ) ψ k = x k k ∈ K (21) ∈IN(v) ψ k = ∈OUT(v) ψ k ≤ x k k ∈ K, v ∈ V \ (s k , d k ) (22) ∈OUT(d k ) ψ k = ∈IN(s k ) ψ k = 0 k ∈ K (23) ψ k ≤ ϕ k k ∈ K, ∈ L (24) ∈L DIST ψ k ≤ REACH DIST k ∈ K (25)
Identify wavelength conflicting paths

θ kk ≥ ψ k + ψ k -1 ∈ L, k, k ∈ K (26) θ kk ≥ ψ k + ϕ k -1 ∈ L, k, k ∈ K (27) θ kk ≥ ψ k + ϕ k -1 ∈ L, k, k ∈ K (28) θ kk ≤ ∈L (w kk + w k k ) k, k ∈ K (29) w kk ≤ ψ k ∈ L, k, k ∈ K (30) w kk ≤ ϕ k ∈ L, k, k ∈ K (31) w kk ≥ ψ k + ϕ k -1 ∈ L, k, k ∈ K (32) Wavelength Assignment λ∈Λ β λk = x k k ∈ K (33)
β λk ≤ β λ λ ∈ Λ, k ∈ K (34) β λk + β λk ≤ 2 -θ kk λ ∈ Λ, k, k ∈ K (35)
Definition of variables

α e ∈ {0, 1} e = {v, v } ∈ E (36) a v ∈ {0, 1} v ∈ V (37) a ∈ {0, 1} ∈ L (38) ϕ k ∈ {0, 1} ∈ L, k ∈ K (39) ψ k ∈ {0, 1} ∈ L, k ∈ K (40) x k ∈ {0, 1} k ∈ K (41) θ kk ∈ {0, 1} k, k ∈ K (42) β λ ∈ {0, 1} λ ∈ Λ (43) β λk ∈ {0, 1} λ ∈ Λ, k ∈ K. (44) w kk ∈ {0, 1} k, k ∈ K, ∈ L. (45)
Constraints ( 9) -( 13) aim to restrict the undirected tree condition. Constraints (9) are the classical subtour elimination constraints in order to guarantee an acyclic graph structure (i.e., a supporting tree for the f-subnet under construction) [START_REF] Nemhauser | Integer and Combinatorial Optimization[END_REF].

As the number of these constraints is exponential, we will enumerate them only implicitly throughout a so-called lazy constraint procedure [START_REF] Aguayo | Solving the single and multiple asymmetric traveling salesmen problems by generating subtour elimination constraints from integer solutions[END_REF]. That is, we first solve the model with [START_REF] Ayoub | Filterless and semi-filterless solutions in a metro-HAUL network architecture[END_REF] restricted to a single constraint, e.g., V = V only. Then, if the set of selected edges is not connected (and so may contain a cycle), we add a constraint [START_REF] Ayoub | Filterless and semi-filterless solutions in a metro-HAUL network architecture[END_REF] for each subset of vertices V corresponding to the vertices of a connected component of the set of selected edges and solve again the model. We repeat this process until all constraints are satisfied. In practice, the number of introduced constraints is much smaller than the exponential number of constraints associated with [START_REF] Ayoub | Filterless and semi-filterless solutions in a metro-HAUL network architecture[END_REF].

Constraint [START_REF] Jaumard | Optimal design of filterless optical networks[END_REF] ensures that the selected edges correspond to the edges of a tree and so that a unique f-subnet will be generated in the entire network. Constraints [START_REF] Jaumard | A two sub-problem decomposition for the optimal design of filterless optical networks[END_REF] and [START_REF] Ayoub | Survivable virtual network mapping in filterless optical networks[END_REF] guarantee the consistency between node and edge variables: an edge is used in the undirected tree if and only if its two endpoints belong to it. Constraints [START_REF] Uzunidis | Bidirectional single-fiber filterless optical networks: modeling and experimental assessment[END_REF] enforce the consistency between link and edge variables: if a link is used, its associated edge belongs to the undirected tree. Here, the notation denotes the link in the opposite direction of link . Constraints ( 14) -( 16) take care of the routing of the requests, including the links hosting the unfiltered channels. Constraints ( 14) and ( 15) ensure that a request can be routed over link if and only if that link belongs to the f-subnet under construction and vice versa, if a link belongs to the output f-subnet structure, then at least one routing uses it. Constraints [START_REF] Nemhauser | Integer and Combinatorial Optimization[END_REF] prevent the use of both and in the routing of a given request, and so ensure that the signal can not be sent back to the sender. Constraints ( 17) -( 20) are the flow constraints for the requests, but include the broadcast effect. Constraints [START_REF] Aguayo | Solving the single and multiple asymmetric traveling salesmen problems by generating subtour elimination constraints from integer solutions[END_REF] states for each request, the total incoming traffic which associated to this request on the destination node should be 1 if this request is routed on this f-subnet, 0 otherwise, so as the total outgoing traffic for the source node of the request. Constraints [START_REF] Mehrotra | A column generation approach for graph coloring[END_REF] ensure that no flow supporting a request k can enter its source node s k .

The next two sets of Constraints ( 19) -( 20) are for taking care of the propagation of the unfiltered channels. Constraints [START_REF] Orlowski | SNDlib 1.0-survivable network design library[END_REF] enforce for every request k and node different from its source node (s k ), that if none of its incoming link is on the route or broadcast effect of request k, then none of the outgoing links can be used for either the routing or the broadcast effect of request k. Constraints (20) make sure on each node, for every request k, that any used outgoing links will carry the traffic for this request if there is incoming traffic of this request. This is to enforce that the broadcast effect will carry over to the traffic going forward. In other words, if ϕ k = a = 1, it will force ϕ k = 1 if a = 1, i.e., if link belongs to the FSN under construction.

Constraints (21) -( 22) are flow constraints that define the routing of request k between its source and destination. Constraints (23) impose no wavelength assignment on the incoming links of the source and on the outgoing links of the destination for the "filtered" routing of k.

Constraints (25) are the reach constraints. That is, for each request the routing distance between source and destination must not exceed a maximum distance (1,500 km). Constraints (24) define the relation between variables ϕ k and ψ k , which indicate ψ k is the subset of ϕ k .

Constraints ( 26) -( 28) are wavelength conflicts constraints expressing that: either k and k share a link between their source and destination nodes, see (26), or that one of the requests is routed to the broadcast part of the other request, see ( 27) and (28). Constraints (29) -(32) are the linearization of the constraints defining variables w kk , that is w kk = ψ k ϕ k . Constraints (33) express that one and only one wavelength is used per granted request, and that no wavelength is used if the request is not granted. Constraints (34) indicate that a wavelength is used in the FSN configuration if any of the granted request uses it. Constraints (35) enforce conflicting requests to use different wavelengths.

B. Relaxed FSN Pricing Problem

The pricing problem described in Section IV-A is a quite complex optimization problem. In order to speed up the solution process, we used a two-step process in the first iterations of column generation to solve the pricing problem.

The two step solution process consists in solving first a relaxation of the FSN pricing problem defined as the FSN pricing problem in which we omit the constraints related to the wavelength assignment, i.e., Constraints (26) to (35). We can then simplify the objective function of the pricing problem, and omit the term λ∈Λ u (2) λ β λ . This relaxation of the pricing problem outputs an FSN backbone, jointly with a routing for the granted requests. Once the relaxed FSN pricing problem is solved, we perform the wavelength assignment, which we reformulate as a graph coloring problem on the wavelength conflict graph. We then use the graph coloring algorithm of [START_REF] Mehrotra | A column generation approach for graph coloring[END_REF] and recompute the reduced cost in order to determine if we end up with an improving FSN configuration.

C. Detailed Solution Process

Figure 3(a) presents the detailed process of the two step solution of the FSN pricing problem. Figure 3(b) shows the overall solution process, using first the two step solution process of the FSN pricing problem, and then the one step process. The advantage of two step pricing model is to reduce the computational burden, especially when we solve data instances with larger networks. After the two step pricing model satisfies the optimality condition, the FSN pricing problem is solved at least once, in order to guarantee the optimal solution. Solution process steps: 1) Generate an initial solution: a) Compute an FSN supported by a spanning tree, which grants all the requests b) Compute the wavelength assignment on the FSN generated by step 1a (with, e.g., the ILP model of Mehrotra et al. [START_REF] Mehrotra | A column generation approach for graph coloring[END_REF]). Note that in the particular case of a single FSN, convergence is very fast, and therefore, we do not need to recourse to an initial solution to speed up the solution process. As we can see in Section V, in practice, either we directly generate an optimal spanning tree, or we find it after a first attempt. When searching for more than one FSN, recourse to an initial solution is needed in order to improve the convergence speed, and reduce the computational times. However, it is enough to generate an initial solution with a unique spanning tree.

V. NUMERICAL RESULTS

In this section, we present our computational results, which we compare both to those of the literature [START_REF] Tremblay | Passive filterless core networks based on advanced modulation and electrical compensation technologies[END_REF] and to our own previous results [START_REF] Jaumard | A two sub-problem decomposition for the optimal design of filterless optical networks[END_REF].

The software used to solve LP and ILP problems is the IBM ILOG CPLEX Optimization Studio version 12.8, and we use C++ to interact with its API. We run the programs on a x64 based processor, Intel(R) CPU E3-1270 v5 @ 3.6 GHz; RAM 32GB, with a Windows 10 operating system.

A. Data Sets

We used the same network topologies as Tremblay et al. [START_REF] Tremblay | Passive filterless core networks based on advanced modulation and electrical compensation technologies[END_REF], i.e., Italy, California, and Germany. We added two networks, also widely used in other papers on optical networking: Cost239 and USA, from [START_REF] Orlowski | SNDlib 1.0-survivable network design library[END_REF]. We recall their characteristics in Table I, i.e., the number of nodes and edges of the network topologies, as well as their number of requests. We consider a first set of requests, with a one unit wavelength requirement for each pair of nodes, as in [START_REF] Tremblay | Passive filterless core networks based on advanced modulation and electrical compensation technologies[END_REF]. We then added a second set of uniformly randomly generated requests on the pairs of nodes, with a random number of units in {0, 1, 2, 3, 4}.

B. Comparative computational results

We summarize our computational results on different number of FSNs in Tables II, III, and IV.

Each table contains the following results:

• z LP : the linear programming value -Model CG FOP: as we limited the number of solutions of the wavelength pricing problem to 500, the resulting LP value is the last computed one, and not necessarily an optimal one. It is denoted by zLP and is not a lower bound. -Model DW FOP: it is the optimal one, and is denoted by z LP .

• # col.: the number of generated configurations for each model

Step • W : an upper bound on the minimum number of required wavelengths in order to provision all requests • CPU: computational times for reaching the ε-optimal solution associated with the W value. 1) Single FSN case: For the single FSN case, we did not recourse to an initial solution, and solution of both models are done starting with an empty initial set of variables/columns.

The first striking result is the reduction of the computation time of 99% on average. We observe that all instances are solved exactly with the new DW FOP model as the lower and upper bounds are equal. This was not the case with the CG FOP model for which the accuracy of the solutions was ranging between 0% and 3.7%. While the number of selected configurations is one in the ILP solution, we observe that the number of generated columns is pretty large for the previous CG FOP model, which is likely to be the key explanatory factor for the huge differences for the computational results.

Results with the second set of requests also have the property that W = z LP . However, computational times are increasing significantly, especially for the last network (USA).

2) Two FSN case: In the case of two FSNs, we use an initial solution made of a single FSN. We generate a single FSN to satisfy all the requests with the objective to minimize the link usage. In addition to the comparison with the CG FOP model [START_REF] Jaumard | A two sub-problem decomposition for the optimal design of filterless optical networks[END_REF], we report some additional computational results of the literature, i.e., the numerical results of Tremblay et al. [START_REF] Tremblay | Passive filterless core networks based on advanced modulation and electrical compensation technologies[END_REF]. As for the single FSN case, we observe that the computation time needed for solving the new DW FOP model is significantly smaller than for the CG FOP model (average speedup of 55%). Moreover, the solutions obtained with the DW FOP model use less wavelengths than the solutions of both [START_REF] Jaumard | A two sub-problem decomposition for the optimal design of filterless optical networks[END_REF] and [START_REF] Tremblay | Passive filterless core networks based on advanced modulation and electrical compensation technologies[END_REF]. Another interesting observation is the reduction of the number of wavelengths when moving from a single FSN to two FSNs. More precisely, we reduce on average the number of wavelengths by 21, and the improvements range from 3 to 50 depending on the instance.

3) Three FSN case: In this experiment, we have changed Constraint [START_REF] Tzanakaki | Analysis of filterless wavelength converters employing cross-gain modulation in semiconductor optical amplifiers[END_REF] in order to select exactly n F = 3 FSNs instead of at most n F FSNs. With this modification, we observe that some instances are no longer feasible, i.e., for some instances it is not possible to partition the network and the set of requests into three FSNs while a solution with two FSNs is possible. Consequently, we analyse only the two instances for which a solution using three FSNs exists, namely California and Cost239.

As for the other experiments, we observe a significant reduction of the resolution time compared to the CG FOP model (92% faster on average). Furthermore, we obtain a solution for the Cost239 instance using only 17 wavelengths instead of 25, which is a significant improvement.

C. Impact of the number of FSNs

In this last experiment, we compare the solutions for the USA network, when moving from one to two FSNs.

Figure 4(a) represent the physical fiber connectivity of USA network. Figure 4(b) depicts the single spanning tree solution of USA network, and Figure 4(c) describes the solution with two FSNs.

In both Figure 4(b) and 4(c), on each link that is used in the network provisioning, we indicate different numbers. Bidirected links indicate that the edge is used is both directions, while uni-directed links mean the edge is used in only one direction. For each direction, the numbers define the number of used wavelengths on the plain links, while they indicate the number of wasted wavelengths (unfiltered channels) on the dash links.

We observe that even in the case of two FSNs, one of them has a spanning tree as a backbone tree.

VI. CONCLUSIONS

In this paper, we proposed a new mathematical model, DW FOP, for the design of filterless optical networks that significantly improves upon the mathematical model and algorithm proposed in [START_REF] Jaumard | A two sub-problem decomposition for the optimal design of filterless optical networks[END_REF], thanks to a decomposition into a single pricing problem combining the definition of a filterless optical sub-network and its full provisioning, including the wavelength assignment. The DW FOP model is much faster to solve and provides better solutions than previous proposals.

In the future, we plan to enhance further the model in order to handle multi-unit requests, without breaking them into single unit requests.

B.

  Jaumard has been supported by a Concordia University Research Chair (Tier I) on the Optimization of Communication Networks and by an NSERC (Natural Sciences and Engineering Research Council of Canada) grant. D. Coudert has been supported by the French National Research Agency (ANR), through the UCA JEDI Investments in the Future project with the reference number ANR-15-IDEX-0001, and the Inria associated-team project EfDyNet (Efficient Dynamic Resource Allocation in Networks). B. Jaumard is with the Department of Computer Science and Software Engineering, Concordia University, Montreal, QC, Canada. Y. Wang is with the Department of Computer Science and Software Engineering, Concordia University, Montreal, QC, Canada. D. Coudert is with Université Côte d'Azur, Inria, CNRS, I3S, Sophia Antipolis, France.

2 )

 2 Figure3(a) presents the detailed process of the two step solution of the FSN pricing problem.Figure 3(b) shows the overall solution process, using first the two step solution process of the FSN pricing problem, and then the one step process. The advantage of two step pricing model is to reduce the computational burden, especially when we solve data instances with larger networks. After the two step pricing model satisfies the optimality condition, the FSN pricing problem is solved at least once, in order to guarantee the optimal solution. Solution process steps: 1) Generate an initial solution: a) Compute an FSN supported by a spanning tree, which grants all the requests b) Compute the wavelength assignment on the FSN generated by step 1a (with, e.g., the ILP model of Mehrotra et al.[START_REF] Mehrotra | A column generation approach for graph coloring[END_REF]). 2) Apply the column generation algorithm in order to solve the linear relaxation of the model of DW FOP: a) Solve restricted master problem with current FSNs b) Solve the two step relaxed FSN pricing problem: i) Compute an FSN with the relaxed FSN pricing problem (see definition in Section IV-B) ii) Calculate the wavelength assignment on this FSN iii) Re-compute the reduced cost on the resulted provisioned FSN iv) If it generates an improving FSN configuration, add it to the RMP and return to Step 2a c) Solve one step full pricing problem with the mathematical model of IV-A. If it generates an improving wavelength configuration, add it to the RMP and return to Step 2a 3) Continuous relaxation of the Master Problem has been solved optimally 4) Solve the last generated restricted master problem with integer requirements for the variables, derive an ILP solution.

TABLE I :

 I Characteristics of network topologies and requests

TABLE III :

 III Comparative results with two FSNs Provisioned filterless sub-networks on two trees Tremblay et al. [8] CG FOP model (Jaumard et al. [11])

								DW FOP model	CPU
		W	z LP	# col.	W	CPU	z LP # col. W	CPU	ratio
	Italy	25	19.98 115/653 23	58m57s 20	340 21	54m26s 7.66 %
	California	-	117.4 23/952 122	1d15h37m 62.5 155 122 20h18m22s 48.74 %
	Germany17	88	62.3 29/1610 73 2d22h41m31s 61.5 106 73 1d8h53m44s 53.47 %
	Cost239	-	22.74 347/1183 28 4d22h17m37s 23	270 25 1d0h47m07s 79.05 %
	USA	-	41.24 117/1011 53	11h51m43s 31.25 302 53	1h40m05s 85.94 %
									average 54.97 %

TABLE IV :

 IV Network parameters for filterless solutions with three trees Provisioned filterless sub-networks on three trees Result of Tremblay et al. [8] CG FOP model (Jaumard et al. [11]) Result DW FOP model Result

									CPU
		W	zLP	# col.	W	CPU	z LP # col. W	CPU	ratio
	California	120	113.4 115/1040 120	5d10h27m29s 41.67 85 120 19h24m37s 85.12 %
	Cost239	-	15.7 357/1386 25	3d19h 9m34s 15.33 201 17	39m44s 99.27 %
								average 92.195 %