Tomoki Asaka
email: asaka.tomoki@e.mbox.nagoya-u.ac.jp

Katsunori Yoshimatsu
email: yoshimatsu@nagoya-u.jp

Kai Schneider
email: kai.schneider@univ-amu.fr

A WAVELET-BASED THREE-DIMENSIONAL CONVOLUTIONAL NEURAL NETWORK FOR SUPERRESOLUTION OF TURBULENT VORTICITY

Keywords: superresolution, wavelets, isotropic turbulence, vortex structures

We develop a wavelet-based three-dimensional convolutional neural network (WCNN3d) for superresolution of coarse-grained data of homogeneous isotropic turbulence. The turbulent flow data are computed by high resolution direct numerical simulation (DNS), while the coarse-grained data are obtained by applying a Gaussian filter to the DNS data. The CNNs are trained with the DNS data and the coarse-grained data. We compare vorticity-and velocity-based approaches and assess the proposed WCNN3d method in terms of flow visualization, enstrophy spectra and probability density functions. We show that orthogonal wavelets enhance the efficiency of the learning of CNN.

INTRODUCTION

In recent years, deep learning has been extensively applied to laminar and turbulent flows. For an overview we refer to review articles, e.g., Refs. [1,2,3]. Turbulent flows have many degrees of freedom, they are multi-scale, and have coherent vortices, e.g., vortex tubes, to which strong intermittency in space and time is attributed. Because direct numerical simulation (DNS) of turbulent flows needs high computational cost, many simulations using turbulence models, in which the small-scale motions are modeled, are performed. To get deep insight into the simulation data and to improve turbulence models, it is beneficial to reconstruct unresolved-scale data. To this end it is important to develop superresolution (SR) techniques which efficiently handle this task. One of the essential topics is hereby superresolution reconstruction of turbulent flow fields from low-pass filtered or coarse-grained turbulent fields, i.e., low-resolution (LR) turbulent fields. SR is originally an image processing technique reconstructing information that cannot be appropriately expressed by LR images. Dong et al. [START_REF] Dong | Computer Vision[END_REF] proposed a SR method using convolutional neural networks (CNN) with deep learning. CNN learns characteristic spatial structures by convolution of training data via compression and extraction of the data. Application of SR using two-dimensional (2D) CNN, so-called CNN2d, to turbulence has been proposed for 2D freely-decaying homogeneous turbulence [6], urban turbulence [5], three-dimensional (3D) forced homogeneous isotropic turbulence (e.g., Refs. [7,8]), and 3D turbulent channel flows (e.g., Refs. [7,8]). We here focus on 3D forced homogeneous isotropic turbulence which is one of the most canonical turbulent flows. In LR fields, small scales in the dissipation range are filtered out and are thus not available. Lui et al. [7] showed that convolutional neural networks using a time series of LR data as an input can indeed reconstruct statistics of velocity better than SR using a single time input data. Generative adversarial networks (GAN) for SR were introduced by Ledig et al. [START_REF] Ledig | IEEE Conference on Computer Vision and Pattern Recognition (CVPR)[END_REF]. The learning by GAN progresses such that the probability density function (PDF) of the generated data is well superimposed on that of the total data. Therefore, GAN can reconstruct the statistics of the total data. Kim et al. [8] showed that though only CNN2d does not reconstruct statistics of velocity and vorticity well, the SR using a kind of GAN, CycleGAN, well preserves these statistics. It is to be noted that their SR based on CycleGAN does not reconstruct strong vortex structures and position of vortices very well (see Fig. 5 in Ref. [8]).

Turbulence exhibits significant spatial and temporal intermittency in particular in the dissipation range. Wavelets, which are well-localized functions in both space and scale, provide a suitable tool to analyse spatially intermittent multi-scale fields [START_REF] Mallat | A Wavelet Tour of Signal Processing: The Sparse Way, Third Edition[END_REF]. Readers interested in the application of wavelets to turbulence can refer to review articles, e.g., Refs. [START_REF] Farge | [END_REF]12]. Coherent vorticity is extracted out of turbulence, by using orthonormal wavelets, which have a fast transform [13,14,15,16]. It was shown that the coherent vorticity well preserves the turbulence statistics, and consists of few percents of the degrees of freedom of turbulence. Furthermore, wavelet analysis is employed in deep learning. For example, Guo et al. [START_REF] Guo | IEEE Conference on Computer Vision and Pattern Recognition Workshops[END_REF] proposed an image SR method using CNN2d and orthonormal wavelets, which well reconstruct edges of structures.

In this paper, we propose a SR method based on the orthonormal wavelet transform and 3D convolutional neural networks (CNN3d) for 3D homogeneous isotropic turbulence. CNN3d is a neural network that learns local area of 3D structure via a 3D convolution filter. Taking the results [13,14,15,16], mentioned above, we conjecture that the use of wavelets enhances the efficiency of learning turbulent data in deep learning. We assess the proposed wavelet-based SR method by using flow visualization, PDFs of vorticity, and enstrophy spectra and compare the results obtained without using wavelets.

The remainder of the manuscript is organized as follows. In section 2 we describe the DNS approach and briefly recall wavelet analysis of turbulent flow. A wavelet-based SR method is presented in section 3. Numerical results are given in section 4 and finally conclusions are drawn in section 5.

DIRECT NUMERICAL SIMULATION AND WAVELET REPRESENTATION

In the following the governing equations with the DNS approach are briefly described and the orthogonal wavelet decomposition of a 3D vector field is presented.

GOVERNING EQUATIONS AND NUMERICAL SIMULATION

We consider 3D incompressible homogeneous isotropic turbulence in a periodic cube Ω = [0, 2π] 3 . The turbulent motion is governed by the Navier-Stokes equations

∂u u u ∂t + (u u u• • •∇)u u u = - 1 ρ ∇p + ν∇ 2 u u u + F F F, (1)
with the divergence-free condition ∇• • •u u u = 0. Here, u u u(x x x,t) is the velocity, p(x x x,t) is the pressure, F F F(x x x,t) is the external forcing, ρ is the constant density, ν is the constant kinematic viscosity, t is the time, The interval t = 1 ≈ 0.48T for the DNS at 512 3 grid points, where T is a large-eddy turnover time defined as T = L/u ′ , L is an integral length scale defined as L = π/(2u ′2 ) k max 0 E(k)/kdk at t = 9, E(k) is the energy spectrum, and u ′ = 2E/3, and E = 0.5 here. The resolution of the DNS are either n 3 g = 512 3 or 1024 3 where k max η ≈ 1.0 and 1.2, respectively. Here, η is the Kolmogorov microscale, and k max is the maximum wavenumber retained by the DNSs.

n 3 g (t, R λ )
training data 512 3 (10, 277), [START_REF] Farge | [END_REF]278) data 1 512 3 (9, 279) data 2 1024 3 (12.5, 373)

x x x = (x, y, z), and ∇ = (∂/∂x, ∂/∂y, ∂/∂z). The arguments x x x and t are omitted at will. The spatial average of u u u was set to zero. We carried out a series of DNSs with 512 3 grid points and 1024 3 grid points, using a Fourier spectral method and a fourth-order Runge-Kutta method for time integration. The aliasing errors are removed by a phase shift method and a spherical cutoff filter. The Fourier modes satisfying k < √ 2n g /3 are retained, where k = |k k k| is the modulus the wavevector, and n g is the number of grid points in each Cartesian direction. For the details on the DNS method, readers may refer to e.g., Ref. [START_REF] Kaneda | [END_REF]. The kinematic viscosity is set to ν = 2.8 × 10 -4 for n g = 512, while ν = 1.5 × 10 -4 for n g = 1024. The initial fields of the DNSs at t = 0 are divergence-free random velocity fields whose energy spectra ∝ k 2 exp(-k 2 /4), and the energy per unit volume, E, is set to 0.5. For forcing, negative viscosity was used as f f f at k < 2.5 such that the energy remains 0.5. We analyse the DNS data after a statistically stationary regime was obtained. Statistics of the different flow simulations are summarized in Table 1.

WAVELET ANALYSIS

Let v(x x x) be a component of a 3D vector field in a sub-domain 3 , where (L 0 , L 0 , L 0 ) is an appropriate position in the periodic box, 0 ≤ L 0 ≤ 2π, and v ∈ L 2 (Ω). The field in Ω s is sampled on 2 j L grid points in each Cartesian direction, where j L (≤ J) is a natural number, and J is the number of octaves in each space direction in the domain [0, 2π] (e.g., J = 9 for 512 grid points). It has to be noted that periodic boundary conditions are present at the boundary of Ω. For the sub-domains Ω s we impose mirror symmetry in each direction of the sub-domain. The wavelet functions ψ µ, j,i i i (x x x) yield an orthogonal basis at scale 2 -j , where µ(= 1, • • • , 7) denotes directions and i i i = (i x , i y , i z ) denotes the position.

Ω s = [L 0 , L 0 + L] 3 ⊆ Ω = [0, 2π]
The field v(x x x) can then be decomposed into an orthogonal wavelet series through the application of a periodization or folding technique [START_REF] Mallat | A Wavelet Tour of Signal Processing: The Sparse Way, Third Edition[END_REF]:

v(x x x) = ⟨v⟩ + J-1 ∑ j=0 v j (x x x), (2)
v j (x x x) = 7 ∑ µ=1 2 j -1 ∑ i 1 ,i 2 ,i 3 =0 ⟨v, ψ µ, j,i i i ⟩ψ µ, j,i i i (x x x), (3)
where ⟨•⟩ = Ω s •dx x x, and ⟨•, •⟩ denotes the L 2 -inner product defined as ⟨g 1 , g 2 ⟩ = Ω s g 1 (x x x) g 2 (x x x)dx x x. Note that ⟨v v v j ⟩ = 0 0 0. Using the scaling function φ j,i i i (x x x) based on φ(x x x), we obtain

j-1 ∑ m=0 v m (x x x) = 2 j -1 ∑ i 1 ,i 2 ,i 3 =0 ⟨v, φ j,i i i ⟩φ j,i i i . (4)
At scale 2 -j we have 7 × 2 3 j wavelet coefficients, ⟨v, ψ µ, j,i i i ⟩, and 2 3 j coefficients, ⟨v, φ j,i i i ⟩. Here we use Coiflet 12 wavelets which have filter length 12, 4 vanishing moments and compact support. In PyWavelets [17] which we use here they are denoted by 'coif2'.

A WAVELET BASED SUPERRESOLUTION METHOD

We develop a SR methodology for 3D homogeneous isotropic turbulence, using CNN3d and orthogonal wavelet analysis. We denote the methodology WCNN3d. In the following, we describe the training data and its preprocessing in section 3.1, and then present the procedure of WCNN3d in section 3.2. In section 3.3 we discuss the output data after the training.

TRAINING DATA

In the training of SR, we use two snapshots of x-components of high-resolution (HR) DNS data and low-pass filtered HR data (LF data). The HR data are the DNS data at 512 3 grid points described in Sec .2, while the LF data computed at 512 3 are obtained by application of an isotropic Gaussian filter to the HR DNS data. The low-pass filter is given by G

(k) = exp -π 2 k 2 /(24k 2 c )
, where we set k c to be 64 for the training data and the data 1 computed at 512 3 . For resolution 1024 3 , data 2 (cf. Table 1), we use k c = 128. We consider here three types of the training data, called type 1, type 2 and type 3: type 1 uses HR and LF data of ω x , the x-component of vorticity ω ω ω(= ∇ × u u u), while type 2 and type 3 use HR and LF data of u x , the x-component of velocity u u u.

In the type 1 of WCNN3d which learns HR and LF data from ω x , we divided the periodic box data with resolution 512 3 into 16 3 subcubes with 32 3 grid points, having side length L = π/8. Each subcube is decomposed into a 3D orthogonal wavelet series by using FWT only on one level. We then obtain wavelet coefficients in seven directions and scaling coefficients on the coarser level. For the analysis of the subcube data we used the option 'symm' for symmetric boundary conditions in PyWavelets. This option needs 5 extra elements of the array in each direction, and thus the number of the elements for the subcubes are finally 21 3 . For simplicity type 1 is denoted here by 'WCNN3d(ω)'.

Type 2 learning uses HR and LF data of u x and we first decompose u x given in a periodic box into 3D orthogonal wavelet series using the option 'periodization'. After multiplying with the weight 2 j at scale 2 -j for all j, we construct a weighted field by using the inverse FWT (IFWT). Then, we divide the weighted velocity into the subcubes at 32 3 grid points, and then follow the procedure of type 1. Type 2 is denoted by 'w-WCNN3d(u)'. The weight 2 j comes from a representative wavenumber of wavelets at scale 2 -j and is motivated by the norm equivalence of wavelets [START_REF] Mallat | A Wavelet Tour of Signal Processing: The Sparse Way, Third Edition[END_REF]. Weighting the wavelet coefficients of the velocity with 2 j yields gradient information of the velocity and is thus comparable with the vorticity-based type 1 approach, which only uses the anti-symmetric part of the velocity gradient.

To compare with w-WCNN3d(u), we also consider type 3 learning HR and LF data of u x where in contrast to type 2 no weight is imposed. The remaining procedure follows the one of type 1. Type 3 is denoted here by 'WCNN3d(u)'.

MACHINE LEARNING MODEL

We illustrate the procedure of WCNN3d in Fig. 1. The training scheme is implemented using the TensorFlow open-source library together with Python 2.7.5 interaction interface. We use Conv3D layer (keras.io) for 3D CNN. Our deep learning model has six layers; one input layer, four middle layers and one output layer. In these layers, we use 3D kernels with their sizes being 3 × 3 × 3 in order to catch local features of the turbulence data. The input and the middle layers consist of thirty-two filters of 3 × 3 × 3 3D kernels. The output layer consists of eight filters of 3 × 3 × 3 kernels. Zero-padding is employed at each layer, to keep the size of the output subcubes' data being the same as those of the input subcubes' data. As an activation function, we use a leaky rectified linear unit (ReLU) defined as

f (x) = x (x ≥ 0), αx (x < 0), (5)
where α = 0.2. The error was measured by the use of the of the difference between the output and correct data. The L 2 -regularization was used to suppress overfitting. Here, the Adam optimizer [20] is used. Our batch sizes are 32 and the number of epochs is set to 30.

To show the influence of the use of the wavelets, we also trained a CNN3d without using wavelets. In this CNN3d, the training is performed in physical space.

TEST DATA AND OUTPUT DATA

After the training was completed, we input the LF data of the DNS data at another time instant (another snapshot), corresponding to data 1 and data 2 in Table 1, to the model described in section 3.2. It has to be noted that we use x, y and z components of either the LF vorticity or the LF velocity, which depends on what type of the learning model we use as input data. In type 1, WCNN3d(ω), and CNN3d without wavelets, we use ω ω ω, while in type 2, w-WCNN3d(u), and type 3, WCNN3d(u), we use u u u. Then, our deep learning model predicts the eight sets of the coefficients, i.e., 7 directions of wavelet coefficients and the scaling coefficient. Using IFWT, we obtain 16 3 subcubes with 32 3 grid points in physical space. All of the subcubes are subsequently merged into the data yielding a flow field with 512 3 grid points. In w-CNN3d(u), after the merging, we apply the FWT, divide the resulting wavelet coefficients by 2 j for all j, and then use IFWT to get the output data. Note that no thresholding is applied in the use of wavelet decomposition. Moreover, we use the DNS data at higher Re λ with 1024 3 grid points, cf. data 2 in Table 1. We input the LF data of data 2 where k c = 128 to our models, WCNN(ω) and w-WCNN(u). The subcube size is 32 3 which is the same as what our model used for data 1. For the high resolution case we show here only the results using w-WCNN(u). WCNN(ω) shows better results than w-WCNN(u).

Tomoki Asaka, Katsunori Yoshimatsu and Kai Schneider Table 2: Mean values m and standard deviation σ for the total data, the LF data, the field predicted by WCNN3d(ω), and the field predicted by w-WCNN3d(u). m σ total 12.8 10.4 LR 8.9 6.37 WCNN3d(ω) 12.5 10.0 w-WCNN3d(u) 12.3 9.66 Figures 2(a),(b),(c) and (d) respectively show the visualization for the total data, i.e., the HR field of data 1, the LF field of data 1, the field predicted by WCNN3d(ω), and the field predicted by w-WCNN3d(u). The visualizations show isosurfaces satisfying |ω| = m + 3σ in a subcube with 256 3 grid points, where m and σ are the mean and the standard deviation of |ω|, respectively. The corresponding values are presented in Table 2. In Fig. 2(a), the total field, we can see well-pronounced vortex tubes. As expected the LF field, shown in Fig. 2(b), exhibits much less tubes, which is due to the application of the low-pass filter to the total field. For the predicted vorticity fields using either WCNN3d(ω) or w-WCNN3d(u), it can be observed in Figs. 2(c) and(d) that the visualizations of the reconstructed fields are in good agreement with those of the total field. These observations are confirmed by 2D cuts of ω x in the plane z = 0, shown in Fig. 3. The effect of the low-pass filtering is nicely seen in Fig. 3(b), where the scale of the coherent structures seems larger than that in Fig. 3(a). In contrast, the structures in the predicted field, shown in Figs. 3(c) and (d), agree well with those in the total field, shown in Fig. 3(a). Hence we can conclude from the 2D and 3D visualizations in Figs. 2 and3 that the deep learning model recovers the structure of the total field. We remark that in our models, we divide the data in the periodic box into subcubes, as mentioned in section 3.2. However, in the visualizations shown in Figs. 2 and3 the influence of the division seems negligible.

To quantify the above observations statistically, we computed the PDFs of ω x and enstrophy spectra. In this paper, we focus on small-scale behavior which can be super-resolved, and hence we consider the statistics of ω ω ω instead of u u u. Figure 4 shows the PDFs of ω x and the enstrophy spectrum Z(k) for the for the total vorticity, the LF vorticity, the vorticity predicted by WCNN3d(ω) using data 1 with resolution 512 3 . For comparison we also performed SR of the LF data using CNN3d without wavelets. In Fig. 4(a), we observe that the PDF of the LF data is much narrower than the total one, which is expected and in accordance with the values of σ in Table 2. We can also see that the PDF of the vorticity reconstructed by WCNN3d(ω) well overlaps with the total PDF. For CNN3d without using wavelets, the PDF of ω x is somewhat skewed, though the PDF seems to overlap the total PDF well. From the enstrophy spectra in Fig. 4(b), we can deduce that WCNN3d(ω) well preserves Z(k) of the total field, though there is a small departure of the predicted spectrum from the total one for kη ≳ 0.7. Concerning Z(k) for CNN3d, it can be seen that the enstrophy spectrum has somewhat smaller values than the total one. It can be noted that the pile up of Z(k) for the total field around kη = 1 is due to the use of the Fourier pseudo-spectral method and its resolution. From the enstrophy spectra we can conclude that the use of the wavelets enhances the deep learning of turbulent flows at least for the case considered here.

The PDFs of ω x and the enstrophy spectra Z(k) for the cases w-WCNN3d(u) and WCNN3d(u) are presented in Fig. 5 together with the total field and LF. We find that w-WCNN3d(u) preserves the total Tomoki Asaka, Katsunori Yoshimatsu and Kai Schneider PDF and the total enstrophy spectrum much better than WCNN3d(u). This shows that the weight 2 j in w-WCNN3d(u) is essential to our SR using the LF velocity field. Small magnitude of spikes in Z(k) at high k range can be seen for w-WCNN3d(u) and WCNN3d(u). These spikes are probably due to the division of the total field into the subcubes. Strictly speaking, the output fields generated by our model are not divergence-free. We thus examine whether the divergence of ω ω ω predicted by WCNN3d(ω) and the divergence of u u u predicted by w-WCNN3d(u) are negligible or not. On the basis of the Helmholtz decomposition, these fields are decomposed into ω ω ω = ∇ × A A A ω + ∇θ ω and u u u = ∇ × A A A u + ∇θ u , where A A A ω and A A A u are vector potentials and θ ω and θ u are scalar potentials. The divergent parts are given by ∇θ ω and ∇θ u . Figure 6 shows the spectrum of ∇θ ω and the k 2 weighted spectrum of ∇θ u . The latter means that the spectrum of ∇θ u is multiplied by k 2 . We find that these spectra have values much smaller than the enstrophy spectrum Z(k) predicted by WCNN3d(ω) and also the spectrum k 2 E(k) predicted by w-WCNN3d(u). Here, E(k) is the energy spectrum predicted by w-WCNN3d(u). Therefore, the influence of the divergence is negligible.

In addition, we observe some spikes in the spectra of the divergence parts. We confirm that the influence of the subcube division is likewise negligible, as already mentioned above. Furthermore, we evaluate whether our model w-WCNN3d(u), which was obtained by using the training data with 512 3 cubes, can well super-resolve the LF data obtained by the application of the Gaussian filter to data 2 computed with resolution 1024 3 (see the end of section 3.3). In Fig. 7(a), we can observe that the PDF of ω x predicted by w-CNN3d(u) is somewhat narrower than that of the total field. However, the enstrophy spectrum is well reconstructed by the SR, as shown in Fig. 7(b). We also confirmed that the divergence issue is not crucial (figure omitted for brevity).

CONCLUSIONS

We have developed wavelet-based 3D CNN (WCNN3d) for SR of coarse-grained data of homogeneous isotropic turbulence. Different approaches have been proposed, vorticity-and velocity-based models using orthogonal wavelets. For the later we also investigated a version with weighted wavelet coefficients which yields velocity gradient information. As input data of the WCNN3d, we used DNS data Tomoki Asaka, Katsunori Yoshimatsu and Kai Schneider of isotropic turbulence in a periodic box and the coarse-grained data are obtained by the application of a Gaussian low-pass filter to the DNS data. We assessed the WCNN3d in terms of 3D visualization of vorticity, PDF of vorticity, and enstrophy spectra. We found that WCNN3d well reproduces vorticity statistics and the positions of the vortices from coarse-grained vorticity fields. For the vorticity-based approach, the use of wavelets enhances deep learning of turbulent flows considered here thanks to the sparsity of the wavelet representation which improves deep learning. For the velocity-based approach, we showed that weighting the wavelet coefficients of velocity, which yields velocity gradient information due to norm equivalence, improves the accuracy and yields results similar to the vorticity-based model. Furthermore, we assessed the divergence issue of the predicted fields and showed that its impact is negligible. We demonstrated the capability to predict a turbulent flow whose Reynolds number is higher than the flows used for the training.

In this paper, our emphasis is put on SR when the intensity of data in wavenumber space is reduced, especially in the high wavenumber range. It would be interesting whether the use of wavelets can improve performance of the SR for the data where the modes in a high wavenumber range are completely cut off. However, this is beyond the scope of this study and remains an issue to be pursued.

Figure 1 :

 1 Figure 1: Illustration of SR network structure.

Figure 2 :

 2 Figure 2: The visualization of the modulus of vorticity |ω| for the isosurface, m + 3σ, for subcubes at size 256 3 : (a) total data, (b) LF data, (c) field predicted by WCNN3d(ω), and (d) field predicted by w-WCNN3d(u).Here, m and σ are the mean value and the standard deviation of the corresponding vorticity modulus |ω ω ω|.

Figure 3 :

 3 Figure 3: Visualization of the x-component of vorticity ω x in a 2D cut at z = 0: (a) the total data, (b) the LF data, (c) the field predicted by WCNN3d(ω), and (d) the field predicted by w-WCNN3d(u).

Figure 4 :

 4 Figure 4: (a) PDFs of ω x and (b) enstrophy spectra Z(k) for the total data, the LF data, the field predicted by WCNN3d(ω) together with the field predicted by CNN3d without wavelets. The wavenumber k is normalized by η of the total field.

Figure 5 :

 5 Figure 5: (a) PDFs of ω x , and (b) enstrophy spectra Z(k) of the field predicted by w-WCNN3d(u), and the field predicted by WCNN3d(u), together with the total data and the LF data. The wavenumber k is normalized by η of the total field.

Figure 6 :

 6 Figure 6: The spectrum of the divergence part of the vorticity field predicted by WCNN3d(ω) (green), and the k 2 -weighted spectrum of the divergence part of the velocity field predicted by w-WCNN3d(u) (red). For comparison, Z(k) for WCNN3d(ω) (blue) and k 2 E(k) for w-WCNN3d(u) (orange) are also plotted. The wavenumber k is normalized by η of the total field.

Figure 7 :

 7 Figure 7: The results using the LF data of data 2 with 1024 3 grid points as input data. (a) PDFs of ω x , and (b) the enstrophy spectra Z(k) for the total data, the LF data, and the field predicted by w-WCNN3d(u). The wavenumber k is normalized by η of the total field of data 2.

Table 1 :

 1 Statistics of the DNS data at different time instants, where R λ is the Taylor microscale Reynolds number.

Acknowledgement The computations presented here were carried out on the supercomputer "Flow" at Information Technology Center, Nagoya University. KS acknowledges partial support by the French Federation for Magnetic Fusion Studies (FR-FCM) and the Eurofusion consortium, funded by the Euratom research and training programme 2021-2022 under grant agreement No 633053. The views and opinions expressed herein do not necessarily reflect those of the European Commission.


Stevan Stanovic
email: stevan.stanovic@ensicaen.fr

Benoit Gaüzère
email: benoit.gauzere@insa-rouen.fr

Luc Brun
email: luc.brun@ensicaen.fr

Maximal Independent Vertex Set applied to Graph Pooling

Keywords: Graph Neural Networks, Graph Pooling, Graph Classification, Maximal Independant Vertex Set

Convolutional neural networks (CNN) have enabled major advances in image classification through convolution and pooling. In particular, image pooling transforms a connected discrete grid into a reduced grid with the same connectivity and allows reduction functions to take into account all the pixels of an image. However, a pooling satisfying such properties does not exist for graphs. Indeed, some methods are based on a vertex selection step which induces an important loss of information. Other methods learn a fuzzy clustering of vertex sets which induces almost complete reduced graphs. We propose to overcome both problems using a new pooling method, named MIVSPool. This method is based on a selection of vertices called surviving vertices using a Maximal Independent Vertex Set (MIVS) and an assignment of the remaining vertices to the survivors. Consequently, our method does not discard any vertex information nor artificially increase the density of the graph. Experimental results show an increase in accuracy for graph classification on various standard datasets.

Introduction

Convolutional neural networks (CNN) have enabled major advances in image classification. An image can be defined as a connected discrete grid and this property enables to define efficient convolution and pooling operations. Nevertheless, social networks, molecules or traffic infrastructures are not represented by a grid but by graphs. Convolution and pooling have been adapted to these structured data into Graph Neural Networks (GNN) [START_REF] Scarselli | The graph neural network model[END_REF]. The adaptation of convolution to graphs can be performed by using learned aggregation functions which combine the value of each vertex with the ones of its neighborhood [START_REF] Hamilton | Inductive representation learning on large graphs[END_REF][START_REF] Kipf | Semi-supervised classification with graph convolutional networks[END_REF]. For graph pooling, the adaptation is mainly performed by selecting vertices [START_REF] Bianchi | Hierarchical representation learning in graph neural networks with node decimation pooling[END_REF][START_REF] Gao | Graph u-nets[END_REF][START_REF] Lee | Self-attention graph pooling[END_REF][START_REF] Zhang | Hierarchical multi-view graph pooling with structure learning[END_REF] or by learning a fuzzy clustering [START_REF] Ying | Hierarchical graph representation learning with differentiable pooling[END_REF]. Fig. 2: Proposed graph pooling. Each node is associated to a score (Section 3.2). Based on this score, a MIVS is computed from which a reduction matrix S is derived. Applying S to both feature and structure leads to a reduced graph G (l+1) .

In this paper, we detail graph convolution and pooling as well as the construction of the reduced graph in Section 2. Section 3 presents our method (Fig. 2) and explains its differences from other pooling methods. We present in Section 4 our experiments where we study different heuristics to select surviving vertices and the mean complexity of our pooling algorithm. We finally compare our method to other pooling strategies on standard datasets using an unified architecture (Fig. 1) in Section 4.

Related Work

GNNs are neural networks applied to graphs. Inspired by CNNs, they are like the latter based on convolution and pooling operations. Let us note that convolution operation should be permutation equivariant. Intuitively, this condition states that a convolution operation may be permuted with a permutation of the graph's nodes.

Convolution operations diffuse the information by a message passing mechanism which allows to learn a representation for each node by aggregating the information of their n-hop neighborhood. Optimized according to a particular task, the resulting nodes' features can be used as input of a node classifier or regressor. According to [START_REF] Balcilar | Analyzing the expressive power of graph neural networks in a spectral perspective[END_REF], current aggregation functions correspond mainly to low-pass filters.

Graph Pooling

Considering graph level tasks like graph classification, we need to compute a graph representation as a fixed sized vector by summing node's representation. This operation is called global pooling. Global pooling realizes an aggregation of all the vertices of a graph and summarizes it in a single vector. Such aggregation must be permutation invariant and is usually performed using basic operators like a sum, mean, maximum or more complex ones [START_REF] Zhang | An end-to-end deep learning architecture for graph classification[END_REF].

However, global pooling performed on all the vertices of a graph leads to sum up an unbounded amount of information (proportional to the graph's size) into a fixed size vector and thus potentially induces a large amount of information loss. Many authors [START_REF] Lee | Self-attention graph pooling[END_REF][START_REF] Ying | Hierarchical graph representation learning with differentiable pooling[END_REF] have proposed to decompose GNNs in several steps alternating convolution and pooling to reduce the graph size while averaging vertices values. Such methods are called hierarchical pooling. The choice of the number of vertices for the reduced graph can be fixed or adapted according to the original graph's size with a ratio. The hierarchical pooling, like the convolution operation should be permutation equivariant.

Notations. For any pooling layer l, we consider a graph G (l) = (V (l) , E (l) ) where V (l) and E (l) are respectively the set of vertices and the set of edges of the graph. Let n l = |V (l) |, we can define A (l) ∈ R n l ×n l the adjacency matrix associated to the graph G (l) where A (l) ij = 1 if it exists an edge between the vertices i and j, 0 otherwise. We also note X (l) ∈ R n l ×f l the feature matrix of the graph G (l) where f l is the dimension of nodes' attributes.

Construction of the set of surviving vertices. It exists multiple methods to reduce the size of a graph within the GNN framework. However, most of these methods lead to the construction of a reduction matrix S (l) ∈ R n l ×n l+1 where n l and n l+1 are respectively the sizes of the original and the reduced graph. This matrix is used to define the attributes and the adjacency matrix of the reduced graph. Each surviving vertex i contributing to its own cluster, we suppose that S (l) i,i = 1. Construction of attributes and adjacency matrix. The reduction matrix S (l) is the basis of the construction of the reduced graph. Based on S (l) , the following two equations allow this construction:

X (l+1) = S (l) X (l) (1)
This last equation defines the attribute of each surviving vertex v i as a weighted sum of the attributes of the vertices v j of G (l) such that S (l) ji = 0.

A (l+1) = S (l) A (l) S (l) (2)
Equation ( 2) can be rewritten as follows for any pair of surviving vertex (i, j):

(S (l) A (l) S (l) ) i,j = n l k,l A (l) k,l S (l) k,i S (l) l,j (3)
Two surviving vertices are therefore adjacent in the reduced graph if they are adjacent in the initial graph (A

(l) i,j = S (l) i,i = S (l)
j,j = 1). Moreover, surviving vertices i and j are adjacent in the reduced graph if it exists a pair of non-surviving adjacent vertices (k, l) assigned respectively to i and j (A

(l) k,l = S (l) k,i = S (l)
l,j = 1). Families of methods. Pooling methods can be divided in two families. First family consists of methods based on the selection of surviving vertices on a given criteria. This criteria can be the result of a combinatorial algorithm [START_REF] Bianchi | Hierarchical representation learning in graph neural networks with node decimation pooling[END_REF] or a learning step like in Top-k methods [START_REF] Gao | Graph u-nets[END_REF][START_REF] Lee | Self-attention graph pooling[END_REF]. The second family regroups methods based on a node's clustering as in DiffPool [START_REF] Ying | Hierarchical graph representation learning with differentiable pooling[END_REF]. Each cluster is associated to a surviving vertex.

Methods of the second group use a fixed number of clusters. Hence, learning S (l) does not allow these methods to take into account the variable size and topology of the graphs. Moreover, due to training, such methods produce dense matrices S (l) with few nonzero values. Equation [START_REF] Borgwardt | Protein function prediction via graph kernels[END_REF] shows that the matrix A (l+1) is then dense and the corresponding graph has a density close to 1 (i.e. will be a complete graph or almost complete graph). Consequently, the structure of the graph is not respected. We say that S (l) is complete.

For Top-k methods, we define the reduction matrix S (l) as the restriction of the identity matrix I n l to the column indices idx corresponding to the surviving vertices:

S (l) = [I n l ] :,idx (4)
Given a surviving vertex i in the original graph, we have thus S

i,φ(i) = 1, where φ(i) denotes the column index of i in the reduced matrix S (l) . All other entries of S (l) are set to 0. The matrix S (l) is called selective since it selects the attributes of the surviving vertices and removes those of non-surviving vertices. Moreover, in this case, rows of S (l) corresponding to non-surviving vertices are equal to 0 and two surviving vertices will be adjacent if and only if they were adjacent before the reduction. This last point may induce the creation of disconnected reduced graphs. Moreover, the drop of non-survivors' features leads to an important loss of information. Let us note that MVPool [START_REF] Zhang | Hierarchical multi-view graph pooling with structure learning[END_REF] increases the density of the graph by considering power 2 or 3 of the adjacency matrix in order to limit the disconnections of the reduced graph. It additionally adds edges to the reduced graph through an additional layer called Structure Learning.

An alternative solution consists to drop Equation ( 2) and to perform a Kron reduction [START_REF] Bianchi | Hierarchical representation learning in graph neural networks with node decimation pooling[END_REF] in order to connect all pairs of surviving vertices adjacent to a same removed vertex. This reduction increases the density of the graph and a sparsification step is required. This last point can creates a disconnected reduced graph. Moreover, the time complexity of the Kron reduction is approximately O(( n l 2 ) 3 ), due to the inversion of a part of the Laplacian matrix of G (l) .

Let us finally note that it exists a graph pooling method which uses an approximation of a Maximum Independent Vertex Set called MEWIS [START_REF] Nouranizadeh | Maximum entropy weighted independent set pooling for graph neural networks[END_REF]. The layer used to compute this approximation significantly increases the complexity of the whole GNN. Moreover, the approximated result provided by this layer does not guarantee that the resulting set of selected vertices is even maximal.

Unlike other methods, we propose to preserve the structure of the original graph as well as the attribute information. We satisfy these properties thanks to a selection of surviving vertices distributed equally on the graph and with an assignment of nonsurviving vertices to surviving vertices. We named our pooling method MIVSPool. 3 Proposed Method

Maximal Independent Vertex Set (MIVS)

Before describing the algorithm, we introduce the notion of Maximal Independent Set and show how this notion can be applied to graph vertices. Maximal Independent Set. Let X be a finite set and N a neighborhood function. A subset J of X is independent if:

∀(x, y) ∈ J 2 : x / ∈ N (y) (5)
J is a subset of X such that for any (x, y) ∈ J 2 , x and y are not neighbors. The elements of J are called the surviving elements.

An independent set J is said to be maximal when no element can be added to it without breaking the independence property, i.e., when we have:

∀x ∈ X -J , ∃y ∈ J : x ∈ N (y) (6)
Equation ( 6) states that each non-surviving element has to be in the neighborhood of at least one element of J . The elements of J are denoted as survivors.

Using Equation ( 5) and ( 6), we have a Maximal Independent Set. We note that it is a maximal but not necessarily a maximum. Indeed, our Maximal Independent Set J is not necessarily those whose cardinality is maximum with respect to all Maximal Independent Sets of X .

If we interpret the construction of a Maximal Independent Set as a subsampling operation, Equation (5) can be interpreted as a condition preventing the oversampling (two adjacent vertices cannot be simultaneously selected) which thus guarantees an uniform distribution of survivors. Conversely, Equation (6) prevents subsampling: Any non-selected vertex is at a distance 1 from a surviving vertex.

Maximal Independent Vertex Set. A Maximal Independent Vertex Set (MIVS) [START_REF] Meer | Stochastic image pyramids[END_REF] of a graph G (l) = (V (l) , E (l) ) is therefore a Maximal Independent Set where the neighborhood is deduced from the edge set E (l) . Adapting Equation ( 5) and ( 6), we select surviving vertices V (l+1) as described below:

∀(v, v ) ∈ (V (l+1) ) 2 : (v, v ) / ∈ E (l) (7) ∀v ∈ V (l) -V (l+1) , ∃v ∈ V (l+1) : (v, v ) ∈ E (l) (8)
Equations ( 7) and ( 8) define the MIVS procedure and state that two surviving vertices cannot be neighbors and a non-surviving vertex must have at least one survivor in its neighborhood. Nevertheless, these two equations don't explain how to select these vertices.

Meer's Algorithm. A simple method has been proposed by Meer [START_REF] Meer | Stochastic image pyramids[END_REF] using an iterative procedure assigning a uniform random variable U([0, 1]) to each vertex. In this procedure, each surviving vertex corresponds at a local maximum of the random variable with respect to its neighbors. According to Equation [START_REF] Kipf | Semi-supervised classification with graph convolutional networks[END_REF], vertices adjacent to these survivors are labeled as non-survivors. Other vertices (not yet labeled) are labeled candidates and the algorithm iterates the selection of surviving and non-surviving vertices on this reduced vertex set (Fig. 3). The algorithm convergence is guaranteed since, at each iteration, at least one candidate is labeled as a survivor. For each vertex v i , its random variable is denoted by x i while Booleans p i and q i are T rue if v i is respectively a survivor and a candidate. The stopping criterion of the algorithm is obtained when all q i are F alse, i.e. when Equation ( 8) is satisfied. At the initialisation of the algorithm, all p i are F alse and all q i are T rue and, for each iteration k, variables p (k) i and q (k) i are updated by:

p (k+1) i = p (k) i ∨ (q (k) i ∧ (x i = max{x j |(v i , v j ) ∈ E (l) ∧ q (k) j })) q (k+1) i = ∧ j|(vi,vj )∈E (l) p (k+1) j (9)
In other words, a survivor at iteration k + 1 is a survivor at iteration k or is a candidate whose random variable (x i ) is greater than the ones of its candidate neighbors. A vertex is candidate at iteration k + 1 if it is not adjacent to a survivor. We note that the neighborhood includes the central vertex. This procedure only involves local computations and is therefore parallelizable.

Adaptation of MIVS to deep learning

Scoring system. Using Meer's algorithm [START_REF] Meer | Stochastic image pyramids[END_REF], the uniform random variable U([0, 1]) plays an important role in the selection of the surviving vertex set. We propose to modify this variable so that x i used in the algorithm is learnt and represents the relevance of vertex v i . We obtain this last property by using an attention mechanism like in SagPool [START_REF] Lee | Self-attention graph pooling[END_REF] where a score vector s ∈ R n l ×1 is returned by a GCN [START_REF] Kipf | Semi-supervised classification with graph convolutional networks[END_REF]. Using this score in our MIVS, we select a set of surviving vertices V (l+1) corresponding to local maxima of the function of interest encoded by s. An uniform distribution of vertices on the graph is guaranteed by the computation of the MIVS.

Assignment of Non-Surviving Vertices. In order to construct our reduction matrix S (l) and take all vertices information into account, we need to assign non-surviving vertices. As a reminder, Equation [START_REF] Lee | Self-attention graph pooling[END_REF] states that each non-survivor has at least a surviving neighbor. Assuming that the score function encodes the relevance of each vertex, we assign each non-surviving vertex to its surviving neighbor with the highest score. We obtain a reduction matrix S (l) with n l+1 clusters corresponding to surviving vertices:

∀v j ∈ V (l) -V (l+1) , ∃!v i ∈ V (l+1) |S (l) ji = 1 with i = argmax(s k , (v k , v j ) ∈ E (l) ∧ p N k ) ( 10
)
where N is the number of iterations required to have a MIVS and s the score vector.

Construction of reduced attributes and adjacency matrix. The construction of attributes of the reduced graph G (l+1) is obtained thanks to an average weighted by s from the set of aggregated vertices to the surviving vertices:

X (l+1) i = 1 j|S (l) ji =1 s j j|S (l) ji =1 s j X (l) j (11)
Equation ( 11) allows to take into account the importance of each vertex in the computation of the attributes of the reduced graph and put more attention on vertices with a high score. As a consequence the learnt vector s can be interpreted as a relevance value. This operation can be achieved by a transformation similar to Equation (1) by substituting S (l) by the matrix

S (l) = D 2 D 1 S (l) with D 1 = diag(s) and D 2 = diag( 1 1 D1S (l)
). The construction of the reduced adjacency matrix A (l+1) is obtained thanks to Equation [START_REF] Bianchi | Hierarchical representation learning in graph neural networks with node decimation pooling[END_REF].

Relaxation of MIVSPool. On some highly dense graphs, MIVSPool may provide a decimation ratio lower than 0.5. In order to correct this point, we additionally introduce MIVSPool comp. which is a MIVSPool where we force the addition of surviving vertices using a Top-k so that the pooling ratio remains always equal to 0.5. Note that Equation ( 7) is then no more valid while Equation (8) still holds.

Experiments 4.1 Datasets

To evaluate our MIVSPool, we test it on a benchmark of four standard datasets: D&D [START_REF] Dobson | Distinguishing enzyme structures from non-enzymes without alignments[END_REF], PROTEINS [START_REF] Borgwardt | Protein function prediction via graph kernels[END_REF][START_REF] Dobson | Distinguishing enzyme structures from non-enzymes without alignments[END_REF], NCI1 [START_REF] Wale | Comparison of descriptor spaces for chemical compound retrieval and classification[END_REF] and ENZYMES [START_REF] Borgwardt | Protein function prediction via graph kernels[END_REF]. The statistics of datasets are reported on Table 1. D&D and PROTEINS describe proteins and the aim is to classify them as enzyme or non-enzyme. Nodes represent the amino acids and two nodes are connected by an edge if they are less than 6 Ångström apart. NCI1 describes molecules and the purpose is to classify them as cancerous or non-cancerous. Each vertex stands for an atom and edges between vertices represent bonds between atoms. ENZYMES is a dataset of protein tertiary structures obtained from the BRENDA enzyme database.

Model Architecture and Training Procedure

The model architecture consists of K blocks made up of a GCN [START_REF] Kipf | Semi-supervised classification with graph convolutional networks[END_REF] convolution followed by a graph pooling. A Readout layer is applied after each block using a concatenation of the average and the maximum of vertices' features matrix X (l) . At the end of our network, the K Readout are summed and the result is sent to a Multilayer Perceptron. The latter is composed by three fully connected layers (respectively 256, 128 and 64 for the number of hidden neurons) and a droupout of 0.5 is applied between the first two. The classification is obtained by a Softmax layer (see Fig. 1).

For the training procedure, we use Pytorch Geometric and we evaluate our neural network with a 10-fold cross validation. We repeat this procedure ten times without setting the seed. The dataset is split in three parts: 80% for the training set, 10% for the validation set and 10% for the test set.

For the hyperparameters, we use the Adam optimizer, set the dimension of node representation, the batch size and the number of epochs respectively at 128, 512 and 1000 and an early stopping is applied if the validation loss did not improved after 100 epochs. A grid search is used for the learning rate within 1e -3 , 1e -4 , 1e -5 , the weight decay within 1e -3 , 1e -4 , 1e -5 and the number of blocks K in [START_REF] Borgwardt | Protein function prediction via graph kernels[END_REF][START_REF] Gao | Graph u-nets[END_REF] to find the best configuration.

Ablation Studies

Scoring function. encodes the relevance of each vertex and as a direct influence on the set of selected surviving vertices by MIVSPool. We vary it using a uniform random variable U([0, 1]) like in Meer [START_REF] Meer | Stochastic image pyramids[END_REF] ( MIVSPool rand ), a trainable normalized projection vector on the features X (l) like in [START_REF] Gao | Graph u-nets[END_REF] (MIVSPool T op-k ), a self-attention mechanism thanks to a graph convolution like in SagPool [START_REF] Lee | Self-attention graph pooling[END_REF] (MIVSPool SagP ool ) and finally a trainable multi-view system across structure and features graph information like in MVPool [START_REF] Zhang | Hierarchical multi-view graph pooling with structure learning[END_REF] (MIVSPool M V P ool ). For trainable variations, Equation [START_REF] Scarselli | The graph neural network model[END_REF] allows to propagate the training to the next convolution and, therefore, our pooling method is endto-end trainable. The results reported Table 2 show that the choice of the score function has a minor influence on the accuracy. Note that, since Top-k and SagPool trainable score functions don't consider structural information, they consequently perform poorly on ENZYMES. However, the use of the MVPool trainable score function allows to obtain the best accuracy on three datasets with comparable standard deviations to other heuristics. In the rest of our article, we choose to take MIVSPool M V P ool as a reference that we denote MIVSPool.

Average number of iterations required by MIVSPool for each pooling step. Table 3 presents the mean number of iterations for each dataset and each pooling step computed over 10 epochs. Note that the number of pooling steps is determined using K-folds for each dataset. Despite an important difference between graph sizes among datasets (Table 1), we note that the number of iterations is comparable between them and less than 5. Knowing that an iteration of MIVS on a graph G = (V, E) requires about |V|d max computations, where d max is the maximum degree of G, the computation of MIVS on one of the graphs of our three datasets is bounded by 5|V|d max . This complexity is less than the one needed by the Kron transformation (Section 2.1).

Comparison of MIVSPool according to other methods

We compare our method to three state-of-art methods: gPool [START_REF] Gao | Graph u-nets[END_REF], SagPool [START_REF] Lee | Self-attention graph pooling[END_REF] and MVPool-SL [START_REF] Zhang | Hierarchical multi-view graph pooling with structure learning[END_REF]. Results in Table 4 show that MIVSPool comp. and MIVSPool obtain the highest or second highest accuracies on D&D, PROTEINS and ENZYMES. This point confirms the efficiency of MIVSPool and the fact that some configurations (highly connected graphs) may induce low decimation rate by MIVSPool which slightly decreases the accuracy. For NCI1, the highest accuracy is obtained by MVPool-SL [START_REF] Zhang | Hierarchical multi-view graph pooling with structure learning[END_REF], MIVSPool comp. being ranked second. Let us note that the mean vertex's degree in this dataset is 2.15 ± 0.11 with many graphs being almost linear. On such simplified topology considering two hops has done by MVPool-SL allows to recover the structure of the reduced graph.

Dataset gPool [START_REF] Gao | Graph u-nets[END_REF] SagPool [START_REF] Lee | Self-attention graph pooling[END_REF] MVPool-SL [START_REF] Zhang | Hierarchical multi-view graph pooling with structure learning[END_REF]

Conclusion

Our graph pooling method MIVSPool is based on a selection of surviving vertices thanks to a Maximal Independent Vertex Set (MIVS) and an assignment of non-surviving vertices to surviving ones. Unlike state-of-art methods, our method allows to preserve the totality of graph information during its reduction.

Fig. 1 :

 1 Fig.1: General architecture of our GNN. Each block is composed of a convolution layer followed by a pooling layer. Features learned after each block are combined to have several levels of description of the graph.

Fig. 3 :

 3 Fig. 3: Evolution of MIVS algorithm on a graph from the NCI1 dataset (Section 4.1). Number inside each vertex corresponds to its score s. Candidate, Survivor and Non-Survivor vertices are respectively denoted by: v f , v f , v f .

Table 1 :

 1 Statistics of datasets

	Dataset	#Graphs	#Classes	Avg |V|	Avg |E|
	D&D	1178	2	284 ± 272	715 ± 694
	PROTEINS	1113	2	39 ± 46	72 ± 84
	NCI1	4110	2	29 ± 13	32 ± 14
	ENZYMES	600	6	33 ± 15	62 ± 26

  Score function MIVSPoolrand MIVSPool T op-k MIVSPool SagP ool MIVSPool M V P ool

	D&D	77.04 ± 0.63	77.10 ± 1.00	75.10 ± 0.76	77.38 ± 0.94
	PROTEINS	75.15 ± 0.44	75.36 ± 0.60	75.11 ± 0.74	75.62 ± 0.47
	NCI1	72.16 ± 0.55	73.82 ± 0.94	73.72 ± 0.71	72.97 ± 0.71
	ENZYMES	45.80 ± 1.35	37.55 ± 1.94	38.68 ± 2.81	46.80 ± 1.53

Table 2 :

 2 Study of MIVSPool according to the score function

	Dataset	Pooling 1	Pooling 2	Pooling 3	Pooling 4	Pooling 5
	D&D	4.1 ± 0.6	3.6 ± 0.6	3.2 ± 0.6	-	-
	PROTEINS 3.1 ± 0.7	2.7 ± 0.6	2.3 ± 0.5	-	-
	NCI1	3.2 ± 0.5	2.9 ± 0.5	2.5 ± 0.5	2.3 ± 0.5	2.1 ± 0.3
	ENZYMES 3.2 ± 0.6	2.7 ± 0.6	2.4 ± 0.5	-	-

Table 3 :

 3 Averaged number of iteration of MIVS for each pooling step.

Table 4 :

 4 MIVSPool MIVSPoolcomp. D&D 75.76 ± 0.82 75.92 ± 0.92 77.26 ± 0.37 77.38 ± 0.94 77.88 ± 0.73 PROTEINS 73.49 ± 1.44 74.30 ± 0.62 75.04 ± 0.67 75.62 ± 0.47 75.81 ± 0.75 NCI1 71.66 ± 1.04 72.86 ± 0.68 74.79 ± 0.58 72.97 ± 0.71 73.44 ± 0.68 ENZYMES 36.93 ± 2.45 35.30 ± 1.80 39.45 ± 2.57 46.80 ± 1.53 45.60 ± 2.37 Comparison of MIVSPool according to other hierarchical methods. Highest and second highest accuracies are respectively in bold and blue.

Acknowledgements:

The work was performed using computing resources of CRIANN (Normandy, France).


Willem Boiten

ECRIN Christian Ohmann

Ayodeji Adeniran

ECRIN Steve Canham

Monica Cano

Abadia " Gauthier

Luisa Chiusano

ERINHA Romain David

EATRIS/IRFMN Maddalena Fratelli

Phil Gribbon

EATRIS/Lygature Jan-Willem Boiten

Gauthier Chassang

EMBRC Maria Luisa Chiusano

EATRIS Rebecca Ludwig

BBMRI), Mihaela Matei (ECRIN Michaela Th Mayrhofer

Mihaela Matei

( Eu-Openscreen

Petr Holub

ELIXIR Arshiya Merchant

ECRIN), Luca Pireddu (CRS4/BBMRI Maria Panagiotopoulou

Alex Sanchez

Irene Schlünder

( Bbmri

EOSC-LIFE WP4 TOOLBOX: Toolbox for sharing of sensitive data -a concept description

Keywords: .2. Workflow 2: requesting and sharing data 6. Navigation and user interface 7. Sustainability & toolbox governance 8 Sensitive data, tags, categorisation, life sciences, toolbox, data sharing, EOSC-Life

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

This project has received funding from the European Union's Horizon 2020 research and innovation programme under grant agreement No 824087.

OBJECTIVE & OVERVIEW OF THE DATA SHARING TOOLBOX

A large segment of the data to be handled and processed by biomedical research infrastructures is sensitive data. Sharing those data is causing an additional layer of complexity on top of the more generic technical and data security issues. The sensitivity of the data may arise from its personal nature (in particular health data), but can also be caused by intellectual property considerations, biohazard concerns, or the Nagoya protocol. In all these cases particular care needs to be taken when sharing or reusing the data, in particular in a cloud context.

The EOSC-Life Toolbox aims to provide guidance to:

1) Researchers or other data providers (e.g., sponsors, institutions or private organisations) wishing to make their sensitive data available for future reuse, i.e., enabling future sharing of data; 2) Researchers and data providers in the case where a researcher wishes to make use of sensitive data made available by the data provider, i.e., enabling actual sharing of data.

The content of toolbox will be derived from two main sources: (1) existing recommendations, procedures, best practices, and links to software (tools) to support data sharing and reuse operations relevant for sensitive data management in the cloud environment, esp. EOSC; (2) guidelines and other useful resources drafted in the context of EOSC-life (e.g., WP7 guidelines on the required capabilities of cloud providers). The toolbox will not be designed de novo; instead, it will help scientists to navigate to previously collected high-quality content available throughout our collective infrastructure landscape. Also, the sustainability of the toolbox will be a design consideration from the start: EOSC-Life is a project with a clear end date. Any toolbox can only be useful and successful if it is maintained and updated beyond that end date.

The creation of the toolbox and filling it with content will be driven by use cases, which may stem from within EOSC-Life (e.g., WP1, WP3, WP6 (Nagoya protocol)) or from the Research Infrastructures' network / other EU-projects. Data sharing issues highlighted by those use cases will prioritize the content collection for the toolbox and will provide test cases for the usability of the toolbox.

The toolbox for data sharing will be the main output for the EOSC-LIFE deliverable 4.3 ("Guidance and policy on standards and tools to facilitate sharing and reuse of multimodal data (including imaging), cohort integration and biosamples").

SCOPE OF THE TOOLBOX

Control of the scope of the toolbox is essential in order to be able to create a useful product within the time and budget available. A high-level overview of this scope in a few bullets:

• Only secondary use of data is considered; in case of primary use only those aspects related to future reuse of the data will be considered.

• The main actor for the toolbox is the researcher with two viewpoints: ○ a researcher who is sharing data or wishes to share the data in the future, or ○ a researcher who is requesting data for secondary use. • The toolbox will not be developed de novo, but assembled from existing resources, preferably by providing pointers to resources rather than copying content. However, the toolbox will provide sufficient context to those pointers to allow users to understand their scope and applicability.

It would be a pitfall to aim for the perfect product from the start. The initial focus should be the "minimum lovable product".

GENERAL OUTLINE OF THE TOOLBOX

The guidance for handling sensitive data in EOSC-Life will be provided via a toolbox where we have three main aspects to consider:

1. Use cases/user stories 2. Content: guidances, tools, and best practice description + some generic background information for novice users 3. Navigation, user interface elements and functionalities to easily access the content The toolbox is a website, where different navigation options should guide users to the right content, not only for simple use cases basically consisting of one data sharing issue, but also for more complex use cases having to deal with multiple issues. Maintenance and sustainability of the website will be an essential design consideration: it should be easy to keep the content upto-date and the responsibility beyond the EOSC-Life project lifetime should be clear from the start. The value of the toolbox highly depends on the reliability and actuality of the content, not only for the toolbox itself, but also for the content from other sources that are referenced. A mechanism to keep control of the actuality of that content should be considered. A disclaimer should be included on the toolbox declining all warranties in case recommendations from the toolbox have caused any damages.

USE CASES & USER STORIES

Prototypic use cases are identified that represent the data sharing needs of our research community. The use cases should be well selected in order to be sufficiently specific to be easily translatable to concrete user needs in practice, while at the other hand the use cases should be sufficiently generic to allow for broad usage. So, it is required to have a representative set of use cases covering a broad range of data sharing problems related to EOSC-Life. The use cases can be derived from the existing and upcoming demonstrators in WP1, 3, and 6, and from specific topics raised by other EOSC-Life work packages.

Concrete data sharing dilemmas will be distilled from these use cases, maintaining the link to the use cases to make the needs tangible. Different approaches will be used to represent the use cases (e.g., textual user stories or forms to be filled in) and to link the use cases to the content.

Currently we have a number of sources of use cases which can be included:

1. The evaluation of the applications for the WP1 open call for use cases. 2. Common use cases with WP2 allowing to iron out the requirements for handling sensitive data in the toolset developed in that work package, in particular for workflows to be used for sensitive data. At this time, only one use case with sensitive data has been identified by WP2: Sharing digital pathology images (led by Petr Holub). 3. The -Beyond One Million Genomes-project (https://b1mg-project.eu/) could offer a valuable source of use cases. Four use case areas have been identified in this project: rare disease, cancer, common complex disease, and infectious disease/COVID-19. Each of these areas need to be more specific to be suited for guiding the toolbox development. 4. A set of use cases previously collected: (Paris; January 23 rd , 2020) as described in a separate working document [START_REF] Boiten | Guidance for sensitive multimodal data in EOSC-Life -minutes of first scoping discussion[END_REF] .

These use cases should be as concrete as possible: real research projects with real data sharing issues. When generalizing too soon we may miss the real-life problems which are often caused by apparently insignificant details.

From these use cases we can derive more generic user stories. This can be done in multiple ways: in the Paris meeting referred to above, we tried to identify the user stories as generic data sharing issues. This has been further prioritized in two sessions as summarized in the minutes 2 from the EOSC-Life WP4 meeting in Brussels (March 4 th , 2020). Alternatively user stories can be summarised in simple statements from the researcher's perspective in an extensive list of examples collected by ECRIN (Steve Canham). These user stories are all formatted in a similar fashion, e.g., "as a researcher, as part of study planning, I want to ensure that consents I obtain will support the agreed data sharing policy".

CONTENT OF THE TOOLBOX

Regarding the content of the toolbox there are two main workflows to be considered which may also be reflected in the user interface:

1. Planning data sharing in a research project, i.e., planning for future sharing 2. Requesting and sharing data for secondary use, i.e., actual data sharing

WORKFLOW 1: PLANNING FOR SHARING

This workflow is summarized in the following workflow diagram: Possible aspects to consider when planning for future sharing:

-learn about data sharing in general (procedures, options, repositories) -check funder/institutional requirements for data sharing -define data sharing procedures in data management plan -provide summary of data sharing procedure in study/research protocol -informed consent -incorporate data sharing summary in study registry (if relevant) -align data sharing plan with all collaborators of the project -learn about ethical and legal frameworks

WORKFLOW 2: REQUESTING AND SHARING DATA

Possible aspects to consider when actually wishing to share data:

-if data sharing is planned prospectively, data sharing plan needs to be reviewed and updated or otherwise a new data sharing plan needs to be developed -decision on de-identification (e.g., anonymisation, pseudonymisation -dependent on explicit consent) -application of de-identification This project has received funding from the European Union's Horizon 2020 research and innovation programme under grant agreement No 824087.

-carry out risk assessment: anticipate risks and plan mitigation measures (e.g., Data

Protection Impact Assessment (DPIA) if requested) -select file formats for data and meta data -explore repository, where the data should be stored (conditions, access tape, costs) -sign data transfer agreement -transfer data to sustainable repository -monitor use of transferred data -submit request (form) to repository -sign data use agreement with repository -download requested data on your computer or use protected analysis environment -describe data sharing tools / conditions / contact / licences / provenances in a data paper -perform secondary analysis -publish secondary analysis -citing the data paper in the secondary publication It is not the intention of the EOSC-Life project to develop new guidelines, best practices, tools/services. Instead, existing content should be assessed according to standardised criteria relevant for the EOSC-Life project. However, in case of certain white spots, EOSC-Life WP4 may organize expert workshops to gain access to specific content for the toolbox. An example of such a workshop was the meeting on pseudonymization/anonymization organized in Paris (January 22-23, 2020). The draft report 3 of this meeting is available.

The BBMRI-ERIC Knowledge Base (as described in a publication 4 ) may provide an important starting point for the design, relevant content, and sustainability aspects of the toolbox: https://www.bbmri-eric.eu/elsi/knowledge-base/. AA demonstration of that Knowledge Base for this working group was organized July 27 th , 2020. A new release of the Knowledge Base was issued in autumn 2020. The next version update is scheduled for Q2 2021.

EOSC-Life WP7 is currently working on a guideline for the required capabilities of cloud providers when wishing to process sensitive data. The resulting guidelines should certainly be included in the toolbox.

NAVIGATION AND USER INTERFACE

The navigation is essential to the use of the toolbox and should be based upon the needs of a user related to a specific problem or use case. EOSC-Life WP9 is supporting the usability aspects of the toolbox implementation ("user journey"). A usability webinar (September 28 th , 2020) has been organized to inform the toolbox team regarding tips and tricks in this realm. It was based on the EBI Training on user experience: https://www.ebi.ac.uk/trainingbeta/online/courses/user-experience-design/what-is-user-experience-design/ One possible methodological approach could be to follow the approach outlined in the book by J.J. Garrett [START_REF] Garrett | The Elements of User Experience: User-Centered Design for the Web and Beyond, Second Edition[END_REF] distinguishing product as a functionality and product as information, covering different planes.

Another example that could inspire us with the further creation of the toolbox is the Data Stewardship Wizard created for the RDA COVID-19 initiative. A demo of this toolbox was organized for the larger EOSC-Life Toolbox working group on August 24 th , 2020.

We will start a small working party that is dedicated to the technology to be used for the toolbox. The selected technology should be based on open-source software to avoid vendor lock-in (whenever possible). Google Analytics is an important tool to be used to obtain an indication of the actual usage of the toolbox. -A hide/show button for each functionality to permit users to personalize /simplify the interface The actual UX design should ideally be tested in practice with representatives from the intended user group. Also, an on-line feedback option could serve as a source of user feedback during the production phase of the toolbox. It may also be a way to capture reports on outdated content. Direct search may be supported by tagging the content. This basic approach to navigation could be supported by a form/questionnaire to be filled in, which characterises the main dimensions of data sharing. The content (e.g., guidances, best practice, tools) could be tagged with these characteristics and thus made identifiable via search. A detailed proposal [START_REF] Ohmann | EOSC-Life WP4 Toolbox: Categorisation system for resources to be referenced in the toolbox for sharing of sensitive data[END_REF] for such a tagging system has been created, as well as a study protocol 7 to test its practical usability. This could be a first higher layer for discovery of content, which may be extended by more specific navigation, if needed. A question to be answered is whether and how the different navigation strategies (wizard, mind-mapping, direct search, form) could work together and whether tagging/annotation of content can be used to serve different navigation strategies. In fact, a hybrid approach combining different techniques will probably be preferable.

SUSTAINABILITY & TOOLBOX GOVERNANCE

The usefulness of the toolbox is not only determined by the usability of the toolbox, but also by the quality, the reliability, and the actuality of its content. It is therefore crucial to have a (lightweight) governance process in place, e.g., by reviewing every toolbox topic by two or three independent experts before publishing it. One of the review criteria should cover the predicted sustainability of the content. Previous experience shows that updates of documents are often not officially released anymore [START_REF] Sariyar | Sharing and Reuse of Sensitive Data and Samples: Supporting Researchers in Identifying Ethical and Legal Requirements[END_REF] . It should be unlikely that the content will be outdated rapidly; nevertheless, some periodic review will be needed to maintain the content of the toolbox over time. Preferably, this should be automated as far as possible, but experience from BBMRI's Knowledge Base tells that this is hard to realise in practice. For that reason, BBMRI is using a manual 3-monthly update cycle for the content on the ELSI service page on the BBMRI-ERIC website.

  This project has received funding from the European Union's Horizon 2020 research and innovation programme under grant agreement No 824087.

Figure 1 :

 1 Figure 1: workflow when planning for future data sharing of sensitive data.

Figure 2

 2 Figure 2 aims to visualize these aspects in a workflow:

Figure 2 :

 2 Figure 2: Key aspects to consider when making sensitive data available for sharing.

A

  few generic considerations for the navigation: -In terms of usability, cross-references should be used with care. When too many clicks are involved users tend to give up. -Management of the content should be independent of the navigation. It should be possible to update the content without breaking the navigation. Multiple navigation approaches need to be supported; possible options are: -Decision tree; it has the clear advantage that steps are executed in the right order -Direct search -Infographic -Classic alphabetical index -Landing page with most recent additions -A toolbox in a nutshell page with most useful / used information This project has received funding from the European Union's Horizon 2020 research and innovation programme under grant agreement No 824087.

  The toolbox can consist of recommendations, procedures, links to tools, and examples/best practices. Preferably the best practices are derived from the use cases.

	• A certain level of expertise is expected from the targeted audience, i.e., "intermediate
	level of expertise" is expected: neither a novice in data sharing nor an ethical/legal
	expert. It should also be considered that users may be experienced in one aspect, while
	being novice in other areas. Although not really in scope for the toolbox, we do intend
	to provide some resources for novice users.
	• The content of the toolbox should be of relevance for cloud usage of sensitive data in
	general, and the European Open Science Cloud in particular. The priority for filling the
	toolbox will be set by either:
	○ representative use cases, or
	○ specific topics raised by other EOSC-Life work packages, or
	○ specific input from the WP4 EOSC-Life toolbox team.
	•

This project has received funding from the European Union's Horizon 2020 research and innovation programme under grant agreement No 824087.

This project has received funding from the European Union's Horizon 2020 research and innovation programme under grant agreement No 824087.