Soumaya Ben-Jabrallah

Florent Malloggi
email: orent.malloggi@cea.fr

Luc Belloni

Luc Girard

Dmitri Novikov

Cristian Mocuta

Dominique Thiaudière

Jean Daillant

Ben Soumaya

Jabrallah

†lions

Electrolytes at interfaces: accessing the first nanometers using X-ray standing waves

published or not. The documents may come

Introduction

The interfacial behavior of ions is of key importance in a number of phenomena and processes ranging from biological systems to microand nanouidics. Ion-surface interactions are important both in a direct way, e.g. in sorption or ion exchange related eects, or indirectly because modifying the eective surface charge they will aect the structure of the electric double layer, hence the interaction between charged particles in aqueous solutions.

1 In this context, particularly interesting are the so-called specic eects, referring to phenomena where ions of the same valency have a dierent eect. 8 The interfa- cial behavior of ions is also of key importance in micro-and nanouidics.

9 Phenomena occurring within the diuse Gouy-Chapman layer indeed become dominant in determining ow properties as channel size, Debye length and slip length become on the same order of magnitude, leading to unusual phenomena.

10 For example, with charged walls, channels can be lled with a unipolar solution of counter-ions, suggesting that the type and concentration of ions can be controlled by the surface charge density of the channel wall. 11 12 In addition, "specic eects" could allow for even ner control over the ion type as specic adsorption will eventually determine the double layer structure.

27,28

We show here that the method can be extended to light ions (down to chlorine or sulfur) on surfaces like silica which are of very broad interest in colloidal science.

Experimental

The x-ray standing wave technique.

30

In order to calculate the uorescence intensity, the sample is divided in as many layers j as necessary to describe the elemental distribution. We choose the layers to be thin enough for the local density of element k in layer j ρ k,j (z)

to be constant. The uorescence intensity for element k at grazing angle of incidence θ is then given by:

I k (θ) = B k s j=0 h j h j+1 ρ k,j (z)T k,j (z) | E j (z) | 2 dz, (1)
where the constant B k depends on the uores- cence yield, detection eciency and geometrical factors. E j (z) is the electric eld in layer j which is calculated in each layer using the matrix method. 31 T k,j (z) is the transmission of the uorescence intensity emitted by element k from layer j to the detector. As we are looking for nanometric layers close to the substrate, T k,j (z) is in fact independent of j. With this approximation,

I k (θ) = B k T k j h j h j+1 | E j (z) | 2 dz. (2)
Writing for a plane wave E j (x, z) = (A + j e ik z,j z + A - j e -ik z,j z )e +i(ωt-k x,j x) ,

where the A ± j coecients are determined using the matrix method following 31 and integrating over layer j, one simply has:

h j h j+1 | E j (z) | 2 dz = |A + j | 2 2k i z [1 -exp(-2k i z δh)] + |A - j | 2 2k i z [exp(2k i z δh) -1] + |A + j * A - j | 2ik r z [1 -exp(-2ik r z δh)] + |A + j A - j * | 2ik r z [exp(2ik r z δh) -1]. ( 4
)
with δh = h jh j+1 .

Though this method is accurate, it does not allow one to develop an intuition of the expected shape of the standing wave curve. For a monolayer located at z, one would have

I(θ, z) = 1 + R(θ) + 2 R(θ) cos(2πz/λ -φ(θ)) (5)
where R(θ) is the intensity of the reected beam and φ(θ) is the phase between the incident and reected beam which varies from π to 0 across the Bragg peak. Eq. ( 5) can help comparing dierent X-ray standing wave curves, even in slightly more complex cases. This is illustrated in Fig. 2 where we have represented the calculated uorescence from a hypothetical 2Å thick layer of ions adsorbed from a 1µm thick lm of solution. The uorescence intensity is rep-resented for dierent locations of the surface layer above the interface when this thin layer contains as many ions as the solution (50% of total number of ions), or only 5% of the total number of ions.

Obviously, the sensitivity is very high for the 50% case. However the dierent location can still be distinguished in a t with good quality data in the 5% case as the shape of the curve is still signicantly dierent. Sample thickness determination using X-ray reectivity.

Reectivity curves were recorded using a scintillator detector and carefully normalized to the incident beam intensity.

They were analyzed using the Reec software (A.Gibaud, G.Vignaud) based on the matrix method.

31

We used a model consisting of the multilayer substrate, a solution layer and the 4µm thick Ultralene R lm which was tted to the data, allowing to determine the thickness of the liquid layer. Obtaining good ts is important in order to constrain the standing waves analysis. The change in refractive index when a liquid layer is present modies the shape of the curve in the q z ≈ 0.1Å -1 region i.e. in between the critical angle for total external reection at the substrate-solution interface and the Bragg peak (Fig. 3). Though the solution layer is generally either too thick or not homogeneous enough to give rise to fringes that could be resolved in the reectivity curves, its thickness can nevertheless be determined from the shape of the reectivity curve in the region below and close to the critical angle (q z ∼ 0.05Å

-1 in Fig. 3).

For a very thick layer, the intensity at the critical angle can indeed be reduced by more than one order of magnitude compared to the empty cell because the incident and reected X-rays will have to travel hundreds of microns in the absorbing solution for grazing angles of incidence in the mrad range, even for a liquid layer thickness < 1µm. One should also notice that the thickness we determine is an average over the approximately 8 mm long and 250µm wide beam footprint on the sample.

X-ray uorescence measurements and analysis.

Fluorescence curves were recorded using a 4element energy sensitive detector and rst accurately normalized to the incident beam intensity. The uorescence curves were analyzed following the procedure detailed in Ref. Representative spectra are given on Fig. 4 for a mixture of 10mM solutions of KCl and KI. The Cl K α line at 2622 eV, the Ar K α line at 2957 eV (Ar is present in air), the K K α line at 3313 eV, the Ca K α and K β lines at 3690 eV and 4014 eV (Ca is an impurity in the Ultralene lm, not aecting the measurements) and the I L α line at 3938 eV and L β line at 4220 eV can in particular be identied on Fig. Standing wave curves.

Typical standing wave curves (uorescence intensity for each element as a function of the grazing angle of incidence across the Bragg peak) are given on Fig. 5.

Standing wave curves for the thicker sample (2.85 µm) have a prole which resembles the Bragg peak one. This is expected as these curves are dominated by the homogeneous ionic distribution in the bulk liquid (Fig. 2). Neglecting uorescence eects close to the interface for a thick solution layer, the uorescence will be excited by the incident and reected beam and will be roughly proportional to the incident beam plus reected beam intensity which is peaked at the Bragg angle. Dierences with respect to this shape give access to surface eects which can be determined and are more prominent for thinner samples (like the 250nm thick sample) also shown in Fig. 5. In such a case, the position and the shape of the standing wave curve is already indicative of the position of the ions with respect to the multilayer which is therefore unambiguously determined.

As mentionned above, we analyzed the standing wave curves using a very simple model consisting only of a homogeneous bulk solution and a Stern layer modelled as a stack of 2Å thick ionic layers, one for each ion. The sample is divided in as many layers as necessary to describe the elemental distribution and the electromagnetic eld is calculated in each layer using the matrix method.

31 Fitting the model (Eq. ( 5))

to the experimental curve by minimizing the χ 2 function, one obtains the composition, con-

Bragg S2 -I S2 -K S2 -Cl S1 -Cs S1 -K S1 -Cl
Grazing Angle of Incidence As can be seen in Fig. 6

d 2 φ dz 2 = - ρ 0 r . ( 6
)
In Eq. ( 6), ρ = i n i ez i is the charge density, with n i the concentration of the ion i of valency z i . e = 1.602×10 -19 C is the elementary charge, 0 = 8.85 × 10 -12 F/m is the vacuum permittivity and r the relative permittivity. Concen- trations in the double layer obey a Boltzmann distribution:

n i = n i∞ exp - z i eφ k B T , (7)
with k B the Boltzmann constant and T the temperature. Insertion of Eq.( 7) into Eq.( 6) gives after simple integration the Poisson-Boltzmann equation for the potential:

dφ dz 2 = 2k B T 0 r i n i∞ exp - z i eφ k B T -1 . (8)
The charge in the double layer can be written:

σ DL = ∞ d ρ(z)dz = - ∞ d 0 r d 2 φ dz 2 = -0 r dφ dz ∞ d , (9)
where d is the location of the outer Helmholtz plane which xes the border between the Stern layer and the diuse layer (Fig. 7). By inserting Eq. ( 8) in Eq. ( 9), one obtains.

- Si0 - Si0 - Si0 - Si0 - Si0 - Si0 - Si0 - Si0 + H + H + H + K + H + H + K + K + H + K - Cl - Cl - Cl + K + K - Cl z=d z=0
σ 2 DL = 2k B T 0 r i z i n i∞ exp - z i eφ(d) k B T -1 , (10)
known as the Grahame equation, which relates the charge in the electric double layer σ 2 DL to the potential at the outer Helmholtz plane φ(d).

On the other hand, the composition of the Stern layer is determined by surface equilibria, related to the dissociation of surface silanol groups for our silica surface. For protons and a monova-lent salt with cation M e + , one would have:

Γ SiO -[H + ] z=0 = Γ SiO -[H + ] exp - eφ(0) k B T = K H Γ SiOH , (11)
and

Γ SiO -[M e + ] z=d = Γ SiO -[M e + ] exp - eφ(d) k B T = K M e + Γ SiOC , (12)
where the concentrations involving silanol groups Γ SiOH , Γ SiO -and Γ SiOC are surface concentrations and it has been made explicit that H + ions are located at height 0 and M e + cations at height d. K H and K M e + are equilibrium constants. The total density of surface sites is

Γ = Γ SiO -+ Γ SiOH + Γ SiOC . (13)
The potential dependance in the Stern layer can be obtained by considering the Stern layer as a molecular condenser of capacitance C Stern where the potential drop is linear:

φ(d) -φ(0) = σ 0 /C Stern . (14)
We now have a full set of equations Eqs. ( 10)- [START_REF] Weber | Photoemission from Aqueous Alkali-Metal-Iodide Salt Solutions Using EUV Synchrotron Radiation[END_REF] which can be solved numerically to determine all surface and bulk concentrations.

A bare surface density of sites Γ = 8 SiO -sites /nm 2 (reduced by protonation and ion complexation, Eq. ( 13 In contrast to the previous estimates of pK K , our measurements provide a direct access to surface concentrations.

We We typically spread 500nl -1 µl droplets containing 3-6×10 14 ions for a mM solution. This crude calculation shows that a range of 0.1 to 10 is expected for the surface to bulk concentration ratio. The ratio of adsorbed ions in the Stern layer with respect to the total number of ions in the bulk is represented on Fig. 6 (bottom). This ratio ranges from 0.1 for thick and concentrated samples to ≈ 5 for thin layers and dilute solutions where a signicant amount of ions go to the surface in agreement with the previous estimate. It also implies that the above modeling using the Poisson-Boltzmann-Stern approach, which requires an innite reservoir, cannot be used without caution, at least for the least concentrated solutions (0.1mM).

A ratio of 0.1 to 10 falls well in the range where the dierent elements can be located as discussed earlier and explains why some curves are bulk-like with standing wave curves having a shape which resembles the Bragg peak (concentrated thick samples) and why surface-like features are more prominent in others.

It should be noted here that as the Stern layer is much thinner and concentrated than the total lm, the standing wave curve is mainly sensi- equal to the surface roughness. As we do not determine absolute concentrations, the distributions have been scaled arbitrarily. See Table 1 in Supporting Information for tting parameters.

Two representative examples are given on Fig. 9 for a 1mM KCl/CsCl mixture in a 50nm thick sample and a 1mM KCl/KI mixture in a 120nm thick sample. The XSW curves have been analysed using the model described above, where the ionic distributions have been convolved with the surface roughness. We rst note that there is a large concentration of both cations and anions in a 4 to 5 Å thick layer, with concentrations ≈ 100 times larger than in the bulk of the lm. Then, looking more deeply into the details, we see that in all cases, cations are located closer to the surface than anions, as expected for a negatively charged surface, and K + and Cs + show an almost equal propensity for the surface.

Regarding the competition between Cl and I -, we nd that the center of the Cl distribution is systematically closer to the surface, compared to the I distribution. In this case the shift is 2.5Å (see top insert of Fig. 9). This is in contrast with the air-water interface where I was shown to have a higher propen- sity for the surface compared to Cl -.

Conclusion

In this paper, we have shown that X-ray standing waves can be used to investigate ionic distributions in nanometer thick liquid layers using

24 TheseFigure 1 :

 241 Figure 1: (a) Schematics of the standing wave experiment. Incident (wave vector k in ) and reected (wave vector k r ) X-ray waves interfere to create a standing wave eld perpendicular to the interface. The species (sketched by blue dots) present in the XSW eld are excited and emit secondary uorescence photons for absorption edges below the incident beam energy. (b) Schematics of the experimental cell. The liquid layer is conned between the multilayer and a exible X-ray transparent membrane.

13 Determining

 13 the interfacial distribution of ions with sucient accuracy is, however, a challenge. Indeed, the ionic concentration deviates from the bulk concentration only in the Stern layer and diuse Gouy-Chapman layer: from 1nm to 100nm depending on the solution concentration. Several surface sensitive techniques like photoemission, 1417 sum-frequency generation (SFG), second-harmonic generation (SHG), 18,19 X-ray uorescence 7 or X-ray re- ectivity 20,21 have been used. More recently Atomic Force Microscopy has also been used to directly probe the electrical double layer. 2224 However, SFG and SHG have surface sensitivity but are not directly sensitive to ions and lack depth resolution, X-ray reectivity or neutron reectivity have no direct element sensitivity. Though photoemission has recently been extended to determine surface potential at the colloid/electrolyte interface, 16 it has no direct sensitivity to the ion distribution. In this study, we used the X-ray standing waves (XSW) technique which has the advantage of having both element and depth sensitivity to probe the rst nanometers at the solid surface. The technique rst developed to locate atoms in crystals or adsorbed at crystal surfaces 25 has also been applied to soft condensed matter. 26 It has in particular been used to investigate the adsorption of heavy ions on crystal surfaces from solutions.

  The X-ray standing waves (XSW) technique is based upon creating an electromagnetic standing wave eld by interference between incident and Bragg reected X-rays (Fig1a).29 The standing wave phase depends on the deviation from the exact Bragg angle and can be varied in a controlled way, thus moving the X-ray eld nodes in a range spanning half of the diraction plane spacing. The variation of secondary uorescence signal coming from the atoms, that get excited by this XSW eld, can be then used to locate the atom positions with (sub)Angstrom resolution. The method allows one to simultaneously investigate positions of dierent atoms, monitoring their specic emission lines. The standing waves are present where the incident and diracted waves overlap, which includes the solid/liquid interface and the bulk solution above. A specic XSW variant especially suited for investigation of liquid/solid interfaces employs XSW generators, based on articially grown multilayer structures (ML). Another advantage of this approach is the possibility to tune the resolution by changing the diracting plane spacing as well as the large choice of ML materials. In particular it would be possible to apply an external bias by embedding a conducting layer in the ML. The ML-generated XSW were used for instance for investigation of Zn 2+ distribution above a charged lipid layer 27 and in studies of electrochemical deposition of iodine on Pt.

  incidence [degree]Grazing angle of incidence[degree]

Figure 2 :

 2 Figure 2: Calculated uorescence curves for a 2Å thick layer in a 1µm liquid lm. In the bottom curves, the layer contains as many ions as the rest of the lm, and only 5% in the top curves. The layer is positioned at 10Å, 12Å, 14Å, 16Å, 18Å and 20Å above the reference plane. In the experiments, there is a top silica layer in between the reference plane and the layer and the interface is located ≈ 10-15 Å above the reference plane marked by a Cr layer in the experiments.

Figure 3 :

 3 Figure 3: Measured X-ray reectivity and best t for an empty cell (no liquid), a 4.7µm and a 200nm liquid layer thickness.

Fig. 3

 3 Fig.3shows experimental reectivity curves (as measured, without any corrections) as well as the corresponding ts before and after solu-

  software with a peak model based on the physics of detection, including a Voigtian peak shape, a tail and a shelf.

34 Figure 4 :

 344 Figure 4: Fluorescence intensity (experiments:symbols -ts:solid lines) recorded at 1.95 • and 2.1 • (Bragg peak) from an empty cell and a 250nm thick cell lled with a 10mM KCl + 10mM KI solution.

4 .

 4 An advantage of uorescence based methods is their high sensitivity to chemical elements in low quantities (traces). Samples with minute concentration (10 -4 M) or very thin liquid lay- ers (a few 10s of nms) could in particular be studied. Integrated over the beam footprint for a 100nm thick sample, the technique allows to resolve 1 femtomole or better. The dierences in the curves recorded at a grazing angle of incidence of 0.25 • below the Bragg peak (1.95 • ) and at the Bragg peak (graz- ing angle of incidence of 2.10 • ) in Fig. 4 directly result from the interfacial distribution of ions. Indeed, shifting the standing wave eld distribution, dierent positions above the interface are probed as shown in Fig. 2. This is however more clearly seen by representing the uorescence for each chemical element as a function of the grazing angle of incidence around the Bragg peak (standing wave curves).

Figure 5 :

 5 Figure 5: Standing wave curves (experiments:symbols -ts:solid lines) for a mixture of 1mM solutions of KCl and CsCl (S1) and 10mM solutions of KCl and KI (S2). Liquid layer sample thickness was 2.85 µm for (S1) and 250nm for (S2) respectively. Curves have been shiftedfor clarity. The Bragg peak is shown in red as reference. See Table1in Supporting Information for tting parameters.

Figure 6 :

 6 Figure 6: Top: Normalized uorescence of all ions as a function of the sample thickness times solution concentration. Bottom: Ratio of the number of ions in the Stern layer to the number of ions in the bulk as a function of product thickness × concentration.

Figure 7 :

 7 Figure 7: Schematics of the Stern layer and the double layer. z=d is the location of the outer Helmholtz plane.

  )), a capacitance of the Stern layer C Stern ≈ 2.9F/m 2 and K H = 10 -7.5 -10 -7 are well established in the literature. 36 The complexation constants for cations are less documented. Davis et al.

  37 report a much weaker specic adsorption for SiO 2 than for the other oxides, with pK K ≈ 0.5. Sonnefeld 38 suggests an even weaker binding. On the contrary, Tao et al.

  Fig. 8 bottom).

Figure 8 :

 8 Figure 8: Top: Γ SiOK as measured in this work and calculated using the Gouy-Chapman-Stern model as explained above with pK K =2. Bot- tom: Surface charge calculated using the same model with pK K =1,2 and 3 together with mea- surements from Dishon et al.

  consistency with our results.Partition between Stern layer and solution.

  tive to the Stern layer and it is more dicult to determine a concentration prole in the rest of the lm. Distribution within the Stern layer. The standing wave technique also allows for the determination of the position of ions relative to the multilayer lattice, and therefore relative to each other. Most interesting is of course the competition between K + and Cs + in KCl and CsCl mixtures and between Cl -

Figure 9 :

 9 Figure 9: X-ray standing wave curves and ionic distributions (insets) for a 1mM KCl/CsCl mixture (bottom) and a 1mM KCl/KI mixture (top). The KCl/CsCl sample was 50nm thick and the KCl/KI sample was 120nm thick sample. The width of the distributions has been set

  this higher propensity of the Cl ions for the surface is also reected in the Stern layer concentration as we nd a ≈ 1.6 ± 0.3 larger concentration of chloride ions compared to iodide ions. These two ndings are not necessarily in contradiction since both the surface chemistry (aecting short range interactions), including interaction with the other ions, and the dielectric function (aecting the dispersion forces) are dierent.

  surfaces and ions of broad interest in colloidal science and nanoscience. Not only the amount of ions in the bulk and Stern layer can be determined, but also their distribution in the direction perpendicular to the surface. Confronting such data to simulations taking into account all molecular interactions 40,41 should lead to a de- tailed understanding of ionic distributions at aqueous interfaces which is of utmost impor-tance for several areas of science.

  SiO -instead of the K + concentration in the Stern layer. Using the Gouy-Chapman-Stern model, we can calculate Γ SiO -

	10mM solution and signicantly less good for
	all concentrations with pK K =3.
	Previous experiments gave access to the eec-
	tive surface charge Γ using the same parameters for dierent concen-
	trations, allowing for an indirect comparison of
	our results to previous experiments	22 .	23 Note
	that in both Refs.	22 and,	23 surfaces were pre-
	pared by plasma cleaning which might lead to
	a dierent surface chemistry. Whereas Ref.	22
	shows a better agreement with pK K =3 Ref.,	23
	shows better agreement with pK K =2 as in our
	case.
	A dierence between the model and our mea-
	surements is that there is a nite reservoir
	of ions in our experiments, in contrast with
	the theory which assumes a xed concentra-
	tion. Using the Gouy-Chapman-Stern model,
	nd Potassium
	surface concentrations in the Stern layer of
	0.14/nm 2 , 0.3/nm 2 and 1.5/nm 2 for 0.1mM,
	1mM and 10mM KCl / KI mixtures respec-
	tively. As an indication, a 1/nm 2 surface con-
	centration would correspond for a 2Å thick
	layer to a concentration of 3mol/l.
	These results can be compared to the Gouy-
	Chapman-Stern model presented above in or-
	der to check whether they are in the expected
	range. With pK K = 2, we obtain Γ SiOK ≈ 0.13/nm 2 for a 0.1mM solution of KCl at pH=7,
	Γ SiOK ≈ 0.5/nm 2 for a 1mM solution and Γ SiOK ≈ 1.13/nm 2 for a 10mM solution, in fairly good agreement with our experimental re-
	sults (Fig. 8 (top)). As shown on Fig. 8 (top),
	the agreement is less good with pK K =1 for the

one expects to have in between ≈ 0.4 × 10 14 and ≈ 6 × 10 14 ions in the Stern layer on a 2 × 2cm 2 substrate (for Γ SiOK ≈ 0.1/nm 2 and Γ SiOK ≈ 1.5/nm 2 respectively).

Acknowledgements: The authors gratefully thank Christian Blot for manufacturing the sample cell and Elodie Barruet for concentration measurements of the solutions using capillary electrophoresis. The Marie Sklodowska Curie initial Training network SOMATAI under EU Grant Agreement No. 316866 is gratefully acknowledged for funding the PhD grant of Soumaya ben Jabrallah.


• pre-/post-conditions • imperative programs

• conclusion

• motivation

• algorithm Motivation 3

• to be fully published in articles or journals

• Coq proofs seem to me unreadable (and by normal human ?)

• algorithms on graphs = a good testbed

• nice algorithms should have simple formal proofs

• how to publish formal proofs ?

• Why3 allows mix of automatic and interactive proofs

• first-order logic is easy to understand A one-pass lineartime algorithm • depth-first search algorithm

• with pushing non visited vertices into a working stack

• and computing oldest vertex reachable by at most a single « back-edge »

• when that oldest vertex is equal to currently visited vertex, a new strongly connected component is in the working stack on top of current vertex.

• then pop working stack until currently visited vertex

The algorithm (2/2) The program (2/3)

• blacks, grays are sets of vertices; sccs is a set of sets of vertices • naming conventions:

x, y, z for vertices; b for black sets; s for stacks; cc for connected components;

Towards imperative program

• implementation of graphs

• vertices as integers in an array

• successors as lists for every vertex

• see http://jeanjacqueslevy.net/why3

Conclusion

• readable proofs ?

• simple algorithms should have simple proofs to be shown with a good formal precision

• compare with other proof systems

• further algorithms (in next talks ?)

• graphs represented with arrays + lists • topological sort, articulation points, sccK, sscT

• Why3 is a beautiful system but not so easy to use !

  version with lists and finite sets• the working stack is a list

  s2 ^y0 6 2 s2 ^edge x 0 y 0 • Coq proof: there exists x 0 , y 0 with y 0 is whitex 0 = x x 0 6 = ximpossible because no black to white impossible because successors are black y 0 2 s3 = stack x 0 = x x 0 6 = x impossible because crossedge impossible because m1  rank y' s1 < rank x s1 y 0 2 sccs impossible because sccs disjoint from stack Pre/Post-conditions (1/3)


Yosra Ben Slimen

Sylvain Allio

Julien Jacques

Model-Based Co-clustering for Functional Data

Keywords: co-clustering, functional data, SEM-Gibbs algorithm, latent block model, ICL-BIC criterion, mobile network, key performance indicators

In order to provide a simplified representation of key performance indicators for an easier analysis by mobile network maintainers, a model-based co-clustering algorithm for functional data is proposed. Co-clustering aims to identify block patterns in a data set from a simultaneous clustering of rows and columns. The algorithm relies on the latent block model in which each curve is identified by its functional principal components that are modeled by a multivariate Gaussian distribution whose parameters are block-specific. These latter are estimated by a stochastic EM algorithm embedding a Gibbs sampling. In order to select the numbers of row-and column-clusters, an ICL-BIC criterion is introduced. In addition to be the first co-clustering algorithm for functional data, the advantage of the proposed model is its ability to extract the hidden double structure induced by the data and its ability to deal with missing values. The model has proven its efficiency on simulated data and on a real data application that helps to optimize the topology of 4G mobile networks.

Introduction

With the introduction of new technologies and services in mobile networks, the complexity of these latter have increasingly grown creating an heterogeneous Email addresses: yosra.benslimen@orange.com (Yosra Ben Slimen), sylvain.allio@orange.com (Sylvain Allio), julien.jacques@univ-lyon2.fr (Julien Jacques)

Preprint submitted to Elsevier December 23, 2016 environment where different architectures (micro-, macro-, pico-, femto-cells)

and different radio access technologies (GSM, UMTS, LTE . . .) coexist. New difficulties have been resulted such as network management, optimization, troubleshooting and planning. Automated networks have also been created such as the self-organizing networks [START_REF] Ramiro | Self-Organizing Networks (SON): Self-Planning, Self-Optimization and Self-Healing for GSM, UMTS and LTE[END_REF] that demand new techniques for self-healing and self-optimization. Therefore, mobile operators need to deal with these new challenges in order to provide a top quality of services without increasing costs [START_REF] Mishra | Fundamentals of Cellular Network Planning and Optimisation[END_REF].

The quality of services is measured by data that are generated from multiple sources including key Performance Indicators (KPI) which are measurements collected from network elements such as transceivers, cells, or sites [START_REF] Mishra | Fundamentals of Cellular Network Planning and Optimisation[END_REF].

They are defined by mathematical formulas derived from different counters and computed periodically from the network with different temporal granularities (weekly, daily, hourly or less). For instance, Figure 1 illustrates a sample of p = 30 KPIs for n = 20 daily observations. From a statistical perspective, these KPIs are considered as functional data [START_REF] Ramsay | Functional data analysis[END_REF] which is a type of data that has recently appealed to researchers since they were for longtime inaccessible for statistics. However, with the advance of modern technology, more and more data are being recorded continuously during a time interval (or intermittently at numerous discrete time points). They become very frequent, not only in the telecommunication field, but in numerous other domains like medicine, economics and chemometrics (see [START_REF] Ramsay | Functional data analysis[END_REF] for an overview). Functional data is the observation (sample path) of a stochastic process X = {X(t), t ∈ T }, where T can be for instance a time interval, or any other continuous subset.

The KPIs may be specific to each radio access technology (GSM, UMTS, LTE,. . .) and to each constructor (Huawei, Ericsson,. . .). Therefore, as the number of technologies, services, cell types, and constructors grows, the KPIs observed by the support team become enormous and they may need to be analyzed over a large period (several weeks or months). Hence, on one hand, observing all the KPIs makes their daily analysis by the engineers a difficult task and also their treatment by the self-organizing networks more greedy in Co-clustering has successfully proven its efficiency in many applications such as recommendation systems [START_REF] Bennett | The netflix prize[END_REF] or text mining [START_REF] Dhillon | Co-clustering documents and words using bipartite spectral graph partitioning[END_REF]. According to [START_REF] Brault | Methods for co-clustering: a review[END_REF], two families of the block co-clustering techniques can be distinguished, namely: (a) the matrix reconstruction based family in which the problem is formulated as a matrix approximation using dissimilarity metrics and a set of constraints (see [START_REF] Banerjee | A generalized maximum entropy approach to bregman co-clustering and matrix approximation[END_REF] for an example) (b) the model-based family that uses probabilistic models in order to define the blocks [START_REF] Govaert | Clustering with block mixture models[END_REF][START_REF] Shan | Bayesian co-clustering[END_REF][START_REF] Wyse | Block clustering with collapsed latent block models[END_REF]. Many types of data have been treated when dealing with co-clustering such as categorical data [START_REF] Keribin | Estimation and selection for the latent block model on categorical data[END_REF], ordinal data [START_REF] Jacques | Model-based co-clustering for ordinal data[END_REF][START_REF] Matechou | Biclustering models for two-mode ordinal data[END_REF],

discrete [START_REF] Govaert | Latent Block Model for Contingency Table[END_REF] or continuous data [START_REF] Nadif | Model-Based Co-clustering for Continuous Data[END_REF]. However, there does not exist, to the best of our knowledge, co-clustering algorithm for functional data even though a lot of clustering algorithms have been proposed for this type of data (see [START_REF] Jacques | Functional data clustering: a survey[END_REF] for a survey).

The purpose of this paper is to propose a co-clustering technique based on the Latent Block Model (LBM, [START_REF] Govaert | Latent Block Model for Contingency Table[END_REF]) that we adapt for functional data. The LBM assumes local independence, i.e. the n×p curves are assumed to be independent once the row and column partitions are fixed. Since the notion of probability distribution for functional data is not well defined [START_REF] Delaigle | Defining probability density for a distribution of random functions[END_REF], a Functional Principal Components Analysis (FPCA, [START_REF] Ramsay | Functional data analysis[END_REF]) is used in order to plug the functional data into a finite-dimensional space. This strategy is frequent in functional data clustering and has proven its efficiency [START_REF] Bouveyron | Model-based clustering of time series in groupspecific functional subspaces[END_REF][START_REF] Chiou | Functional clustering and identifying substructures of longitudinal data[END_REF][START_REF] Jacques | Funclust: a curves clustering method using functional random variable density approximation[END_REF][START_REF] Jacques | Model-based clustering of multivariate functional data[END_REF]. Once each curve is being identified by its principal components, the probability distribution of these latter can be modeled by a multivariate (Gaussian) distribution with block-specific parameters.

The paper is organized as follows. Section 2 introduces the notations, the transformation of the discrete observations of curves into functional data and functional principal components analysis. The definition of the latent block model for functional data is described in Section 3. Given the numbers of rowclusters and column-clusters, a Stochastic EM algorithm embedding a Gibbs sampling (SEM-Gibbs) is proposed in Section 4 for the estimation of the model parameters. Since the number of row-and column-clusters must be estimated in practice, a strategy based on an ICL-BIC criterion is proposed in Section 5 that may be used to determine these numbers. The behavior of the model is studied on simulated data in Section 6. Finally, Section 7 presents an application of the co-clustering model on real data of mobile networks extracted within Orange Labs, France.

From discrete data to functional principal components

The data under study are a sample of n observations. Each observation is described by a set of p curves (functional features). The statistical model underlying data, represented by multivariate curves, is a stochastic process with continuous time:

X = {X(t)} t∈[0,T ] with X(t) = (X 1 (t), . . . , X p (t)) ∈ R p , p ≥ 2.
One way to explore such data is to perform a clustering of multivariate curves [START_REF] Jacques | Model-based clustering of multivariate functional data[END_REF]. In the present work, the goal is to go further by also clustering together the functional features having the same behavior.

As previously mentioned, the latent block model that we propose assumes a probability distribution on the functional principal components of the curves.

This section explains how to transform the discrete observations into functional data and then, how to perform a FPCA.

Transformation of the observed discrete curves

The main source of difficulty when dealing with functional data consists in the fact that these latter belong to an infinite-dimensional space, whereas in practice, data are generally observed at discrete time points and with some noise. Thus, in order to reflect the functional nature of data, a smoothing may be considered. Smoothing methods consider that the true curve belongs to a finite-dimensional space spanned by some basis of functions such as trigonometric functions, B-splines or wavelets (see [START_REF] Ramsay | Functional data analysis[END_REF] for a detailed study). Smoothing assumes that each observed curve x ij (1 ≤ i ≤ n, 1 ≤ j ≤ p) can be expressed as a linear combination of basis functions {φ j } =1,...,Mj :

x ij (t) = Mj =1 a ij φ j (t), t ∈ [0, T ], (1)
where {a ij } =1,...,Mj are the basis expansion coefficients. These coefficients can be estimated by least square smoothing for instance [START_REF] Ramsay | Functional data analysis[END_REF]. In this work, due to the nature of the KPIs under study, the same basis {φ } =1,...,M is used for all the functional features. The choice of the basis as well as the number of basis functions strongly depends to the nature of data. Hence, they can be set empirically.

Principal components analysis for functional data

From the set of functional data, it is interesting to have optimal representation of curves into a functional space of reduced dimension. The main tool to answer this request is the principal components analysis for functional data (FPCA, [START_REF] Ramsay | Functional data analysis[END_REF]). It consists in computing the principal components C h and principal factors f h of the Karhunen-Loeve expansion:

X(t) = µ(t) + h≥1 C h f h (t), t ∈ [0, T ]. (2)
When curves are assumed to be decomposed into a finite basis of function [START_REF] Ramiro | Self-Organizing Networks (SON): Self-Planning, Self-Optimization and Self-Healing for GSM, UMTS and LTE[END_REF], FPCA consists in a usual PCA of the basis expansion coefficients using a metric defined by the inner products between the basis functions. In theory, the number of principal components are infinite. However, in practice, due to the fact that the curves are observed at discrete time points and that they are approximated on a finite basis of functions, the maximum number of components one can compute is equal to the number M of basis functions used for approximation.

In this work, in order to project all the data onto the same FPCA space, functional principal components analysis is applied on the whole data set of curves

x, without distinction between curves from different observations or curves from different features. Moreover, in order to reduce the dimensionality of the problem, only the first m ≤ M principal components are considered. This number m is fixed empirically so that the principal components express a given part of the total variance.

Latent block model for functional data

This section describes the co-clustering technique for functional data based on the latent block model.

Let x = (x ij ) 1≤i≤n,1≤j≤p be the matrix of curves (x ij : x ij (t), t ∈ [0, T ])
whose rows are observations and whose columns are the functional features.

After functional principal components analysis, each curve x ij is summarized by its principal components

c ij = (c h ij ) 1≤h≤m . Let c = (c h ij ) 1≤i≤n,1≤j≤p,1≤h≤m
denotes the set of all the principal components.

The objective of co-clustering is to divide the data into K r row-clusters and K c column-clusters. The clusters are mixed in varying proportions denoted by α kr for the row-mixing proportion of the row-cluster k r and by β kc for the column-mixing proportion of the column-cluster k c . Hence, the data are summed up in K r × K c blocks. Each block contains data belonging to the same Gaussian distribution N (µ krkc , Σ krkc ).

Let v be the row-clustering matrix v = (v ikr ) i=1...n,kr=1...Kr with v ikr = 1 if row i belongs to the cluster k r , 0 otherwise. Let p(v ikr ) be the probability of row i to belong to the cluster k r , with the constraint that Kr kr=1 p(v ikr ) = 1. In the same way for the column partitioning, let w = (w jkr ) j=1...p,kc=1...Kc be the column-clustering matrix and let p(w jkc ) be the probability of column j to belong to the cluster k c . In the following, the straightforward ranges for i, j, h, k r and k c will be omitted for simplicity of notations.

The latent block model for functional data is defined by its density:

p(x; θ) = v∈V w∈W p(v; θ)p(w; θ)f (c|v, w; θ) (3)
where,

• V is the set of all possible partitions of the rows into K r groups, W is the set of all possible partitions of the columns into K c groups,

• p(v; θ) = ikr α v ikr kr ; p(w; θ) = jkc β w jkc kc , • f (c|v, w; θ) = ijkrkc p(c ij ; µ krkc , Σ krkc ) v ikr w jkc with:
p(•; µ krkc , Σ krkc ) is the m-variate Gaussian density with mean µ krkc = (µ hkrkc ) 1≤h≤m and covariance matrix Σ krkc ,

• θ = (α kr , β kc , µ krkc , Σ krkc ) 1≤kr≤Kr,1≤kc≤Kc , is the whole set of the model parameters.

The parameter θ of the latent block model is to be estimated, and we propose to use maximum likelihood inference.

Inference via a SEM-Gibbs algorithm

The objective of this section is to determine the latent blocks of the model by estimating the parameter θ using maximum likelihood inference. The maximum likelihood estimation of θ results in an optimization of the observed loglikelihood l(θ; x) = ln p(x; θ) where p(x; θ) is defined by [START_REF] Ramsay | Functional data analysis[END_REF]. The usual way used to maximize the log-likelihood in the presence of missing observations is the EM algorithm [START_REF] Dempster | Maximum likelihood from incomplete data via the EM algorithm[END_REF]. However, since the LBM involves a double missing structure, v and w, the maximum likelihood inference is computationally infeasible with an EM algorithm. For instance, with a data matrix of size 20 × 20 and with

K r = K c = 2, computing (3) requires K n r × K p c ≈ 10 12 terms.
For this work, we choose to use a stochastic version of the EM algorithm in which the missing data simulation is performed without the need of computing the whole missing data distribution thanks to a Gibbs sampler [START_REF] Keribin | Estimation d'un modèle à blocs latents par l'algorithme SEM[END_REF]. A second advantage of SEM-Gibbs algorithm is that it is expected to be insensitive to its initial values. Starting from an initial value of the parameter θ (0) and of the missing data w (0) , the q th iteration of the partial SEM-Gibbs alternates the following SE-Gibbs and M steps. SE-Gibbs step. Execute a small number (at least 1) of successive iterations of the two following steps:

1. Generate the row partition v (q+1) ikr |c, w (q) for all 1 ≤ i ≤ n, 1 ≤ k r ≤ K r according to: p(v ikr = 1|c, w (q) ; θ (q) ) = α (q) kr f kr (c i |w (q) ; θ (q) ) k r α (q) k r f k r (c i |w (q) ; θ (q) )
where

c i = (c h ij ) j,h and f kr (c i |w (q) ; θ (q) ) = jkc p(c ij ; µ (q) krkc , Σ (q)
krkc ) w (q) jkc .

Generate the column partition w

(q+1) jkc |c, v (q+1) for all 1 ≤ j ≤ p, 1 ≤ k c ≤ K c according to: p(w jkc = 1|c, v (q+1) ; θ (q) ) = β (q) kc f kc (c j |v (q+1) ; θ (q) ) k c β (q) k c f k c (c j |v (q+1) ; θ (q) )
where

c j = (c h ij ) i,h and f kc (c j |v (q+1) ; θ (q) ) = ikr p(c ij ; µ (q) krkc , Σ (q) krkc ) v (q+1) ikr . M step. Estimate θ (q+1) given v (q+1) , w (q+1) .
For the mixing proportions, the estimations are:

α (q+1) kr = 1 n i v (q+1) ikr and β (q+1) kc = 1 p j w (q+1) jkc ,
and for the means and the covariances of each block:

µ (q+1) krkc = 1 n (q+1) krkc i j c v (q+1) ikr w (q+1) jkc ij and Σ krkc (q+1) = 1 n (q+1) krkc S (q+1) krkc
where n

(q+1) krkc = ij v (q+1) ikr w (q+1) jkc and , S (q+1)
krkc = ij ((c ij -µ (q+1) krkc ) t (c ij -µ (q+1) krkc )) v (q+1) ikr w (q+1) jkc .
Choosing the parameters estimation and the final partition. The SE-Gibbs and M steps are iterated for a given number of iterations. After a burn-in period, the final estimation θ is defined by the mean of the sample distribution. For the final partitions, since v and w depend to each other, a new sampling of v, w| θ is simulated by successive SE-Gibbs steps with θ = θ. Then, every v ikr (respectively w jkc ) is obtained by computing the mode of the new sampling distribution.

Choice of the number of clusters

In the previous section, the numbers of clusters, K r in rows and K c in columns, are supposed to be known. However, in real applications, it may not be easy to precisely guess these numbers. Therefore, we propose to use the ICL-BIC criterion developed for continuous data in [START_REF] Lomet | Sélection de modèle pour la classification croisée de données continues[END_REF]:

ICL-BIC(K r , K c ) = log p(x, v, ŵ; θ) - K r -1 2 log n - K c -1 2 log p - ν 2 log(np)
where:

• v, ŵ and θ are the respective estimations of the row partition, the column partition and the model parameters obtained at the end of the estimation algorithm,

• The complete-log likelihood is given by: log p(x, v, ŵ; θ) = repeat

{K old r , K old c } = {K r , K c }; if ICL-BIC(K r + 1, K c ) ≥ ICL-BIC(K r , K c + 1) and K r < K max r then K r = K r + 1; else if ICL-BIC(K r + 1, K c ) < ICL-BIC(K r , K c + 1) and K c < K max c then K c = K c + 1; end end until (K r ≥ K max r and K c ≥ K max c ) or ICL-BIC(K r , K c ) ≤ ICL-BIC(K old
r , K old c ); Algorithm 1: Greedy search for exploring the numbers of clusters

Numerical experiments

The aim of this section is to evaluate, through simulation studies, the robustness of the co-clustering model and to determine its strengths and limitations.

The first experiment consists in producing a set of simulated data to test the behavior of SEM-Gibbs algorithm in terms of convergence, in order to choose the best number of iterations. The second experiment verifies the quality of the parameter estimations and of the final partitions resulted by SEM-Gibbs algorithm. The third experiment checks if the greedy search algorithm can detect the right number of row-and column-clusters.

Experimental setup

The experimental setup is composed of four steps. Let's assume that the data set is divided into 9 blocks (K r = K c = 3), the first step consists in assigning, to each block, a mean curve with discrete values. The 9 mean curves are chosen among the daily evolution of the KPIs under study, presented in the next section. They are also chosen according to two scenarios so that two families of experiences are considered with different levels of difficulty. The first family (F amily1) is set such that the co-clustering task is less difficult: only 2 couples among the 9 mean curves are close. The second family (F amily2) provides a more challenging situation: the mean curves of the blocks are set such that 4 couples among them are close.

In the second step, each mean curve is smoothed using B-splines as follows:

µ krkc (t) = M l=1 a krkcl φ l (t)
, where the number of basis functions is empirically set to M = 10. Figure 3 (respectively Figure 4) illustrates the F amily1 (respectively F amily2) and their corresponding smoothing.

The third step consists in simulating curves for each block. By using the basis expansion coefficients, each curve is simulated as follows : x(t) = l a l φ l (t) where a l ∼ N (a krkcl , 10). These simulations are used to generate a number s of data sets of size n × p, where n, p ∈ {50, 100, 500}.

In the fourth step, a FPCA is applied for each data set. The number of principal components m is chosen so that they cover at least 80% of the information. Figure 5

Tuning of the number of iterations

This experiment aims to choose the length of the SEM-Gibbs chain. It allows to verify that after a burn-in period, the simulations are stable which means that the SEM-Gibbs chain has achieved its stationary distribution. Consequently, it 265 allows to choose the number of iterations. The convergence is tested with q iterations of α

(q) kr , β (q)
kc and µ (q) krkc where q ∈ [0.

.100] over s = 20 simulations of data sets with different sizes (s simulations for each size). The number of blocks are set to K r = K c = 3 which are the right numbers of row-and columnclusters. Figure 6 illustrates the parameters convergence of all the simulations 270 where each curve represents the values of the corresponding parameter over q iterations and per simulation. The data sets are of size 50 × 50, 100 × 100 and 500 × 500, where the functional data are derived from the F amily1.

We notice that SEM-Gibbs algorithm has achieved its stationary distribution. Since by the 50 th iteration, the algorithm has already converged, for the rest of the experiment series, the number of iterations is set to q = 50 with a burn-in period equal to 20.

Validation of the SEM-Gibbs algorithm

In this second experiment, the purpose is to validate the parameters estimations resulted from SEM-Gibbs algorithm through s = 50 simulations of data sets with different sizes. Let K r = K c = 3 and the functional data are derived from both families. The comparison between the real parameters and their estimation is held in terms of the following metrics:

• ARI r (respectively ARI c ): the adjusted rand index of the real row-(respectively column-) clustering and its estimation,

• ∆α: the distance between the real α and its estimation, where ∆α = Kr kr =1 |α kr -αkr | Kr

, and similarly for ∆β,

• ∆µ: the distance between the real µ and its estimation, where ∆µ =

Kr kr =1 Kc kc =1 h l=1 |µ l kr kc -μ l kr kc | Krkcl .
Table 1 presents the results over the 50 simulations. Each metric is composed of two values: a mean distance (or ARI) followed by a standard deviation in parenthesis. As shown in Table 1, ARI r and ARI c are close to 1 which proves that the final partitions estimated by the model, v and ŵ, are very close to the real partitions v and w. α and β are two parameters that belong to [0, 1].

295

The fact that the metrics ∆α and ∆β are close to zero proves that these two parameters are successfully estimated. As for the parameter µ, the possible values depend to the data set. In our simulations, the range of values belongs to [-120, 100] for F amily1 and to [-120, 14000] for F amily2, and the data sets are not normalized. That is why, knowing the range of values, the results of ∆µ are acceptable enough to conclude that the estimation of the different parameters is close to the reality which proves the efficiency of our model. 2. We notice that the right number of clusters {3, 3} is selected most of the time even when the blocks have close behaviors (i.e. the functional data belong to F amily2). The results of both methods are equivalent but with the Greedy search algorithm, the execution time is faster, which proves the efficiency of the 315 Greedy search algorithm using the ICL-BIC criterion for choosing the number of blocks.

Application to mobile networks monitoring

This section presents an application of the proposed co-clustering algorithm to mobile network monitoring. Due to the huge number of KPIs captured contin-320 uously from the network, analyzing all of them daily is impossible for engineers.

As a result, they are forced to ignore most of the KPIs and to lose the information behind, which may affect the network performance and consequently the quality of services offered by the mobile operators. Even for the automated networks, the enormous amount of data treated by the algorithms may impose a 325 huge computational capacities and a time consuming tasks which is not helpful especially in the case of real-time applications.

In this part, we are interested in the study of the relationship between the behavior of different daily-captured KPIs for a specific geographical area. KPIs with similar behaviors will be grouped in the same block. In this way, the information induced by KPIs can be summarized which offers a simple representation for engineers as well as self-organizing networks.

Description of the real data

The data are collected using an internal tool of Orange Labs, France. We focus our study on Long Term Evolution (LTE) sites located in Lyon as illustrated in Figure 7. LTE is a standard technology that brings the cellular communication to the fourth generation (4G) era [START_REF] Rinne | Lte, the radio technology path towards 4g[END_REF]. In wireless telephony, a cell [START_REF] Mishra | Fundamentals of Cellular Network Planning and Optimisation[END_REF] is the geographical area covered by a cellular telephone transmitter.

The transmitter facility itself is called the cell site After pre-processing, the data set is composed of 1287 × 170 curves, each curve is identified by 8 principal components. After digging in the data of each block, it can be noticed that the proposed approach has succeed to discern the row-clusters in terms of days. As illustrated in Figure 10, the first row-cluster mainly contains, Monday, Tuesday, Wednesday and Friday. The second row-cluster mainly contains Thursday. The third rowcluster mainly contains Saturday and Sunday. It also contains in less proportions the week days. However, when digging into the cells of each row-cluster as described in Table 3, we notice that the third row-cluster exclusively contains the cells of two sites "Ile Roy" and "Fontaine saône" that happens to be isolated sites situated next to the "Saône" river of Lyon and to train railways. Having this information, we can conclude that these cells have atypical traffic due to their locations which makes the algorithm to consider their activities on work days similar to the activity of the other cells on week-ends. The fourth rowcluster has a similar behavior as the first row-cluster but it differs with the sites as presented in Table 3.

390

Regarding the column-clusters, the algorithm has dissociated the KPIs according to their families. As illustrated in Figure 10, we can notice that the The rest of the sites x x x x KPIs of type "Quality" are mainly in the first and second column-clusters.

The KPIs of type "Accessibility" and "Retainability" are mainly in the fourth column-cluster and the second column-cluster.

With this initial study on mobile networks, we notice that our model may help to dissociate the data in terms of cells and days which will help the mobile operators for a better planning of the topological structure of their networks.

The network topology [START_REF] Mishra | Fundamentals of Cellular Network Planning and Optimisation[END_REF] is the arrangement of the various elements of a network in a geographical area and the presented work may help with its optimization which is a very known problematic in the telecommunication field. However, in order to have more interesting results for self-healing and self-optimization, a deeper study is considered (particularly, with bigger K max r and K max c ).

Conclusion and Future works

While functional data analysis is widely used in many real applications, the co-clustering of functional data has never been proposed. In this paper, a modelbased co-clustering for functional data is introduced. In the presented latent block model, each block of curves is identified by the multivariate Gaussian distribution of the principal components of the curves. The model parameters are estimated using a SEM-Gibbs algorithm and the number of row-and columnclusters can be chosen by using a greedy search algorithm based on an ICL-BIC criterion. In addition, the model can also be used for clustering of multivari-ate functional data [START_REF] Jacques | Model-based clustering of multivariate functional data[END_REF], by simply setting the number of column clusters to p. Finally, through a simulation study, the proposed algorithm has proven its efficiency and the application of the model on mobile network monitoring has shown interesting results.

As future work, we will held a deeper study on mobile network troubleshooting, optimization and topology planning using the LBM for functional data.

Inspired by the clustering model presented in [START_REF] Jacques | Model-based clustering of multivariate functional data[END_REF], the proposed model can be improved by considering block-specific FPCA. Moreover, following the strategy of [START_REF] Banfield | Model-based Gaussian and non-Gaussian cluster-495 ing[END_REF][START_REF] Celeux | Gaussian parsimonious clustering models[END_REF], several parsimonious sub-models can be introduced by imposing constraints on the covariance matrix Σ krkc such as: (1) full covariance matrix common to all blocks; (2) diagonal block-specific covariance matrices; [START_REF] Ramsay | Functional data analysis[END_REF] diagonal covariance matrix common to all blocks. The study of the possible parsimonious sub-models is an interesting alternative in order to reduce the execution time. Finally, we are currently working on an R package that may be an interesting option for the prospective users of the model.

Figure 1 :

 1 Figure 1: An example of functional data set composed of 20 observations and 30 KPIs

Figure 2 :

 2 Figure 2: Binary data set (1), data re-organized by partition on I (2), by partitions on I and J simultaneously (3) and summary matrix (4)

2 .

 2 p(c ij ; μkrkc , Σkrkc ),• ν is the number of continuous parameters of the LBM:ν = K r K c m + m(m+1)Strategy of exploration. In the co-clustering context, exploring all the possiblecombinations of values of {K r , K c } (with K r ≤ K max r and K c ≤ K max c) becomes rapidly computationally demanding. Inspired by the strategy developed in[START_REF] Robert | Latent block model and model selection with application in pharmacovigilance[END_REF], a greedy search algorithm is proposed. It allows to only explore a relevant subspace of possible combinations of {K r , K c } (see Algorithm 1). At each step, the algorithm consists in computing the ICL-BIC criterion of the models obtained with one additional cluster, either in row or in column. The solution with the best ICL-BIC criterion is retained and the previous step is repeated until the ICL-BIC criterion does no longer increase. initialization: {K r , K c } = {2, 2};

  illustrates the mean of FPCA variance proportions after s = 50 simulations of data sets of size 100 × 100. For all the samples' sizes and the families of simulation, 80% of the information is summarized by the first 3 principal components. Consequently, for what follows, the number of principal components is set to m = 3. 00:00 05:45 11:30 17:15 23

Figure 3 : 13 00Figure 4 :Figure 5 :

 31345 Figure 3: Discrete mean curves belonging to F amily1 and their smoothing

6. 4 .

 4 Choice of K r and K c The aim of this section is to verify that the greedy search algorithm presented by Algorithm 1 with the ICL-BIC criterion can detect the right numbers of rowand column-clusters. Therefore, s = 50 simulations are generated over data sets of different sizes where the functional data are derived from F amily1 or from F amily2. Given K max r = K max c = 4, we compute the number of selections of each couple {K r , K c }. The results are compared to the exhaustive method in which all the possible combinations of the couples {K r , K c } are tested as 310 presented in Table

  . The cell provided by a cell site can be from one mile to twenty miles in diameter, depending on terrain and transmission power. The data are extracted from 99 cells, for 170 KPIs and for 13 days: seven days among them correspond to an atypical week in the telecommunication field, distinguished by the end of summer holidays (from August, 25 to August 31, 2016) and the other six days correspond to a typical work days (from September, 16 to September, 21, 2016).The KPIs are extracted with a granularity of 15 minutes (therefore, each daily KPI contains 96 values). The extracted KPIs belong to different families of indicators. The first family, labeled "Quality", refers to the network's ability to address its supported services regarding its characteristics. The second family, labeled "Accessibility", refers to the network's ability to meet with the users demands for accessing to the different services. The third family, labeled "Retainability", refers to the network's ability to maintain its users connections regardless the quality. In this extract, 85 KPIs belong to "Quality", 60 KPIs belong to "Accessibility" and 25 KPIs belong to "Retainability".

Figure 7 :

 7 Figure 7: The geographical area used for the data extraction

}> 4 and K max c > 4 .

 4 = {4, 4}, the greedy search algorithm suggests to dissociate the data set in 4 row-clusters and 4 column-clusters, since they maximize the ICL-BIC criterion. In this first analysis, we are restricted to {4, 4} for an 20 easier interpretation, but since the Greedy search algorithm chose the maximum values, a deeper analysis should be tackled with K max r By fixing K r = K c = 4, this experiment aims to apply the LBM for functional data on the data set. Let the number of iterations of the SEM-Gibbs algorithm be equal to 50 with a burn-in period equal to 20. Figures 8 and 9 illustrate the estimated parameters per block i.e. αkr , βkc and μkrkc when 1 ≤ k r ≤ 4 and 1 ≤ k c ≤ 4. Let's notice that the third column-cluster only contains 2% of the KPIs.

Figure 8 :

 8 Figure 8: Estimated mixing proportions: α (left) and β (right)

Figure 9 :Figure 10 :

 910 Figure 9: Estimated mean per block

Table 1 :

 1 Validation of the parameters and the final partitions estimated by SEM-Gibbs algo-

	rithm over data sets of different sizes
			F amily1			F amily2
	Size of the data sets	50 × 50 100 × 100 500 × 500	50 × 50	100 × 100 500 × 500
	ARIr	0.93	0.97	0.98	0.95	0.96	1
		(0.18)	(0.11)	(0.09)	(0.15)	(0.14)	(0)
	ARIc	0.92	0.92	0.95	0.9	0.99	0.97
		(0.19)	(0.19)	(0.15)	(0.2)	(0.08)	(0.11)
	∆α	0.028	0.01	0.008	0.02	0.02	0
		(0.07)	(0.04)	(0.04)	(0.06)	(0.06)	(0)
	∆β	0.03	0.03	0.02	0.04	0.004	0.01
		(0.07)	(0.07)	(0.07)	(0.08)	(0.032)	(0.05)
	∆µ	2.14	13.2	0.92	590.13	33	103.13
		(4.01)	(30.02)	(3.12)	(1248.48)	(229.52)	(412.38)

Table 2 :

 2 The results of the exhaustive method and the Greedy search algorithm over 50

	simulations of data sets of different sizes
		F amily1		F amily2
	Size of the data sets	50 × 50	100 × 100	50 × 50	100 × 100
	Kr /kc
	Exhaustive method

Table 3 :

 3 Row-clusters repartition in terms of sites

	row-cluster	No.1 No.2 No.3 No.4
	Brignais	x		x
	Ile Roy			x
	Fontaine Saône			x
	Dugusclin	x	x	x
	Fosses Ours	x	x	x
	Saint Gryphe	x	x	x