Serge Iovleff
email: serge.iovleff@univ-lille.f

Seydou Nourou Sylla
email: seydou.sylla@pasteur.sn

Cheikh Loucoubar
email: cheikh.loucoubar@pasteur.sn

Block clustering of Binary Data with Gaussian Co-variables

The simultaneous grouping of rows and columns is an important technique that is increasingly used in large-scale data analysis. Motivated by the analysis of an epidemiological and genetic data set, we present a novel coclustering method using co-variables in its construction. It is based on a latent block model taking into account the problem of grouping variables and clustering individuals by integrating information given by a set of co-variables. Numerical experiments on simulated data sets and results on the genetic data set highlight the interest of this approach.

Introduction

Clustering is a data analysis method that aims to group together a set of observations into homogeneous classes. It plays an increasingly important role in many scientific and technical fields. Its aim is the automatic problems solving by decision-making based on the observations and to define rules for classifying objects depending on qualitative or quantitative variables. Clustering is the most popular technique for data analysis in many disciplines. In recent years, co-clustering has been increasingly used.

Unlike classical clustering, which groups similar objects from a single collection of objects, co-clustering or bi-clustering [START_REF] Madeira | Biclustering algorithms for biological data analysis: a survey[END_REF] aims to group two types of entities simultaneously, based on similarity of their pairwise interactions. It is most often used with bipartite spectral graphing partitioning methods in the field of extracting text data [START_REF] Dhillon | Co-clustering documents and words using bipartite spectral graph partitioning[END_REF] by simultaneously grouping documents and content (words). It is used for analyzing huge corpora of unlabeled documents [START_REF] Xu | Co-clustering analysis of weblogs using bipartite spectral projection approach[END_REF] in order to simultaneously understand aggregates of subsets of web users sessions and information from the page views. Co-clustering algorithms have also been developed for computer vision applications: it is used for grouping images simultaneously with their low-level visual characteristics and for content-based search [START_REF] Guan | Spectral images and features co-clustering with application to content-based image retrieval[END_REF].

In this paper we extend co-clustering methods allowing simultaneous detection of associations between variables and individuals by taking into account covariables. Our method is to be used when one want to co-cluster a set X of (binary) variables, and individuals in coherence with independent (continuous) variables Y measured on these same individuals.

This co-clustering approach is motivated by a malaria data set from Senegalese populations [START_REF] Trape | The rise and fall of malaria in a west african rural community, dielmo, senegal, from 1990 to 2012: a 22 year longitudinal study[END_REF]). We want to co-cluster a set X of binary variables, i.e. presence of genetic variants in SNPs (Single-Nucleotide Polymorphism), and individuals (patients) in coherence with a quantitative variable Y measured on these same individuals. In the biological application we have in mind, the additional measure of interest Y should be taken into account by the co-clustering process in order to obtain significant results.

In biology, the selection in a set of variable X those associated with a quantitative outcome Y , is generally done using genome-wide association studies (GWAS). In a typical GWAS there is measure of hundreds of thousands, or millions, of genetic variants (SNPs), and the attempt is to identify regions harboring SNPs that affect some phenotype of interest, see for example [START_REF] Timmann | Genome-wide association study indicates two novel resistance loci for severe malaria[END_REF]) for a GWAS applied to severe falciparum malaria in patients and controls from Ghana, West Africa. This goal can naturally be cast as a variable selection regression problem, with the SNPs as the covariates in the regression. Using the co-clustering presented here, the problem is reversed: we find cluster of individuals and SNPs using the phenotype Y as prior and then find the most influential genetic variants by inspecting the posterior distribution (see parts 2.6 and 3.2).

A second goal in our study is to dichotomize the initially quantitative outcome Y into a binary variable reflecting two categories of risk or two different levels of severity of the disease. Usually, the median value on Y or the value 0 are used as threshold [START_REF] Loucoubar | Detecting multi-way epistasis in family-based association studies[END_REF]).

Using, the co-clustering approach presented here, the outcome variable (Y ) is automatically clustered in two groups and this allow us to find the optimal cut off to its binarization, i.e. the partition of Y that maximize its association with determinant features.

The paper is organized in two parts. In the first part (section 2) we develop the models, formulas and algorithms. We explain the principle of block mixture models through section 2.1. The latent block model for binary variable taking into account continuous co-variables and the model parameters estimation are proposed in section 2.2. The parameter estimation method is described in section 2.3 and 2.4. The choice of the optimal number of blocks and the measure of influence of each variable on the co-variable Y are presented in section 2.5 and 2.6. In the second part (section 3) the method is illustrated on simulated data (section 3.1) and on real genetic data (section 3.2).

2 Block mixture models

Classical latent block model

Let x be a data set doubly indexed by a set I with n elements (individuals) and a set J with m elements (variables). We represent a partition of I into g clusters by z = (z 11 , . . . , z ng ) with z ik = 1 if i belongs to cluster k and z ik = 0 otherwise, z i = k if z ik = 1 and we denote by z .k = ∑ i z ik the cardinality of row cluster k. Similarly, we represent a partition of J into d clusters by w = (w 11 , . . . , w md ) with w j = 1 if j belongs to cluster and w j = 0 otherwise, w j = if w j = 1 and we denote by w . = ∑ j w j the cardinality of column cluster .

The block mixture model formulation is defined in [START_REF] Govaert | Clustering with block mixture models[END_REF] and [START_REF] Bhatia | blockcluster: An r package for model-based co-clustering[END_REF] (among others) by the following probability density function

f (x; θ ) = ∑ u∈U p(u; θ ) f (x|u; θ )
where U denotes the set of all possible labels of I × J and θ contains all the unknown parameters of this model. By restricting this model to a set of labels of I × J defined by a product of labels of I and J, and further assuming that the labels of I and J are independent of each other, one obtains the decomposition

f (x; θ ) = ∑ (z,w)∈Z ×W p(z; θ )p(w; θ ) f (x|z, w; θ ) (1)
where Z and W denote the sets of all possible labellings z of I and w of J. Equation (1) defines a Latent Block Model (LBM).

2.2 LBM for binary variables with co-variables: General formulation

From now, we assume that x is a binary data set. Let y represent a data-set (covariables) of R p indexed by I. In order to take into account this set of co-variables, the classical block model formulation is extended to propose a block mixture model defined by the following probability density function

f (x, y; θ ) = ∑ (z,w)∈Z ×W p(z; θ )p(w; θ ) f (x, y|z, w; θ ). (2)
By extending the latent class principle of local independence to our block model, data pairs (x i j , y i ), i = 1, . . . , m, j = 1, . . . , n are independent once z and w are fixed. Hence we have

f (x, y|z, w; θ ) = ∏ i, j f (x i j , y i |z i , w j ; θ ).
We choose to model the dependence between x i j and y i using the canonical link for binary response data f (x i j |y i , β z i w j ) = logis(β 0,z i w j + β T z i w j y i ) x i j 1logis(β 0,z i w j + β T z i w j y i )

1-x i j = e

x i j (β 0,z i w j +β T z i w j y i )

1 + e β 0,z i w j +β T z i w j y i

(3) with (β 0 , β k,l ) ∈ R p+1 and logis(x) = e x /(1 + e x ). Data points y i , i = 1, . . . , m are independent once z is fixed. In the examples presented in section 3, we choose

f (y|z; θ ) = ∏ i φ (y i ; µ z i , Σ z i )
with φ denoting the multivariate Gaussian density in R p . In order to simplify the notation, we add a constant coordinate 1 to vectors y i and write β k,l in the latter rather than (β 0,k,l , β k,l ).

The parameters are thus θ = (π, ρ, β , µ, Σ), where π = (π 1 , . . . , π g ) and ρ = (ρ 1 , . . . , ρ d ) are the vectors of probabilities π k and ρ that a row and a column belong to the kth row component and to the th column component respectively, β = (β kl ) are the coefficients of the logistic function, µ and Σ are the means and variances of the Gaussian density. In summary, we obtain the latent block mixture model with pdf

f (x, y|θ ) = ∑ (z,w)∈Z ×W ∏ i, j π z i ρ w j e x i j (y T i β z i w j ) 1 + e y T i β z i w j φ (y i ; µ z i , Σ z i ).
Using the above expression, the randomized data generation process can be described by the four step row labelling (R), column labelling (C), co-variable data generation (Y) and data generation (X) as follows:

(R) Generate labels z = (z 1 , . . . , z n ) according to the distribution π = (π 1 , . . . , π g ). (C) Generate labels w = (w 1 , . . . , w m ) according to the distribution ρ = (ρ 1 , . . . , ρ d ). (Y) Generate for i = 1, ..., n vector y i according to the Gaussian distribution

N p (µ z i , Σ z i ).
(X) Generate for i = 1, ..., n and j = 1, ..., m a value x i j according to the Bernoulli distribution f (x i j |y i ; β z i w j ) given in (3).

Model Parameter Estimation

The complete data is represented as a vector (x, y, z, w) where unobservable vectors z and w are the labels. The log-likelihood to maximize is

l(θ ) = log f (x, y; θ ) (4)
and the double missing data structure, namely z and w, makes statistical inference more difficult than usual. More precisely, if we try to use an EM algorithm as in standard mixture model [START_REF] Dempster | Maximum likelihood from incomplete data with the em algorithm (with discussion)[END_REF] the complete data log-likelihood is found to be

L C (z, w, θ ) = ∑ k z .k log π k + ∑ w . log ρ + ∑ i, j,k, z ik w j log f (x i j , y i ; θ k ). (5)
The EM algorithm maximizes the log-likelihood l(θ ) iteratively by maximizing the conditional expectation Q(θ , θ (c) ) of the complete data log-likelihood given a previous current estimate θ (c) and (x, y):

Q(θ , θ (c) ) = E L C (z, w, θ ) x, y, θ (c) = ∑ i,k t (c) ik log π k + ∑ j, r (c)
j log ρ + ∑ i, j,k, e (c) ik j log f (x i j , y i ; θ k ) where t (c) ik = P(z ik = 1|x, y, θ (c) ), r (c) jl = P(w j = 1|x, y, θ (c) ), e (c)
ik j = P(z ik w j = 1|x, y, θ (c) ). Unfortunately, difficulties arise due to the dependence structure in the model, in particular to determine e (c) ik j . The assumed independence of z and w in (1) is not conserved by the posterior probability.

To solve this problem an approximate solution is proposed in [START_REF] Govaert | Clustering with block mixture models[END_REF] using the [START_REF] Hathaway | Another interpretation of the em algorithm for mixture distributions[END_REF] and [START_REF] Neal | A View of the Em Algorithm that Justifies Incremental, Sparse, and other Variants[END_REF] interpretation of the VEM algorithm. Consider a family of probability distribution q(z ik , w j ) verifying q(z ik , w j ) > 0 and the relation q(z ik , w j ) = q(z ik )q(w j ), for all i, j, k, l. Set t ik = q(z ik ) and r jl = q(w j ), t = (t ik ) ik for i = 1, . . . , n, k = 1, . . . , g and r = (r jl ) jl for j = 1, . . . , m and l = 1, . . . , d. One easily shows that

l(θ ) = FC (t, r; θ ) + KL(q(z, w) p(z, w|x, y, θ )) (6)
with KL(q p) denoting the Kullback-Liebler divergence of distribution p and q,

FC (t, r; θ ) = ∑ k t .k log π k + ∑ r . log ρ l + ∑ i, j,k, t ik r j log f (x i j , y i ; θ k ) + H(t) + H(r) (7)
and H(t), H(r) denoting the entropy of t and r, i.e.

H(t) = ∑ ik t ik logt ik , H(r) = ∑ jl r jl log r jl .
FC is called the free energy or the fuzzy criterion. As the Kullback-Liebler divergence is always positive, the fuzzy criterion is a lower bound of the log-likelihood and is used as a replacement for it. Doing that, the maximization of the likelihood l(θ ) is replaced by the following problem argmax t,r,θ

FC (t, r, θ ).
This maximization can be achieved using the BEM algorithm detailed hereafter.

Block expectation maximization (BEM) Algorithm

The fuzzy clustering criterion given in (7) can be maximized using a variational EM algorithm (VEM). We here outline the various expressions evaluated during E and M steps.

π k = t .k n and ρ = r . m . (8)
Also, for t, r fixed, the estimate of model parameters β is obtained by maximizing

β kl = argmax β ∑ i j t ik r jl log f (x i j |y i ; β ), k = 1, . . . , g, l = 1, . . . , d. (9)
Details are given in appendix B. Finally parameters of the Gaussian density are given by the usual formulas

µ k = 1 t .k ∑ i t ik y i and Σ k = 1 t .k ∑ i t ik (y i -µ k )(y i -µ k ) T . (10)
Putting everything together, we obtain the BEM algorithm.

BEM algorithm: Using the E and M steps defined above, BEM algorithm can be decomposed as follows:

Initialization Set t (0) , r (0) and θ (0) = (π (0) , ρ (0) , β (0) , µ (0) , Σ (0) ).

(a) Row-EStep Compute t (c+1) using formula c+1) , Σ (c+1) using equations ( 8) and ( 10) and estimate β (c+1/2) by solving maximization problem ( 9).

t (c+1) ik = π (c) k ∏ jl f (x i j |y i ; β (c) kl )φ (y i ; µ (c) k , Σ (c) k ) r (c) jl ∑ k π (c) k ∏ jl f (x i j |y i ; β (c) kl )φ (y i ; µ (c) k , Σ (c) k ) r (c) jl . ( 11
) (b) Row-MStep Compute π (c+1) , µ (
(c) Col-EStep Compute r (c+1) using formula r (c+1) jl = ρ (c) l ∏ ik f (x i j |y i ; β (c+1/2) kl ) t (c+1) ik ∑ l ρ (c) l ∏ ik f (x i j |y i ; β (c+1/2) kl ) t (c+1) ik . ( 12
)
Observe that r jl does not depend of the density of y.

Selecting the number of blocks

BIC is an information criterion defined as an asymptotic approximation of the logarithm of the integrated likelihood [START_REF] Schwarz | Estimating the dimension of a model[END_REF]). The standard case leads to write BIC as a penalised maximum likelihood:

BIC = -2 max θ l(θ ) + D log(N)
where N is the number of statistical units and D the number of free parameters and l(θ ) defined in (4). Unfortunately, this approximation cannot be used for LBM, due to the dependency structure of the observations (x, y). However, a heuristic have been stated to define BIC in [START_REF] Keribin | Model selection for the binary latent block model[END_REF] and [START_REF] Keribin | Estimation and selection for the latent block model on categorical data[END_REF]. BIC-like approximations ICL lead to the following approximation as n and m tend to infinity

BIC(g, d) = -2 max θ log f (x, y; θ ) + (g -1) log n + (d -1) log m + λ log n + gd(p + 1) log(mn)
with λ the number of parameters of the y distribution. For LBM, the intractable likelihood f (x, y; θ ) is replaced by the maximized free energy FC in (7) obtained by the BEM algorithm.

Measuring the Influence of a Variable

Let j be fixed (a column of the matrix x). We would like to measure the effect of the variable x j = (x i j ) n i=1 on y. It is possible to obtain a measure of this effect by looking to the posterior probability of y.

Lemma 1 Let (x, z, w) be fixed. For l = 1, . . . , d let m l denotes the number of columns with label l, i.e m l = #{w jl = 1, j = 1, . . . m} and for a row i fixed let m il denotes the number of elements such that w jl = 1 and x i j = 1, i.e. m il = #{w jl x i j = 1, j = 1, . . . n}. The posterior probability of the co-variable y is

f (y|x, z, w, θ ) ∝ n ∏ i=1 d ∏ l=1 π z i ρ m l l logis(y T i β z i l ) m il 1 -logis(y T i β z i l ) m l -m il φ (y i ; µ z i , Σ z i ) ∝ n ∏ i=1 π z i φ (y i ; µ z i , Σ z i ) d ∏ l=1 ρ m l l e m il y T i β z i l 1 + e y T i β z i l m l
Alternatively, for k = 1, . . . , g, let n k denotes the number of rows with label k, i.e. n k = #{z ik = 1, i = 1, . . . , n}. The posterior probability of the co-variable y is

f (y|x, z, w, θ ) ∝ m ∏ j=1 ρ w j g ∏ k=1 π n k k ∏ i:z i =k e x i j (y T i β z i w j ) 1 + e y T i β z i w j φ (y i ; µ k , Σ k ).
The proof of this lemma is straightforward and therefore omitted. Assuming z and w known, we measure the influence of a variable using its contribution to the posterior probability. Fixing j, taking the logarithm and eliminating terms independent of x j , we obtain the influence measure criteria

I( j) = log ρ w j + n ∑ i=1 x i j y T i β z i w j -log(1 + exp(y T i .β z i w j )) . (13)
Equation ( 13) is interpreted as the log-contribution to the posterior distribution of the variable x j . Replacing the unknown labels w j and z i by their MAP estimators ŵ j and ẑi , we are able to sort the variables from the most to the less influential.

Numerical Experiments

Simulated data

In this part, we fix g = 2 (the number of cluster of the individuals) and we study the number of correctly classified rows/columns for different values of m, n, d. This study is only partial and is focused to values "near" the dimensions of the real data set which motivate this work (section 3.2).

We also conduct experiments about the computational time when the number of columns is huge. These experiments were performed for a future work and are only given for information. The real data set illustrating this paper has a much less number of columns.

Error rate

We simulate 80 times data set and average the number of columns correctly classified The cluster of a column is estimated using the maximum a posterior (MAP) estimator ŵ j = arg

d max l=1 r jl .
Comparison between the graphics in figure 1 shows that the number of incorrectly classified columns labels increases as d increases while it remains relatively constant with m. An other salient feature is that when the number of individuals m is greater, this error rate is lower. The number of correctly classified rows is stable near 0.9 for all tested configurations of the parameters and is not displayed.

Computational time

We compute 80 times the elapsed time of the estimation procedure for various configurations of the parameters on a HP Zbook G3 with Intel Core i7-6700HQ (2.60 GHz, 2133 MHz, 6 MB L3 cache, 4 cores, 45W). The (averaged) computing time as a function of m when g = 2 for different values of m (the number of columns) and when d (the number of cluster in columns) take values 2 and 6 is plotted in figure 2 below We can observe that as n grows the elapsed time grows linearly, but that the slope increases as d (the number of class in columns) is increased.

Real Data Analysis

Here, we study data from an epidemiological and genetic survey of malaria disease in Senegal. Data were collected between 1990 to 2008. We worked on a dataset including n = 885 individuals with measured malaria risk score (phenotype) and with genotypes available on several candidate genes for susceptibility/resistance to the disease. A total of m = 45 Single Nucleotide Polymorphisms (SNPs) were considered across these genes and were used as genetic variables. The malaria risk score was a quantitative measure normally distributed and was considered as a covariable for this co-clustering method. The SNPs are coded in dominant effect on the disease risk. Using the BIC criteria (figure 3), we choose to focus on the model with d = 2 groups of individuals and g = 11 groups of SNPs.

Analysis for phenotype data

We used the malaria risk score as phenotype. The choice of a mixture model or not depends on the application context. In the case of genetic epidemiology, we are often interested in the comparison between the susceptible and the resistant to a given disease. Here, we looked for genes that could explain differences between susceptible and resistant, justifying the use of a mixture model on this phenotype.

After block-clustering, we found that individuals were divided into two groups: Individuals with lowest malaria risk scores (-2.101 to 0.040) defined as resistant and individuals with highest malaria risk scores values (0.042 to 3.509) defined as susceptible.

In figure 4 we observe on the left a bi-modal distribution of the phenotype after clustering (i.e. when y is conditioned by (x, z)). On the right we observe how SNP with high/low levels of mutations are grouped together.

Analysis for genotypes data

We looked at the SNPs to determine which ones would potentially be involved in malaria susceptibility / resistance.

The proposed methodology allowed the selection of the most significant SNPs according to the influence measure proposed in section 2.6. The most influential SNPs belong to class 1 and class 9. It is noted that the SNPs of these classes have been shown in the literature to have a high significance effect on malaria. Most Glucose-6-Phosphate Dehydrogenase (G6PD) and hemoglobin SNPs are grouped into these 2 classes. Reviews from exiting literature give us: G6PD deficiency is 2015), Toure, Konate, Sissoko, Niangaly, Barry, Sall, Diarra, Poudiougou, Sepulveda, Campino, Rockett, Clark, Thera, Doumbo, and in collaboration with The MalariaGEN Consortium (2012). Hemoglobins S and C (HbS and HbC respectively) are well known to be two variant forms of normal adult hemoglobin (HbA) resulting from distinct mutations in the β -globin gene. The protective effect of HbS against Plasmodium falciparum malaria has been shown by several authors [START_REF] Beet | Sickle cell disease in the balovale district of northern rhodesia[END_REF], Allison (1954b,a), [START_REF] Williams | How do hemoglobins s and c result in malaria protection?[END_REF]. In the case of HbC, the protection is highest in homozygous individuals with HbC.

The results found with our co-clustering method confirm the link between malaria and sickle cell polymorphism (Hbs) and G6PD.

Association between phenotype and genotypes

The most common approach used to screen association between genetic data and phenotypic data is the GWAS method (Genome Wide Association Studies). GWAS usually performs linear regressions of a quantitative phenotype on each of the genotype variables.

Here, we proposed a co-clustering method, that can state as a step prior to association analysis. The method reorganizes data by simultaneously identifying optimal partitions of the phenotype and clusters of markers so that this partition of the phenotype would be better explained by a given set of markers.

In the co-clustering method, the phenotype is used as co-variable to find optimal clusters on the variables as well as on the individuals. In our Senegalese malaria data, we obtained a dichotomy of the phenotype. This dichotomy allowed us to divide individuals into two categories: susceptible and resistant. In this part, we compare results of GWAS studies between classical approaches and the one using co-clustering as a first step.

• Classical approaches:

testing association between the raw phenotype and genotypes by simple linear regression testing association between the binarized phenotype (by choosing a cutoff, e.g. value zero) and genotypes by logistic regression. • Using co-clustering followed by association tests: testing association between the binarized phenotype (using the optimal cutoff identified by the co-clustering algorithm, here value 0.042) and genotypes by logistic regression.

The association signal between the phenotype and each of the genotype variables (here the SNPs) is represented by the corresponding p-value. Lower the p-value, better the association signal. Number of good signals (i.e. SNPs with p-values less than 0.05) is compared between each approach (liner regression, logistic regression using 0 as cut off to binarize Y, logistic regression using the optimal cut off identified through co-clustering of the data). It is presented in Figure 6 at different level of significance of the signal ( p-value 5e-4 ; 5e-4 ¡ p-value 5e-3 ; 5e-3 ¡ p-value 5e-2). Figure 6 shows that there is more signals at the 5% threshold for the phenotype from co-clustering compared to the two other phenotypes.

In summary the proposed methodology allows to detect more significant SNPs compared to classical methods by increasing detection power through optimal clustering of the data.

Conclusion

In this article, our main contribution is to develop a co-clustering model taking into account a (mixture of) Gaussian co-variables. Applications have been made on simulated and real data sets. Our preliminary results are confirmed in previous studies in Africa. The method offers good classification performance on complex data sets (large number of variables and classes). This method can be useful in a wide variety of classification problems with Gaussian predictors and will allow us to discover new patterns of genes allowing to understand and evaluate the mechanism existing between genetics and malaria in an African population particularly in a Senegalese rural area. Further analysis could be done with more SNPs in another paper in preparation. Estimation is performed using a R package (with computational part in C++) that will be soon be available on the CRAN website https://cran.r-project.org/. Meanwhile the package is available on demand to the authors.

A Computing the (rows and columns) E-Step

For the E-Step t ik value maximize the fuzzy criterion given in equation ( 7). Derivative with respect to t ik gives ∂ FC (t, r; θ ) ∂t ik = log π k + ∑ j, r j log f k (x i j , y i ; θ )logt ik -1.

Equating this equation to zero, taking exponential and recalling that ∑ k t ik = 1, we obtain that t ik is updated as

t (c+1) ik = π (c)
k ∏ j,l f (x i j , y i ; θ (c) ) r (c) jl ∑ k ∏ j,l f (x i j , y i ; θ (c) ) r (c) jl

.

For numerical reason, we prefer to compute the logarithm of this expression which is log(t

(c+1) ik ) ∝ log(π (c) k ) + ∑ j,l r (c)
jl log f (x i j , y i ; θ (c) ).

Recall that (see equation 3) log f (x i j |y i ; β .

Observe that the Gaussian distribution does not depend of j nor l. This term become constant when summing over i and k and disappears when r jl values are normalized.

B Computing the M-Step

For the M-Step, we use a Newton-Raphson algorithm in order to solve the equation ( 9). For each pair (k, l) the function to maximize can be written k,l (β ) = ∑ i, j r jl t ik x i j y T i βr jl t ik log(1 + exp(y T i .β ))

The first derivative with respect to the d-th coordinate β d is

∂ k,l (β ) ∂ β d = ∑
i, j r jl t ik x i j y i,dr jl t ik y i,d exp(y T i β ) 1 + exp(y T i β )

giving the following expression for the gradient giving the following expression for the hessian

∇ β k,l (
H β = -Y t DWY with W = diag r .l exp(y T i .β ) (1 + exp(y T i β )) 2 = diag (r .l µ i (1 -µ i ))

  (d) Col-MStep Compute ρ (c+1) using equations (8) and estimate β (c+1) by solving maximization problem (9). Iterate (a)-(b)-(c)-(d) until convergence.

Figure 1 :

 1 Figure 1: Rates of correctly classified columns when the number of rows is 400 and 800. The number of columns (m) is between 40 and 80. The number of cluster (d) is between 6 and 12. There is only two groups of rows.

Figure 2 :

 2 Figure 2: computational elapsed time for m = 2000, 6000, 10000, 14000 and 18000 (in minutes) and for various values of n.

Figure 3 :

 3 Figure 3: BIC computation for different values of d and g. We observe that it is minimal for g = 2 and d = 11 among tested d values (2, . . . , 14) and g values (2, 3, 4) .

Figure 4 :

 4 Figure 4: Co-Clustering Results: phenotype and SNPs. (a) -Empirical Distribution of the phenotype (histogram) -Distribution of the susceptible (red) -distribution of the resistant (green) -mixing distribution (grey). (b) Array with the presence/absence of mutations before and after block-clustering

Figure 5 :

 5 Figure 5: Representation of each block variable according to the influence measure. Here, block 1 has the highest influence measures (417 to 553) while block 7 has the lowest (1 to 6)

Figure 6 :

 6 Figure 6: Number of significant p-values of each method

,

  β ) = Y T D(Xµ) with Y = [y i ] N i=1 , X = ∑ j r jl x D = diag(t ik ) N i=1The second derivative with respect to β d and β d is

  E-Step: we compute either the values of t (respectively r) with r (respectively t) and θ fixed (formulas (11), (12) hereafter). Details are given in appendix A.M-Step: we calculate row proportions π and column proportions ρ. The maximization of FC w.r.t. π, and w.r.t ρ, is obtained by maximizing ∑ k t .k log π k , and ∑ r . log ρ respectively, which leads to

* This work was supported by LIRIMA,


Ana Alonso Rodríguez

Francesca Rapetti
email: francesca.rapetti@univ-cotedazur.fr

Elena Zappon
email: elena.zappon@polimi.it

Ana Alonso Rodríguez

New degrees of freedom for high-order Whitney approximations of Darcy's flows

Keywords: Darcy equation, Whitney finite elements, high order, new degrees of freedom, graph theory, exterior differential calculus Mathematics Subject Classification (2010) 65N30, 76M10, 05C05

We consider the flow of a viscous incompressible fluid through a rigid homogeneous porous medium. We propose a high-order discretisation based on Whitney finite elements, namely, Raviart-Thomas finite elements of degree r + 1 for the discharge and discontinuous piecewise polynomial finite elements of degree r for the pressure, with r ≥ 0. We comment on the use of new degrees of freedom that have a clear physical meaning, the so-called weights on the small simplices, for the involved discharge and pressure fields. We describe a new numerical strategy to solve the discrete problem based on a tree-cotree block-decomposition of the unknowns, that is natural when considering these new degrees of freedom. Preliminary numerical tests in two dimensions confirm the stability of the adopted method and the effectiveness of the new degrees of freedom.

Introduction

In this paper, we consider Darcy's problem in the framework of the finite element exterior calculus and propose new degrees of freedom for the high order mixed finite element discretisation of this problem on simplicial meshes.

Darcy's problem describes the flow of a viscous incompressible fluid in a porous medium. It is a well-known problem from the mathematical point of view (see, e.g., [START_REF] Bernardi | Discrétisations variationnelles de problèmes aux limites elliptiques[END_REF], [START_REF] Quarteroni | Numerical approximation of partial differential equations[END_REF]), and in scientific journals, it is among the most studied ones, often coupled with either the heat equation or fracture mechanics. Darcy's problem is chosen on purpose here, the idea underlying this work is indeed to show that the documented mathematical pitfalls characterising this problem can be straightforwardly avoided by relying on compatible approaches as those suggested by the finite element exterior calculus (FEEC).

The FEEC uses tools from differential geometry and homological algebra to define, or suggest, the finite element discretisation that is the most compatible with the geometric and topological structure of the PDE problem under exam (see, e.g., [START_REF] Arnold | Finite Element Exterior Calculus[END_REF], [START_REF] Arnold | Finite element exterior calculus, homological techniques, and applications[END_REF], [START_REF] Arnold | Finite element exterior calculus: from Hodge theory to numerical stability[END_REF] and the references cited therein). By reformulating the continuous problem in terms of differential forms, one gets some indications on the discrete formulation with a choice of finite element spaces that is natural, only dictated by the nature of the involved fields and by the mathematical structure of the considered problem. More specifically, the choice of the degrees of freedom that are suggested by the FEEC approach reflects the nature (and global regularity properties) of the fields they represent. When structural features of the continuous problem are preserved at the discrete level, as it is the case with compatible discretisations, stability and consistency properties required for the convergence are guaranteed. This is detailed in the present work for the considered problem and supported by numerical results.

Natural finite element spaces suggested by the FEEC approach are Whitney's differential k-forms P - r+1 Λ k (T ) of polynomial degree r + 1 [START_REF] Arnold | Finite element exterior calculus: from Hodge theory to numerical stability[END_REF], that are the Raviart-Thomas RT r+1 (for k = 2) and the discontinuous DP r (for k = 3) FE spaces, respectively, when passing to vector formalism [START_REF] Quarteroni | Numerical approximation of partial differential equations[END_REF]. Here we rely on these FE spaces with a new set of degrees of freedom on simplices, firstly introduced in [START_REF] Rapetti | Whitney forms of higher degree[END_REF]. For a k-form w ∈ P - r+1 Λ k (T ), classical degrees of freedom are the traditional moments (see, e.g., [START_REF] Arnold | Finite element exterior calculus, homological techniques, and applications[END_REF]). The new degrees of freedom are integrals of w on small k-simplices, a set of subsimplices obtained by suitable contractions of the original mesh elements. These contractions are homothetic transformations on the space, namely the new subsimplices have smaller size while keeping the same shape and orientation of the original mesh elements. This allows for adopting degrees of freedom that reflect the nature of the fields they are associated with and to preserve the meaning of degrees of freedom as cochains, for r > 0, in the spirit of Whitney's original work [START_REF] Whitney | Geometric integration theory. 400[END_REF] dating back to 1957. For r = 0, this is known since the 80s, when Bossavit made the connection between Whitney's forms P - 1 Λ k (T ) used by geometers and the first kind of Nédélec H(curl) and H(div) conforming elements used (see, e.g., [START_REF] Bossavit | Computational electromagnetism: variational formulations, complementarity, edge elements. 375[END_REF]). By choosing the weights as degrees of freedom for finite element discretisations, we can adopt, in fluid mechanics, techniques that are well-known in electro-magnetism. Here, it is the case of a tree-cotree block-wise decomposition, as the one proposed in [START_REF] Alotto | Mixed finite element methods and tree-cotree implicit condensation[END_REF] for r = 0, to impose the divergence free constraint and generalise it to the case r > 0 of high order approximations. This new tree-cotree approach results in a smaller system to be solved for the discharge unknown.

The paper is organised as follows:

-In Section 2, we formulate Darcy's problem in terms of differential forms. This is a delicate step where the involved physical phenomenon matters to determine the nature of the fields at play. The FEEC approach then guides the discrete formulation. If n denotes the dimension of the physical domain Ω where the Darcy's problem is set and T is a simplicial mesh over Ω, the pair of spaces P - r+1 Λ n-1 (T ), of the Raviart-Thomas H(div)-conforming finite elements of degree r + 1, for the approximation of the discharge, and P - r Λ n (T ), of the discontinuous piecewise polynomial elements of degree r for the L 2 -approximation of the fluid's pressure pop up naturally.

-In Section 3, we write, in vector formalism, the discrete formulation suggested by the FEEC approach of the previous section. We recall the stability and consistency of the FE approach and give its matrix formulation. -In Section 4, we introduce the new degrees of freedom, the weights, defined in [START_REF] Rapetti | Whitney forms of higher degree[END_REF], and the basis functions for the high order mixed finite element spaces involved in the discretisation of Darcy's equations. -The next three sections are dedicated to numerical results. In Section 5, we focus on the numerical approximation of the inf-sup constant. We check, both in h and in r, its convergence to the continuous value for a model problem. The choice of degrees of freedom (and thus of the basis functions in duality with them) does not influence its value but the way how the divergence free constraint is imposed. This latter aspect is discussed in Section 6 where a tree-cotree decomposition is used to address this issue. Section 7 includes some considerations about the computational complexity of the proposed approach, that could be introduced in existing software packages to reduce the cpu time for the solution of the discrete problem. -Finally, in Section 8, we recap the methodology adopted in the paper and draw some conclusions on the basis of the results of the numerical tests.

Darcy's problem governing equations

We recall the governing equations of Darcy's stationary flows and comment on the nature of the involved fields (velocity, pressure, ...) in the FEEC framework.

Let Ω be a bounded Lipschitz polyhedral domain of R n , with n ≥ 1. The experimental law formulated in 1856 by the French hydraulic engineer Henry P. G. Darcy (1803Darcy ( -1858) ) in [START_REF] Darcy | Les fontaines publiques de la ville de Dijon[END_REF] can be written in terms of the instantaneous discharge rate q ([L][T ] -1 ) of the fluid through Ω as

q = - κ ν grad (hg) ,
where κ ([L] 2 ) is the intrinsic permeability of the medium in the domain, ν

([L] 2 [T ] -1
) the kinematic viscosity of the fluid flowing through the domain, g

([L][T ] -2
) is the acceleration due to gravity and h ([L]) the hydraulic head.

In this form it is very similar to Ohm's law for conductors J = -σ grad V where the electric conductivity σ is in fact a Hodge operator from the space of differential 1-forms Λ 1 (Ω) to the space of differential 2-forms Λ 2 (Ω). The hydraulic potential (hg) (resp., the discharge rate q) acts like the voltage V (resp., the current J) and the constant κ/ν is the hydraulic conductivity that describes the fluid/medium as the electric conductivity σ characterises the conductor. We now replace the hydraulic potential by the expression hg = p ρ + zg where p

([M ][L] -1 [T ] -2 ) is the pressure at the chosen point, ρ ([M ][L] -3 )
is the density of the fluid and z ([L]) is the elevation of the point above a reference plane, with the positive z-direction pointing upward (so in the direction opposite to the gravitational acceleration). We thus get

q = - κ ν grad p ρ + κ ν f
where f = -grad(zg) is associated with the additional external potential (zg) due to the elevation z. Equivalently one has

ν κ q = -grad p ρ + f .
The problem's formulation in terms of differential forms, as explained in Fig. 1, makes use of two well known operators: one is the exterior derivative d, necessary to compute derivatives, like gradients, divergences, or curls; the other is the Hodge star * that takes into account the properties of the specific media. The fluid discharge q is a 2-form as it indicates the mass of fluid transported across a given surface area per unit time. The hydraulic potential is a 0-form, a scalar function that can be evaluated at any point of the medium. The key point is that these two forms do not live on the same complex. Once introduced the operators d and * , we can organise the information in a Tonti's type diagram [START_REF] Tonti | On the mathematical structure of a large class of physical theories[END_REF], as the one in Fig. 1. The primal complex (resp. dual complex) is drawn vertically on the left (resp. right) and we move up-down (resp. downup) with the operator d. The two complexes are linked by horizontal lines, one line for each Hodge operator. In Darcy's law, the two forms q and dV are related by the Hodge operator * ν κ . Moving vertically upward, from q, by applying d we get s, namely d q = s. Darcy's problem thus reads

* ν/κ q + d( * 1/ρ p) = f , d q = s . (1)
Note that the d contained in the left-hand equation acts on a 0-form whereas the d appearing in the right-hand equation acts on a (n -1)-form. In R 3 with the standard metric they correspond to the grad operator on scalar fields and the div operator on vector fields respectively. From (1), we have

Ω * ν/κ q + d( * 1/ρ p) ∧ v = Ω f ∧ v , (2)
Fig. 1 Tonti's diagram for Darcy's law in terms of differential forms. On the vertical left and right sides of the diagram, we indicate, for n = 2, the primal and dual proxy complexes, i.e., the two de Rham sequences of functional spaces (see, e.g., their definitions in [START_REF] Bernardi | Discrétisations variationnelles de problèmes aux limites elliptiques[END_REF]), associated with Darcy's law in vector formalism. When n = 2, the curl operator is generally denoted either rot when applied to a vector field

v = (v 1 , v 2 ) T (we have rot v = ∂x 1 v 2 -∂x 2 v 1 )
or rot when applied to a scalar field v (we have

rot v = (∂x 2 v , -∂x 1 v) T ). Ω d q ∧ * 1/ρ ψ = Ω s ∧ * 1/ρ ψ . (3)
Here v and ψ are a (n -1)-form and a n-form respectively. Moreover, since ( * 1/ρ p) is a 0-form, the integration by parts' formula yields

Ω d( * 1/ρ p) ∧ v = - Ω * 1/ρ p ∧ dv + ∂Ω Tr( * 1/ρ p) ∧ Tr v .
Hence equations ( 2) and (3) (for s = 0) read

Ω * ν/κ q ∧ v - Ω * 1/ρ p ∧ dv = Ω f ∧ v - ∂Ω Tr( * 1/ρ p) ∧ Tr v , (4)
Ω d q ∧ * 1/ρ ψ = 0. (5)
The FEEC approach guides the discrete formulation. We decompose Ω into a finite set T of n-simplices. For r ≥ 0, let P - r+1 Λ k (T ), with k ∈ {0, 1, . . . n}, be the spaces of Whitney's differential k-forms [START_REF] Whitney | Geometric integration theory. 400[END_REF] of polynomial degree r + 1. For 0 ≤ k ≤ n, it is possible to connect these spaces of differential forms with finite element spaces of scalar and vector functions [START_REF] Boffi | Mixed finite element methods and applications[END_REF], thank to the pioneering work, done by Bossavit in the early 80s (see [START_REF] Bossavit | Computational electromagnetism: variational formulations, complementarity, edge elements. 375[END_REF] details). When n = 3, the set P - r+1 Λ 1 (T ) and P - r+1 Λ 2 (T ) corresponds, respectively, with the H(curl) and H(div) conforming finite element spaces of degree r +1, belonging to the first Nédélec's family [START_REF] Nédélec | Mixed finite elements in R 3[END_REF]. For n = 2, the set P - r+1 Λ 1 (T ) correspond with the H(div)-conforming finite element space of degree r + 1, firstly introduced in [START_REF] Raviart | A mixed finite element method for 2nd order elliptic problems. Mathematical aspects of finite element methods[END_REF]. The pair of spaces P - r+1 Λ n-1 (T ), of the Raviart-Thomas H(div)-conforming finite elements of degree r + 1, for the approximation of the discharge, and P - r+1 Λ n (T ), of the discontinuous piecewise polynomial elements of degree r for the L 2 -approximation of the fluid's pressure pop up naturally among other possible pairs.

We have naturally obtained the mixed finite element discretisation of Darcy flow following the dual complex. Following the primal complex, Darcy flow can be rewritten in terms of the hydraulic potential and this potential discretised in P - r+1 Λ 0 (T ) that corresponds to using Lagrange finite elements. Both approaches involve only the original (primal) mesh T of Ω. An alternative discretisation adopting also the dual mesh has been proposed in [START_REF] Hirani | Numerical method for Darcy flow derived using discrete exterior calculus[END_REF] using the discrete exterior calculus (DEC).

Whitney finite elements for the Darcy problem

We recall in this section some well-known results concerning the mixed finite element approximation of the Darcy's problem solution (see, e.g., [START_REF] Quarteroni | Numerical approximation of partial differential equations[END_REF]). We consider the following boundary value problem:

       ν κ q + grad[ 1 ρ p] = f in Ω div q = 0 in Ω q • n = 0 on Γ N 1 ρ p = H D on Γ D (6)
where n is the unit outward normal to ∂Ω and Γ N ,

Γ D are open disjoint subsets of ∂Ω such that ∂Ω = Γ N ∪ Γ D . The function f ∈ (L 2 (Ω)) n is known, being L 2
(Ω) the space of measurable, square integrable functions on Ω, together with

H D ∈ H 1/2 (Γ D ).
For the latter, we set H 1 (Ω) = {q ∈ L 2 (Ω), grad q ∈ (L 2 (Ω)) n } and H 1/2 (∂Ω) denotes the image in L 2 (∂Ω) of the trace operator γ acting on functions of H 1 (Ω). The space H 1/2 (Γ D ) contains the restrictions to Γ D of the functions belonging to H 1/2 (∂Ω). Let us also consider the following functional spaces:

H(div; Ω) = {v ∈ (L 2 (Ω)) n : div v ∈ L 2 (Ω)} , H 0,N (div; Ω) = {v ∈ H(div; Ω) : v • n = 0 on Γ N } . Remark 1 If Γ N is a subset of ∂Ω with positive measure, the space H 1/2 00 (Γ N ) consists of those elements ξ ∈ H 1/2 (Γ N ) whose trivial extension ξ 0 of ξ by zero to all ∂Ω belongs to H 1/2 (∂Ω) [26], i.e. H 1/2 00 (Γ N ) = {ξ ∈ H 1/2 (Γ N ) | ξ 0 ∈ H 1/2 (∂Ω)}, ||ξ|| H 1/2 00 (Γ N ) = ||ξ 0 || H 1/2 (∂Ω) . The condition v • n = 0 on Γ N means that v • n, ξ 0 ∂Ω = 0 for all ξ ∈ H 1/2 00 (Γ N )
, where •, • ∂Ω stands for the dual pairing. In the following we will use this notation: if v ∈ H 0,N (div; Ω) then

Γ D H D v • n := v • n, H D ∂Ω being H D ∈ H 1/2 (∂Ω) any extension of H D ∈ H 1/2 (Γ D ).
The variational formulation of (6) reads: find (q, p) ∈ H 0,N (div; Ω)×L 2 (Ω) such that

Ω ν κ q • v -div v 1 ρ p = Ω f • v -Γ D H D v • n Ω div q 1 ρ ψ = 0 , for all (v, ψ) ∈ H 0,N (div; Ω) × L 2 (Ω).
As usual, we consider the bilinear forms a(w, v)

= Ω ν κ w•v and b(v, ψ) := Ω div v 1 ρ ψ defined on H 0,N (div; Ω) × H 0,N (div; Ω) and H 0,N (div; Ω) × L 2 (Ω) respectively. If Γ D = ∅ it is easy to check that a(•, •
) is continuous and coercive on the space V 0 defined as

V 0 = {v ∈ H 0,N (div; Ω) : b(v, ψ) = 0 ∀ ψ ∈ L 2 (Ω)} = {v ∈ H 0,N (div; Ω) : div v = 0} ,
and that the inf-sup condition is satisfied, namely such that

inf ψ∈Q sup v∈V b(v, ψ) v H(div;Ω) ψ L 2 (Ω) =: β > 0 , (7)
being

V = H 0,N (div; Ω) and Q = L 2 (Ω). If Γ D = ∅ then the coercivity and inf-sup condition hold replacing L 2 (Ω) by L 2 0 (Ω) = {ψ ∈ L 2 (Ω) : Ω ψ = 0}.

Remark 2

The mathematical tools introduced in this section are strictly related to those involved in the FEEC formulation of the previous section. The first two equations in ( 6) are the L 2 version of (1), with s = 0. The functional spaces L 2 (Ω), H(div; Ω), H 1 (Ω) introduced above, are those of the Hilbert complex related with the de Rham sequences appearing on Fig. 1. The variational formulation of ( 6) is given by, respectively, equations ( 4) and ( 5) in vector formalism, where ∂Ω = Γ D ∪ Γ N together with Tr( * 1 ρ p) = H D on Γ D and Tr(v) = 0 on Γ N . Obtaining a well-posed numerical scheme for this problem requires to reproduce the de Rham sequence at the discrete level. This is detailed here below.

From now on, we assume that Ω is covered by a simplicial mesh T with boundary faces located on either Γ N or Γ D . For the sake of clarity we will assume that Γ D = ∅. We consider the discretisation of the weak form of (6) involving Raviart-Thomas finite elements RT r+1 of degree r + 1 and discontinuous piecewise polynomial functions DP r of degree r, with r ≥ 0. Let us recall that RT r+1 (T ) = P - r+1 Λ n-1 (T ) and DP r (T ) = P - r+1 Λ n (T ). Let us introduce the space RT 0,N r+1 (T ) = RT r+1 (T ) ∩ H 0,N (div; Ω), so that the discrete problem can be formulated as follows: find (q h , p h ) ∈ RT 0,N r+1 (T ) × DP r (T ) such that

Ω ν κ q h • v h -div v h 1 ρ p h = Ω f • v h -Γ D H D v h • n, Ω div q h 1 ρ ψ h = 0 . ( 8
) for all (v h , ψ h ) ∈ RT 0,N r+1 (T ) × DP r (T ). It is easy to check that a(•, •) is coercive in the space V 0 h = {v h ∈ RT 0,N r+1 (T ) : b(v h , ψ h ) = 0 ∀ ψ h ∈ DP r (T )}. In fact if v h ∈ V 0 h , taking ψ h = div v h ∈ DP r (T ) we have in particular 0 = b(v h , div v h ) = Ω 1 ρ (div v h ) 2 , hence V 0 h = {v h ∈ RT 0,N r+1 (T ) : div v h = 0} ⊂ V 0 .
It is also well-known (see, e.g. [START_REF] Bernardi | Discrétisations variationnelles de problèmes aux limites elliptiques[END_REF]) that the discrete inf-sup condition is satisfied, namely,

inf ψ h ∈Q h sup v h ∈V h b(v h , ψ h ) v h H(div;Ω) ψ h L 2 (Ω) =: β h > 0 , (9)
being Q h = DP r (T ) and V h = RT 0,N r+1 (T ). Moreover this condition holds uniformly with respect to h, i.e., β h ≥ β 0 > 0.

The matrix problem corresponding to the variational formulation ( 8) is of the following type A -B B 0

Q P = F -G 0 . (10)
For given bases {w h,j } N0 j=1 and {π h,i } M i=1 of RT 0,N r+1 (T ) and DP r (T ), respectively (being N 0 = dim RT 0,N r+1 (T ) ≤ dim RT r+1 (T ) = N and M = dim DP r (T )), the entries of the matrices A ∈ R N0×N0 and B ∈ R M×N0 are

A k,j = a(w h,j , w h,k ) , B i,j = b(w h,j , π h,i ) .
It is worth noting that the discrete inf-sup condition ( 9) is equivalent to B being full rank. Similarly the components of the vectors F and G ∈ R N0 are

F k = Ω f • w h,k , G k = Γ D H D w h,k • n .
The components of the vectors Q ∈ R N0 and P ∈ R M are the (unknown) coefficients of q h ∈ RT 0,N r+1 (T ) and p h ∈ DP r (T ) in the basis {w h,j } N0 j=1 and {π h,i } M i=1 , respectively. It is time to build up the bases {w h,j } N0 j=1 and {π h,i } M i=1 for a suitable choice of unisolvent degrees of freedom for the fields in the spaces RT 0,N r+1 (T ) and DP r (T ). This is the goal of the next section.

Weights on small simplices and bases in duality

We are seeking for FE bases in duality with a set of degrees of freedom different from the classical moments, the weights on small simplices introduced in [START_REF] Rapetti | Whitney forms of higher degree[END_REF], that have a direct physical interpretation. These degrees of freedom for a kform w ∈ P - r+1 Λ k (T ) are integrals of w on k-chains. This definition fits, for r > 0, with the meaning of degrees of freedom as cochains, in the spirit of Whitney's original work.

We need to rely on multi-index notation α = (α 0 , . . . , α n ) ∈ Z n+1

n i=0 [( 1 n+1 + α i )v i ]/(r + 1
) and 1/(r + 1)-homothetic to T . The latter expression signifies that {α, T } has the same shape and orientation of T and that its dimensions are reduced by a factor 1/(r + 1) with respect to those of T . For k ∈ {0, 1, . . . , n} let us denote ∆ k (T ) the set of k-subsimplices of T . For any f ∈ ∆ k (T ), the k-simplex {α, f } belongs to the small triangle {α, T }, it is parallel to f and 1/(r + 1)-homothetic to f . For the finite element approximation of Darcy problem we use the spaces RT r+1 (T ) = P - r+1 Λ n-1 (T ) and DP r (T ) = P - r+1 Λ n (T ).

Definition 1 For r ≥ 0, in each n-simplex T ∈ T h , the weights of v h ∈ RT r+1 are the functionals that associate v h with its fluxes through the small (n -1)faces of T , {α, f }, namely

ξ f α : v h → {α,f } v h • n T k with α ∈ I n-1 r and f ∈ ∆ n-1 (T ) , (11)
and the weights of ψ h ∈ DP r are the functionals that associate ψ h with its integrals over the small n-simplices {β, T }, namely

ξ T β : ψ h → {β,T } ψ h with β ∈ I n r . (12)
In general the weights of elements of P - r+1 Λ k (T ) are integrals on the small k-simplices, the so called k-weights. In [START_REF] Christiansen | On high order finite element spaces of differential forms[END_REF] it has been proved that these degrees of freedom are unisolvent. However it is well know that the number of small n-simplices is equal to the dimension of P - r+1 Λ n (T ) but the number of small (n -1)-simplices is bigger than the dimension of P - r+1 Λ n-1 (T ). In [START_REF] Alonso Rodríguez | Minimal sets of unisolvent weights for high order Whitney forms on simplices[END_REF] a minimal set of unisolvent k-weights is identified. In the particular case of k = n -1 a minimal set of (n -1)-weights is that obtained by eliminating those weights associated with the small (n -1)-faces at the interior of any element T ∈ T that are parallel to one of the faces of T , for instance, parallel to the face opposite to the first vertex of T .

The spaces RT r+1 (T ) and DP r (T ) are constructed locally, on one element T ∈ T , and then local contributions are assembled by relying on local-to-global numbering of the degrees of freedom, as it is usual in finite elements. The dimensions of the spaces RT r+1 (T ) and DP r (T ) (see e.g. [START_REF] Quarteroni | Numerical approximation of partial differential equations[END_REF], with

x ∈ R n ) are N loc = dim RT r+1 (T ) = dim(P n r ⊕ x P r ) = (n+r-1)! (n+r+1) (n-1)!r! , M loc = dim DP r (T ) = dim P r = (n+r)! n! r! .
We now introduce a basis {w h,j } Nloc j=1 ⊂ RT r+1 (T ) and a basis {π h,i } Mloc i=1 ⊂ DP r (T ), that are in duality with suitable degrees of freedom for v h|T ∈ RT r+1 (T ) and q h|T ∈ DP r (T ), respectively. We recall that, given a basis {φ j } j of a finite-dimensional vector space P and a set {σ i } i of the dual space P , we say that {φ j } j is in duality with {σ i } i if σ i (φ j ) = δ i,j , being δ .. the Kronecker's symbol. Starting from a generic basis {ψ j } j of P , we can build up the basis {φ j } j in duality with the selected functionals {σ i } i , by solving a linear system with matrix V ij = σ i (ψ j ), known in the literature as generalised Vandermonde matrix. For r > 0, the definition of FE bases {w h,j } Nloc j=1 ⊂ RT r+1 (T ) (resp. {π h,i } Mloc i=1 ⊂ DP r (T )) in duality with suitable unisolvent degrees of freedom {σ i } i involves the multiplication by the inverse of the corresponding generalised Vandermonde matrix, as explained in [START_REF] Bonazzoli | High order finite elements in numerical electromagnetism: degrees of freedom and generators in duality[END_REF]. We now detail the construction of the FE bases in duality with the weights on the small simplices.

The generalised Vandermonde matrix for an element T ∈ T is computed from a known basis of the polynomial space associated with P - r+1 Λ k (T ), for instance, a basis with elements of the form λ α ω f with α ∈ I n r and f ∈ ∆ k (T ), being ω f ∈ P - 1 Λ k (T ) the differential Whitney k-form of polynomial degree 1 associated with the k-subsimplex f (see [START_REF] Rapetti | Whitney forms of higher degree[END_REF]). There is a natural correspondence between these polynomial k-forms and the k-simplices {α, f }.

The set {λ α ω f : α ∈ I n-1 r , f ∈ ∆ n-1 (T )} generates P - r+1 Λ n-1 (T )
but it is clear that it is not a basis because, as mentioned before, the number of small (n-1)-simplices is bigger than the dimension of P - r+1 Λ n-1 (T ). It is well known (see, e.g., [START_REF] Arnold | Finite element exterior calculus, homological techniques, and applications[END_REF]) that to obtain a basis it is enough to eliminate those elements with f = T \ v T 0 , being v T 0 the first vertex of T , namely the vertex of T which carries the smallest number, among the n + 1 local vertices of T , in the global numbering of the mesh nodes, and α 0 = 0. Notice that this selection criterium corresponds to the one adopted to define a minimal set of unisolvent weights. In the following, we will just write v 0 instead of v T 0 . On a generic n-simplex T ∈ T , the generalised Vandermonde matrix for the Raviart-Thomas finite elements is a square matrix of size N loc and has entries

(V RT ij =) V RT (α,f ),(α ,f ) = {α,f } λ α ω f • n f .
Here α, α ∈ I n-1 r , the subsimplices f , f ∈ ∆ n-1 (T ), and α 0 = 0 if f = T \v 0 and α 0 = 0 if f = T \ v 0 . For the discontinuous piecewise polynomial finite elements, the generalised Vandermonde matrix is square, of size M loc and has entries

(V DP m =) V DP (α,T ),(α ,T ) = {α,T } λ α , with α, α ∈ I n r .
The inverse matrix (V RT ) -1 (resp. (V DP ) -1 ) gives, column by column, the coefficients to express the functions {w h,j } Nloc j=1 (resp. {π h, } Mloc =1 ) of the dual basis of RT r+1 (T ) (resp. of DP r (T )) as a linear combination of the basis {y h,j } Nloc j=1 (resp. {z h, } Mloc =1 ) with elements of the form y h,j = λ α ω f (resp.

z h, = λ α ω T ) that is w h,j (x) = Nloc i=1 [(V RT ) -1 ] i,j y h,i (x) π h, (x) = Mloc s=1 [(V DP ) -1 ] s, z h,s (x) .
If Q is the vector with entries the weights of a discrete function q h|T and {w h,j } Nloc j=1 ⊂ RT r+1 (T ) is the basis in duality with the set of weights, one has

q h|T (x) = Nloc j=1 Q j w h,j (x) = Nloc j=1 Q j Nloc i=1 [(V RT ) -1 ] i,j y h,i (x) = Nloc i=1 Nloc j=1 [(V RT ) -1 ] i,j Q j y h,i (x) = Nloc i=1 (Q T • [(V RT ) -1 ] i,• )y h,i (x).
In [START_REF] Bonazzoli | High order finite elements in numerical electromagnetism: degrees of freedom and generators in duality[END_REF], it has been stated that the independence of a degree of freedom from the metric of the subsimplex of T that supports it makes the entries of the Vandermonde matrix, such as V RT or V DP , computable once on a generic triangle T and valid in any other triangle T different from T , up to a suitable renumbering of the indices. The inverse matrix (V RT ) -1 (resp. (V DP ) -1 ) of size N loc N 0 (resp. M loc M ) is thus computed only once for a generic T but it is valid for all T ∈ T , up to a suitable renumbering of the indices. The value q h (x) of the discrete field q h at x ∈ Ω is computed as q h (x) = T ∈T q h|T (x) where stands for assembling the local contributions q h|T (x) associated with the elements T ∈ T .

Numerical evaluation of the inf-sup condition

Before solving the Darcy's problem, we wish to estimate numerically the constant of the inf-sup condition, using the matrices computed with the bases defined in the previous section. This is a first important difficulty to face: we thus check that the new proposed FE approach is stable and that the computed quantity (the inf-sup constant) converges to the analytic one (a known value stated by theoretical results), for h decreasing or r increasing.

Coming back to the linear system [START_REF] Arnold | Finite element exterior calculus: from Hodge theory to numerical stability[END_REF], it is easy to check that A is invertible. It is in fact symmetric definite positive. The vector Q can be eliminated yielding

SP = -BA -1 (F -G) being S = BA -1 B T .
The matrix S is called the Schur complement of A. Let us denote M ∈ R M×M the matrix with entries M i,j = Ω π h,j π h,i , We recall the following result (see the proof in [START_REF] Ern | Finite elements II: Galerkin Approximation[END_REF]).

Proposition 1 Assume that the bilinear form a(•, •) is symmetric continuous with constants C a and coercive on V h with constants α and that the bilinear form b(•, •) is continuous with constants C b and that the inf-sup condition (9) holds. Then the matrices S and M are spectrally equivalent, i.e., the following holds for all

Φ ∈ R M β 2 h C a ≤ Φ T SΦ Φ T M Φ ≤ C 2 b α .
Notice that if we set the physical parameters equal to one then we have C a = C b = 1. We have:

Sup v h ∈RT 0,N r+1 Ω div v h ψ h v h (L 2 (Ω)) n = Sup V∈R N 0 Φ T BV (V T AV) 1/2 (13)
Taking in particular V = A -1 B T Φ and recalling that S = BA -1 B T one has

Sup V∈R N 0 Φ T BV (V T AV) 1/2 ≥ Φ T SΦ (Φ T SΦ) 1/2 = (Φ T SΦ) 1/2 .
Let λ * and λ * be respectively the smallest and the biggest eigenvalues of the following generalised problem: find a scalar λ and a vector Φ = 0 such that

S Φ = λM Φ . (14)
Then

λ * = Inf Φ∈R M Φ T SΦ Φ T M Φ and λ * = Sup Φ∈R M Φ T SΦ Φ T M Φ Hence β 2 h = λ * ≤ Φ T SΦ Φ T M Φ ≤ λ * = 1 α .
So we can estimate β h by computing the smallest eigenvalue of the generalised eigenvalue problem [START_REF] Bossavit | Computational electromagnetism: variational formulations, complementarity, edge elements. 375[END_REF] (see [START_REF] Fortin | Stability Analysis of Discrete Generalized Stokes Problem[END_REF]; see also [START_REF] Quarteroni | Numerical approximation of partial differential equations[END_REF], [START_REF] Boffi | Mixed finite element methods and applications[END_REF]).

Let us consider the infinite-dimensional counterpart. If Γ D = ∅ then for any ψ ∈ L 2 (Ω) there exists a unique

u ψ ∈ H 1 0,Γ D (Ω) := {ζ ∈ H 1 (Ω) : γξ |Γ D = 0} such that Ω grad u ψ • grad ξ = - Ω ψ ξ, ∀ ξ ∈ H 1 0,Γ D (Ω).
Let us denote C Ω the constant in the Poincare's inequality

ξ L 2 (Ω) ≤ C Ω grad ξ (L 2 (Ω)) n , ∀ ξ ∈ H 1 0,Γ D (Ω). ( 15
)
Clearly Fig. 2 Convergence of the inf-sup constant β h to the theoretical value √ 2π (left) and π (right) for h → 0, varying the approximation degree r (see Tables 1 and2, resp., for the computed values of β h ).

grad u ψ 2 (L 2 (Ω)) n ≤ ψ L 2 (Ω) u ψ L 2 (Ω) ≤ ψ L 2 (Ω) C Ω grad u ψ (L 2 (Ω)) n , h = 1/J r + 1 1/4 1/
hence grad u ψ (L 2 (Ω)) n ≤ C Ω ψ L 2 (Ω) . Denoting z ψ = grad u ψ , since we have div z ψ = ψ one gets Sup v∈H 0,N (div;Ω)

Ω div v ψ v (L 2 (Ω)) n ≥ Ω div z ψ ψ z ψ (L 2 (Ω)) n = ψ 2 L 2 (Ω) grad u ψ (L 2 (Ω)) n ≥ 1 C Ω ψ L 2 (Ω) .
Table 1 shows the square of the smallest eigenvalue of ( 14) set in the domain Ω = [0, 1] 2 with Γ D = ∂Ω, varying the polynomial approximation degree r and the mesh size h. It converges to √ 2 π ≈ 4.442883 that is the infinite-dimensional counterpart. In fact, it is known (see, e.g., [START_REF] Kuznetsov | Sharp constants in the Poincaré, Steklov and related inequalities (A survey)[END_REF]) that if 2 the smallest positive eigenvalue of the Dirichlet Laplacian.

Ω = (0, a) × (0, b) ⊂ R 2 and Γ D = ∂Ω then the sharp constant in (15) is C Ω = σ -1/2 D being σ D = (π/a) 2 + (π/b)
Table 2 shows the square of the smallest eigenvalue of ( 14) set in the domain Ω = [0, 1] 2 with Γ D equal to the vertical sides of Ω and Γ N equal to the horizontal ones, varying the polynomial approximation degree r and the mesh size h. Figure 2 shows the convergence of the inf-sup constant β h to the theoretical value π for h → 0, varying the approximation degree r. For r = 0 a rather fine mesh is necessary to have β h close to π. As soon as r > 0, the value β h is a good approximation of π already on a very coarse mesh.

Solving the indefinite linear system

A second difficulty to face with the considered Darcy's problem is the imposition of the divergence free constraint on the fluid discharge. We propose a new technique for solving the linear system [START_REF] Arnold | Finite element exterior calculus: from Hodge theory to numerical stability[END_REF]. It is inspired on the tree-cotree decomposition of the unknowns, presented, for low order Raviart-Thomas FEs, in [START_REF] Alotto | Mixed finite element methods and tree-cotree implicit condensation[END_REF]. According to [START_REF] Alonso Rodríguez | Small trees for high order Whitney elements[END_REF], using the weights as degrees of freedom for fields in the spaces RT r+1 (T ) and DP r (T ), with r ≥ 0, the matrix associated with the divergence operator D is the incidence matrix of a directed graph G * = (V * , E * ). This fact yields a generalisation to the high-order case r > 0 of a procedure analogous to the one proposed in [START_REF] Alotto | Mixed finite element methods and tree-cotree implicit condensation[END_REF] for r = 0. By choosing a spanning tree of this graph, it is possible to identify an invertible square submatrix of the divergence matrix. Then, by reordering the unknowns associated with the arcs of the tree and of the cotree, a solution of the linear system (10) can be found using a simple elimination procedure. The same decomposition can be applied when using the classical moments as degrees of freedom (see [START_REF] Alonso Rodríguez | A tree-cotree splitting for the construction of divergence-free finite elements: the high-order case[END_REF]). Note that moments and weights coincide for r = 0. When r > 0, adopting the weights on the small simplices as degrees of freedom, the tree-cotree decomposition turns out to be more natural.

Remark 3 This is related with the construction of Raviart-Thomas divergencefree basis. The problem for elements of degree one has been studied in different papers (see, e.g., [START_REF] Hecht | Construction d'une base de fonctions P1 non conforme à divergence nulle dans R 3[END_REF], [START_REF] Dubois | Discrete vector potential representation of a divergence-free vector field in three-dimensional domains: numerical analysis of a model problem[END_REF], [START_REF] Scheichl | Decoupling Three-Dimensional Mixed Problems Using Divergence-Free Finite Elements[END_REF], [START_REF] Rapetti | Discrete vector potentials for non-simply connected three-dimensional domains[END_REF], [START_REF] Alonso Rodríguez | Graphs, spanning trees and divergence-free finite elements in domains of general topology[END_REF]). For Raviart-Thomas finite elements of higher degree, the problem has been studied in [START_REF] Alonso Rodríguez | A graph approach for the construction of high order divergence-free Raviart-Thomas finite elements[END_REF] and [START_REF] Alonso Rodríguez | A tree-cotree splitting for the construction of divergence-free finite elements: the high-order case[END_REF] using moments as degrees of freedom. Here, we test high order Raviart-Thomas with weights as degrees of freedom, which makes the construction of the matrix of the divergence operator quite simple.

The graph G * = (V * , E * ) has a node for each small n-simplex plus an additional node corresponding to the exterior of the domain, and an arc for each small (n -1)-simplex associated to the minimal set of (n -1)-weights. In Fig. 3, we show the arcs of the graph G * for n = 2 on a partial configuration with three elements [0, 1, 2], [0, 2, 3], [START_REF] Ainsworth | A Bernstein-Bézier basis for arbitrary order Raviart-Thomas finite elements[END_REF][START_REF] Alonso Rodríguez | Minimal sets of unisolvent weights for high order Whitney forms on simplices[END_REF][START_REF] Alonso Rodríguez | A tree-cotree splitting for the construction of divergence-free finite elements: the high-order case[END_REF] selected in the mesh, each of them having one side on ∂Ω, and an example of spanning tree. When using a basis in duality with the weights, the divergence operator D is the incidence matrix of G * with root the node v ext corresponding to the exterior of the domain.

It is well-known that the submatrix of an incidence matrix of a connected directed graph corresponding to a spanning tree is invertible [START_REF] Thulasiraman | Graphs. Theory and algorithms[END_REF]. In particular if one decomposes the divergence matrix as D = [D st , D ct ] being D st the submatrix of D with columns corresponding to edges in F * (the spanning tree), and D ct the submatrix with the columns left (the cotree), then D st is invertible. This can be extended also to the matrix B of system (10) which is defined as

B i,j = b(w h,j , π h,i ).
Notice that the set {w h,j } N0 j=1 does not include the basis functions corresponding to the small edges on Γ N . Using the same decomposition as for D, one gets

B = [B st , B ct ]. (16)
Lemma 1 B st is invertible.

Proof Let V st be such that B st V st = 0. Then, taking

V = V st 0 ,
follows that BV = 0. This means that div v h = 0 being v h the element of RT 0,N r+1 with weights given by the coefficients of V. Then DV = 0 and in particular, D st V st = 0. Since D st is non-singular it follows that V st = 0. Thus B st is non-singular.

Similarly, one can decompose the matrix A as follows

A = A st,st A st,ct A ct,st A ct,ct .
Thus, we have an indefinite linear system equivalent to [START_REF] Arnold | Finite element exterior calculus: from Hodge theory to numerical stability[END_REF], namely

  A st,st A st,ct -B st A ct,st A ct,ct -B ct B st B ct 0     Q st Q ct P   =   F st -G st F ct -G ct 0   . (17)
It remains only to eliminate Q st and P in the following way

Q st = -B -1 st B ct Q ct , (18)
P = -B - st (F st -G st -A st,ct Q ct + A st,st B -1 st B ct Q ct ) . (19)
Then, the system to solve is

ÃQ ct = F , (20)
where

à = A ct,ct -A ct,st B -1 st B ct + B ct B - st [-A st,ct + A st,st B -1 st B ct ] (21) F = F ct -G ct + B - st (F st -G st ). ( 22
)
Note that if A st,ct = A ct,st , then à is symmetric since we can write

à = A ct,ct + (B -1 st B ct ) A st,st (B -1 st B ct ) -A ct,st B -1 st B ct -(A st,ct B -1 st B ct ) .

Numerical results

Numerical tests are done for the Darcy's problem [START_REF] Alonso Rodríguez | Small trees for high order Whitney elements[END_REF] in the particular case of a two dimensional domain and a homogeneous medium with ν/κ = 1 and ρ = 1. We have considered a square domain Ω = [0, 1] 2 and assumed that the flow across the bottom and top of Ω, is zero, namely, q • n = 0 on Γ N = (0, 1) × {0, 1}. (Note that in the two-dimensional problem, x 2 is the elevation.)

On

Γ D = {0, 1} × [0, 1] we imposed p = H D being H D (x) = 1 if x 1 = 0 0 if x 1 = 1
In order to check the convergence rate we construct a synthetic problem with pressure p ex = 1 -x 1 (that clearly verifies the boundary condition on Γ D ), and discharge (see a visualisation in Fig. 4)

q ex = π 4 x 2 1 x 2 sin π 2 x 2 -1 2 x 2 1 cos π 2 x 2 x 1 x 2 cos π 2 x 2 (23)
that verifies the boundary conditions on Γ N and is divergence-free.

For the finite element approximation of the discharge and the pressure we use the space RT 0,N r+1 (T ) and DP r (T ) respectively with r = {0, 1, 2, 3} over successive uniformly refined oriented meshes T of the domain Ω. The meshes T are constructed by dividing Ω into J 2 , J = 2, . . . , 10, rectangles and then, each of them, in two triangles, to finally get 2J 2 (big) triangles. We use bases in duality with the weights on the small simplices that are associated with the principal lattice of order r in each element of the mesh.

Numerical experiments consist in computing the discrete solution solving the linear systems [START_REF] Arnold | Finite element exterior calculus: from Hodge theory to numerical stability[END_REF] and (20) by using the direct solver of Matlab. The choice of a direct solver is done to avoid the CPU uncertainty related to fixing precision thresholds or choosing preconditioning techniques for iterative algorithms.

In Table 3 we report the approximation error for both the discharge and the pressure in terms of the mesh size h and of the polynomial degree r+1. It is the L ∞ -norm of the fields q ex -q h and p ex -p h , evaluated at the barycentre of the small triangles. Let us recall that q h = N0 j=1 Q j w h,j p h = M j=1 P j π h,j . It is worth noting that for a fixed number of degrees of freedom, one gets smaller error using high order methods. Figure 5 shows graphically these data. It is easy to check that the convergence rate of the discharge in terms of h is equal to the degree of the polynomials in the Raviart-Thomas space while the convergence rate of the pressure is even better than the expected one.

Using the tree-cotree approach, it is possible to decrease significantly the dimension of the system to be solved. In the case of a two dimensional domain with connected boundary and Γ N = ∅ the dimension of the space of discontinuous piecewise polynomial functions of degree r is M = (r + 2)(r + 1)n T /2 3 and the dimension of the space of Raviart-Thomas finite elements of degree r + 1 is N = N 0 = (r + 1) n F + r(r + 1) n T being n T (resp., n F ) the total number of triangles (resp., of faces, that are edges in the two dimensional case) of the mesh. The indefinite linear system (10) has dimension N 0 + M, whereas the matrix à in system [START_REF] Ern | Finite elements II: Galerkin Approximation[END_REF] has the dimension of the cotree, namely, N 0 -M. However the computation of à involves the inverse of the matrix B st that has the dimension of the tree M.

The tree-cotree decomposition has been implemented in Maltab using the commands graph, dfsearch and bfsearch. The last two commands define a tree on a graph using a depth-first search (DFS) or a breadth-first search (BFS) respectively. For the study conducted in this work, the two searches are equivalent and have similar computational times.

Table 4 shows the computational time (in millisecond) for the two spanning trees when considering, in the top part, a fixed approximation order r + 1 = 3 on different meshes, in the bottom part, a fixed mesh size h = 1/4 with solve the discrete problem when relying on the new tree-cotree decomposition. We have R = (1 -t r /t c ), with t r (resp. t c ) the cpu time to solve the discrete problem with the new method (resp. with the classical method). By using the values presented in Table 4 for the case r + 1 = 3 and h ≤ 1/5 with DFS, we have R increasing from 0.2 up to 0.7 (thus a cpu time reduction of 70%).

In Table 5 we report the computational time (in millisecond) for the two spanning trees when considering a number of degrees of freedom similar to the ones associated with the values in boxes on Table 3. For the tree-cotree decomposition the computational time is divided into the time to read the mesh and compute the spanning-tree (ST), the time to compute B -1 st (IB st ) and the time to solve the linear system with matrix à (I Ã).

Conclusions

The term finite element exterior calculus refers to the finite element approximation of differential k-forms. In this framework the elements of the finite element spaces are studied as polynomial differential k-forms. It is then natural to use, as degrees of freedom, integrals on k-chains because, for differential k-forms, they have a clear physical meaning. The set of degrees of freedom considered in this work, the weights, are based on this idea. Using the weights many procedures and techniques that are well understood in the case of finite elements of minimal degree can be extended in a natural way to finite elements of higher degree. This advantage of the weights has been exploited in electromagnetic problems but, from our knowledge, not yet in fluid dynamics. In this work, taking advantage of the use of the weights, we implement a treecotree decomposition for the solution of the linear system resulting from the high-order Whitney approximation of Darcys flows.

We have presented the Darcy's problem in terms of differential forms in the FEEC framework, in order to show how naturally the continuous and discrete functional spaces can be selected when going back to vector and scalar fields. The two physical unknown quantities, the discharge and the pressure, are approximated using high order Whitney differential (n -1)-and n-forms, that correspond to Raviart-Thomas finite elements of degree r+1 and discontinuous piecewise polynomial finite elements of degree r, respectively. Then the new degrees of freedom, the weights on the small simplices, have been considered. We have numerically estimated the inf-sup constant and shown that it is independent from h and r, thus confirming the stability of the proposed method. Successively, to solve the discrete problem taking into account the divergence free constraint on the discharge, we have relied on a tree-cotree block-wise decomposition of the unknowns based on graph theory that is very natural when considering the weights on the small simplices as degrees of freedom. This because, being the weights of an (n -1)-form integrals on (n -1)-faces of n-simplices, the divergence matrix is the transpose of a boundary matrix also when r > 0. This property is well-known and exploited in the case r = 0

The proposed approach was meant for three-dimensional applications [START_REF] Alonso Rodríguez | A graph approach for the construction of high order divergence-free Raviart-Thomas finite elements[END_REF] (see for example [START_REF] Alonso Rodríguez | Small trees for high order Whitney elements[END_REF], Fig. 1.4, for the construction of the dual graph in three dimensions). However for simplicity we have implemented it in Matlab in the two-dimensional case. Some numerical tests validate the interest in the proposed approach. When the computational mesh becomes fine, the use of the tree-cotree decomposition reduces significantly the computational time to solve the linear system. This method could be adopted in any FEM package that implements high order mixed finite elements.

+

  and we define |α| = n k=0 α k . Then I n r = {α ∈ Z n+1 + : |α| = r} . For any α ∈ I n r the small simplex {α, T } is the n-simplex with barycentre at the point of coordinates

Fig. 3

 3 Fig.3An example of a primal graph (left) and of a spanning tree (right) on the small mesh of three elements [0, 1, 2],[START_REF] Ainsworth | A Bernstein-Bézier basis for arbitrary order Raviart-Thomas finite elements[END_REF][START_REF] Alonso Rodríguez | Minimal sets of unisolvent weights for high order Whitney forms on simplices[END_REF][START_REF] Alonso Rodríguez | A tree-cotree splitting for the construction of divergence-free finite elements: the high-order case[END_REF], [0, 2, 3] of T h , with edges [0, 1],[START_REF] Ainsworth | A Bernstein-Bézier basis for arbitrary order Raviart-Thomas finite elements[END_REF][START_REF] Alonso Rodríguez | A tree-cotree splitting for the construction of divergence-free finite elements: the high-order case[END_REF], [0, 3] ∈ ∂Ω, whereas[START_REF] Ainsworth | A Bernstein-Bézier basis for arbitrary order Raviart-Thomas finite elements[END_REF][START_REF] Alonso Rodríguez | Minimal sets of unisolvent weights for high order Whitney forms on simplices[END_REF] is also on the boundary of[START_REF] Ainsworth | A Bernstein-Bézier basis for arbitrary order Raviart-Thomas finite elements[END_REF][START_REF] Alonso Rodríguez | Minimal sets of unisolvent weights for high order Whitney forms on simplices[END_REF][START_REF] Alonso Rodríguez | A tree-cotree splitting for the construction of divergence-free finite elements: the high-order case[END_REF]] and [0, 2] on that of [0, 2, 3]. The doted arcs are not on the graph because they are associated to weights that are not in the minimal set.

Fig. 4

 4 Fig. 4 Computed pressure (left) and discharge (right) for the considered synthetic Darcy's problem in Ω = [0, 1] 2 with physical parameters ν/k = 1, ρ = 1.

Fig. 5

 5 Fig.5Approximation error in terms of the mesh size h for different values of r (left) and in terms of the polynomial degree for differente values of h (right). On the top the discharge, on the bottom the pressure. Marked values correspond to the values in boxes in Table3

Table 1

 1 Inf-sup constant β h varying the polynomial approximation degree r and the mesh size h over the domain Ω = [0, 1] 2 with Γ D = ∂Ω.

			6	1/8	1/10	1/20	1/30
	1	4.478674 4.459351 4.452257 4.448915 4.444339 4.443977
	2	4.444860 4.443288 4.443012 4.442936 4.442936 4.442936
	3	4.442914 4.442886 4.442883 4.442883 4.442883 4.442883
	4	4.442883 4.442883 4.442883 4.442883 4.442883 4.442883
					h = 1/J
	r + 1	1/4	1/6	1/8	1/10	1/20	1/30
	1	3.114585 3.129624 3.134863 3.137286 3.140516 3.141114
	2	3.141637 3.141601 3.141595 3.141594 3.141593 3.141593
	3	3.141593 3.141593 3.141593 3.141593 3.141593 3.141593
	4	3.141593 3.141593 3.141593 3.141593 3.141593 3.141593

Table 2

 2 Inf-sup constant β h varying the polynomial approximation degree r and the mesh size h over the domain Ω = [0, 1] 2 with Γ D = ∅ and Γ N = ∅.

		4.49				r=1		3.15
						r=2
		4.48				r=3 r=4
								3.14
	inf-sup constant	4.45 4.46 4.47					inf-sup constant	3.13
								3.12
		4.44										r=1 r=2
												r=3
		4.43						3.11				r=4
		4	6	8	10	20	30	4	6	8	10	20	30
					1/h						1/h

Table 3

 3 Approximation error for the discharge (top) and the pressure (bottom) varying the approximation order r + 1 and the mesh size h. Values in boxes correspond to problems with a total number of degrees of freedom close to 1000. It is worth noting that for a fixed number of degrees of freedom, one gets a smaller error using high order methods.

			1/6	1/8	1/10
	1	1.582e-01	1.115e-01	8.469e-02	6.797e-02
	2	1.169e-02	5.295e-03	3.006e-03	1.934e-03
	3	7.101e-04	2.395e-04	1.090e-04	5.813e-05
	4	3.043e-05	6.086e-06	1.938e-06	7.985e-07
			h = 1/J
	r + 1	1/4	1/6	1/8	1/10
	1	6.707e-03	3.706e-03	2.269e-03	1.515e-03
	2	3.132e-04	9.632e-05	4.302e-05	2.278e-05
	3	4.164e-06	8.987e-07	2.980e-07	1.266e-07
	4	3.713e-07	5.189e-08	1.332e-08	6.535e-09

Acknowledgements This research was supported by the program MathIT financed by the ANR-15-IDEX-01 of the Universite Côte Azur (UCA) in Nice, France. E. Zappon warmly thanks the Università degli Studi di Trento (Italy) for the possibility of studying as ERAS-MUS fellow at the Université Côte Azur, where this work began. The first author was partially supported by the project PRIN 201752HKH8.

The authors wish to dedicate this work to the memory of Christine Bernardi, research director at the CNRS. Numerical methods applied to Darcy's model was one of her favorite subjects.

 4Computational time for the direct solver and tree-cotree method, when considering either the same approximation order r + 1 = 3 on different meshes (top part) or the same mesh size h = 1/4 and different approximation orders. In the second column it is reported the number of degrees of freedom of each test case and, in the third column the time (in millisecond) to solve the linear system with matrix M using the Matlab backslash command.

In the last two columns, we indicate the time (in millisecond) to solve the linear system relying on the tree-cotree decomposition with the depth-first and the breadth-first searches.

In these two last columns, we also report the reduction index (R). Table 5 Computational time for the direct solver and tree-cotree method, when considering a given number of degrees of freedom (see the values in boxes of Table 3). The approximation order together with the mesh size and the number of degrees of freedom of each test case are reported in the first and second columns, respectively. In the third column the time (in millisecond) to solve using the Matlab backslash command the linear system with matrix M . The last eight columns concern the tree-cotree decomposition with the depth-first search and the breadth-first search. ST indicates the time to read the mesh and compute the spanning tree, IBst indicates the time to invert the matrix Bst and I Ã indicates the time to invert Ã. Tot indicates the total time, namely, the sum of the values ST, IBst and I Ã contained in the three previous columns.

Direct

different approximation orders. It also includes the computational time to solve the complete linear system with a direct method. The computational time is similar for the two choices of spanning tree. The BFS tree seems to be slightly more efficient. The solution of the discrete problem based on the presented tree-cotree decomposition is clearly faster than the direct solution of the linear system when the total number of degrees of freedom increases. One can indeed estimate the index of reduction, say R, on the cpu time necessary to


Christina Katsamaki
email: christina.katsamaki@inria.fr

Fabrice Rouillier
email: fabrice.rouillier@inria.fr

Elias Tsigaridas
email: elias.tsigaridas@inria.fr

PTOPO: A Maple package for the topology of parametric curves

published or not. The documents may come

Topology of Parametric Curves

The study of parametric curves is a classical topic in computational algebra and geometry Sendra and [START_REF] Sendra | Algorithms for Rational Real Algebraic Curves[END_REF]; Boissonnat and Teillaud (2006); [START_REF] Busé | Implicitizing rational curves by the method of moving quadrics[END_REF]. The interest for computing with parametric curves has been motivated, among others, by the omnipresence of parametric representations in computer modeling and computer aided geometric design, e.g., [START_REF] Rida | Geometric Design Using Space Curves with Rational Rotation-Minimizing Frames[END_REF].

PTOPO allows the computation of the topology of a real parametric plane curve in maple. In particular, it computes an abstract graph that is isotopic (Boissonnat and Teillaud, 2006, p. 184) to the curve in the embedding space.

Let C be an algebraic curve over C 2 , parametrized by the map

φ : C C t → φ 1 (t), φ 2 (t) = p 1 (t) q 1 (t) , p 2 (t) q 2 (t) , (1)
where p i , q i ∈ Z[t] are of size (d, τ ) for i = 1, 2 , and C is the Zariski closure of Im(φ). We call φ(t) a parametrization of C. We study the real trace of C, that is

C := C ∩ R 2 .
The theoretical background of the PTOPO package, together with a detailed complexity analysis, is given in [START_REF] Katsamaki | On the geometry and the topology of parametric curves[END_REF].

PTOPO computes all points necessary for representing the geometry of the curve, as well as for producing a certified visualization of plane curves. We call them special points and these are the cusps, the multiple, the extreme and the isolated points. On top of that, PTOPO constructs a planar graph whose vertices correspond to points on the curve connected accordingly. The drawing of the graph offers a visualization of the curve in an R 2 -area that contains all its special points.

The computation of the parameters of special points is done by solving over R a bivariate system of polynomials and univariate real solving. Using SLV ? and RootFinding[Isolate], we have a certified implementation of general purpose exact computations with one and two real algebraic numbers, like comparison and sign evaluations. For a parametrization φ, we consider the following system of bivariate polynomials:

h i (s, t) = p i (s)q i (t) -q i (s)p i (t) s -t , for i = 1, 2. (2)
We partition the roots of the bivariate system according to the kind of special point they correspond to, denoting the sets by T C (cusps), T M (multiple points), T E (extreme points), T I (isolated points) and T R P (real poles). Then, we compute the set of all intersections of the curve with a box containing all special points, denoted by T B . Within this box, called characteristic box, the topology of the curve is fully described.

Let T be the union of the above sets of parameters. Add a vertex v i for every distinct point corresponding to a parameter in T and let λ(v i ) be its label set (containing the parameter values giving the point). If for two consecutive elements t 1 < t 2 in T there exists a pole s ∈ T R P such that t 1 < s < t 2 , then we split T into two lists: T 1 containing the elements ≤ t 1 and T 2 containing the elements ≥ t 2 . Continue recursively for T 1 and T 2 , until there are no poles between any two elements of the resulting list. This procedure partitions T into T 1 , . . . , T .

Consider each T i with more than one element separately. For any two consecutive elements t 1 < t 2 in T i , with t 1 ∈ λ(v i,1 ) and t 2 ∈ λ(v i,2 ), we add the edge {v i,1 , v i,2 }. If p ∞ exists, we add an edge to the graph connecting the vertices corresponding to the last element of T and the first element of the T 1 .

For plane curves, if N = max{d, τ }, the complexity of computing the topology becomes O B (N 6 ). However, visualizing the curve on top of the abstract graph construction, increases the bound to O B (N 7 ). We tested PTOPO for many parametric curves from the literature. In practice, for all examples, PTOPO computes the topology in less than a second.

The PTOPO package can be obtained from https://webusers.imj-prg.fr/ ~christina.katsamaki/ ptopo/. PTOPO requires the SLV package.

> r e s t a r t ; LIBPATH := " . . . t h e path t o ptopo . . . " : r e a d ( c a t (LIBPATH, "SLV . mpl " ) ) ; r e a d ( c a t (LIBPATH, " ptopo . mpl " ) ) ; > with (PTOPO) ;

The easiest way to see the topology of a curve using PTOPO is to call the draw function with arguments the parametrization of the curve φ = p 1 q 1 , p 2 q 2 . The draw function can be called in two ways. It can take as arguments a list of numerators and a list of denominators.

> p := [ t ˆ2 + 1 , 1 ] ; q := [ t ˆ4 + 1 , t ˆ3 ] ; PTOPO:-draw ( p , q ) ; Alternatively, the argument can be a list of rational functions.

> p h i := normal ( [ s e q ( p [ i ] / q [ i ] , i = 1 . . 2 ) ] ) ; PTOPO:-draw ( p h i ) ;
It is possible to use different variable names, the defaults are 't' and 's'. The number of points to be drawn within the smooth arcs is controlled by the last argument. If the argument is n ,then 2n + 1 points are computed. Increasing n increases the computational cost, but improves the drawing. Note that it does not affect the topological correctness of the drawing.

> PTOPO:-draw ( p , q , ' t ' , ' s ' , 5 ) ;

Although the draw command is enough, if a picture of the topology is needed, PTOPO provides commands for obtaining the topologically interesting points separately. In that case, we initialize by giving a parametrization, and just as like for the draw command, the parametrization can be given as a pair of lists, the first containing the numerators and the second the denominators, or as a list of rational functions.

> p := [ t ˆ2 + 1 , 1 ] ; q := [ t ˆ4 + 1 , t ˆ3 ] ; PTOPO:-p a r a m e t r i z a t i o n ( p , q ) ; > p h i := [ ( 3 * t ˆ4+ 4 * t ˆ3 + 32 * t ˆ2 + 28 * t + 9 9 ) / ( t ˆ2+ t +7) * ( t ˆ2 + 1 ) , ( ( t ˆ2 + t + 7 ) ˆ3 ) / ( ( t +6) * ( t ˆ2 + 1 ) ˆ2 ) ] ; PTOPO:-p a r a m e t r i z a t i o n ( p h i ) ;

There is a point at infinity. We reparametrize....

----****----The parametrization is proper ----****----...

To compute the topology of the given parametric curve, we call The output contains points in 2D (except for poles) and the corresponding parameter values (when applicable). The labels are i for isolated points, p for poles, e for extreme points, d for double points, c for cusps, b for boundary points.

Alternatively, can get separately isolated points, point at infinity, poles, extreme points, double points and cusps. Moreover, we can see points intersecting the boundary of the box containing everything topologically interesting.

Conclusion

The PTOPO package is practical and easy to use. It computes the topology of a curve and can also visualize it. It is fast enough to make its use practical for most applications.

Since the algorithmic foundations for plane and space curves are not different, we plan to extend the implementation to handle space curves in the future.

  After the topology is computed we can see a summary of the points of interest.

	> PTOPO:-t o p o l o g y ( ) ;
	> PTOPO:-point summary ( ) ;
	* Isolated points:
	: [ -0.618034, -0.485868]	i
	: [ -0.618034, +0.485868]	i
	* Point at infinity:
	: [ +0.020816, +0.002915]
	* Poles:
	: t: -0.142857	p
	* Extreme points:
	: [ +1.207107, +3.751142],	t: -0.866081	e
	: [ +1.207107, -3.751144],	t: +0.458575	e
	* Boundary points:
	: [ +1.201470, +5.000000],	t: -0.793994	b
	: [ +1.201470, -5.000000],	t: +0.407871	b
	* Cusps:
	: [ +0.000000, +0.000000],	t: +7.000000	c
	* Double points: