Philippe Abrahami

Michel • Cécile

Publié

Evosheep -Cuneiform Data on Sheep Husbandry in Mesopotamia from the third millennium BC to the second millennium BC

). The research focuses on the region of sheep domestication in Southwest Asia, specifically the Near East, to understand the onset of differentiation of early domestic sheep into multiple breeds. Part of this research is based on the study of cuneiform texts.

The sheep plays a central role in the economy of ancient Mesopotamia through the exploitation of its fleece. The majority of the textiles produced in the Ancient Near East are made of wool. Sheep were the main provider of meat and milk. Other by-products such as leather, fat and sinews, also had an important role in crafts. Sheep organs, in particular the liver, were also used in divination (Postgate 2009(Postgate -2011)). As part of the Evosheep project, we looked at the data in cuneiform texts relating to the different varieties of sheep, where they came from and how they were reared.

Cuneiform sources from the third and early second millennia provide much evidence on the varieties of sheep and sheep husbandry. They originate from major cities of southern Mesopotamia such as Ur, Uruk, Umma, Girsu, and Isin. Important data also come from the Middle Euphrates and Central Anatolia. These texts, written in Sumerian and Akkadian languages, include mainly letters, contracts and administrative documents produced by palaces, temples and family households. Lexical lists, especially from Nippur, provide a whole categorization of sheep according to their origins, age, gender, colour, as well as practical vocabulary related to their exploitation and the management of flocks OB Nippur ura 3. (Fig. 1 Old Babylonian Nippur ura3 lexical list CDLI P228711). Varieties of sheep are designated according to their physical characteristics, especially the tail and the fleece, or to their geographical origins (Steinkeller 1995: 51-54). Several words describe the shape and length of the tail. The most common is the fat-tail sheep (udu gukkal, gukkalum)(Fig. 2 Fat-tail sheep Awassi sheep) These animals are well-attested already during the 26 th century BCE, and they are not associated with a specific colour. Large flocks of this type of sheep are mentioned in an Old Babylonian letter from Mari (18 th century BCE) (ARM 28, 179: 18-21):

Take from your people one thousand fat-tailed sheep, (that is) 500 male fat-tailed sheep and 500 fat-tailed ewes, and bring them to me. The fleece is the second important key feature which distinguishes the different breeds. The most frequent is the aslum sheep (Ur III), also referred to as 'sheep with long wool fibres' (udu siki-sud/gid2). The 'woolly sheep' (udu-siki) appears already in the Early Dynastic texts (26 th cent.) in reference to sheep producing greater quantity or quality of wool, but these 'woolly sheep' do not form a special breed. A text from Girsu suggests that they were selected from a herd of ordinary sheep during an inspection (Nik 1, 155, Selz 1989: 377-378) Another feature that characterises breed varieties is the colour of the sheep fleece. The most numerous evidence relates to black sheep (udu ge6). Sheep qualified as white (udu babbar) are less common presumably because sheep are by default white, so they may be qualified as such by opposition to black ones. The great majority of black sheep originate from western Iran (Šimaški) in Ur III administrative texts (HSS 68, 15: 8-9). In a few instances, sheep are said to be spotted (gun3-a). Their present-day name, 'Jacob sheep' (Fig. 5 Spotted sheep Jacob sheep), originates from Genesis 30: 25-43, which explains the breeding method used by Jacob to obtain his wages from the flock of Laban, which consisted of black, spotted, speckled and striped sheep and goats. Wild sheep are regularly attested as breeding animals, the most commonly used to create hybrids being the 'mountain sheep' (udu hur-sag). This animal is often mentioned together with other wild animals. The 'foreign sheep' (udu kur-ra) is also considered as a wild variety of sheep even though it could be raised, and there were shepherds in charge of 'foreign sheep'. According to literary texts, the wild ram (atudum) lived in the mountain or in the forest, but it also occurs in some flocks in Central Anatolia. The mouflon (tišānum) was considered precious enough to be offered as a diplomatic gift on the occasion of a royal marriage (ARMT 26/1, 11:24) .

Other varieties of sheep were distinguished according to their geographical origin, whether a country or a town. Sumerian sheep (udu eme-gi7) are opposed to 'foreign sheep'. Sheep could come from different parts of Iran, named broadly as Elamite sheep (udu elam). Among these are those from Šimaški, north-west of Iran, in the Zagros moutains. Less numerous are the sheep from DuhduhNI, Sabum, Anšan, and Hupum.

In the Ur III text, sheep may come from the Lullubu, which is located in the Northern Zagros. There are also Šarumiyum sheep, presumably from the Diyala; most of these are recorded as breeders. According to the early II nd millennium texts, several sheep are from Sutû in north Babylonia, from Amurru west of Mesopotamia, or from Šubaru, north-east of Mesopotamia. It is not always clear whether or not some of these geographical characterizations are to be understood as specific breeds.

The terminology of sheep as regards to age and sex of the animal is especially welldeveloped. The generic Sumerian term udu (immerum) includes both male and female sheep. Male sheep (udu-nita2) are supposedly castrated animals because they appear in great numbers in herd compositions. The ewe (u8 immertum) can be pregnant (erītu).

The adult animal (udu-gal) is distinguished from the old one (udu libir). Regarding the young sheep, the terminology is very detailed. There are two expressions to refer to the newborn animal: the newborn sheep (udu gibil) and the just born lamb (sila4 u3-tu-da), then comes the suckling lamb, male or female (sila4/kir11 ga), the half-weaned lamb (sila4/kir11 gab-ba), the lamb of the spring (sila4-nim), the age of the animal may be given in the expression 'a sheep of x years old' (udu mu-x), or according to the number of times the animal was plucked: A 'lamb in age to be plucked' (sila4-ur4), or a 'lamb plucked for the first time' (sila4-1-ur4 ) (Heimpel 1993: 134;Van de Mieroop 1993: 162-165).

There are very large numbers in the overall figures attested in texts already from the early III rd millennium, some of which exceed ten thousand animals. Numbers are even higher during the Ur III period with, for example, a total of 347,394 sheep among a large variety of other domestic and wild animals (RA 063, 102:9). According to Maeda (1989, 84), some 77,000 sheep and goats would have been transferred through the royal animal park of Puzriš-Dagan every year.

The sheep and goats belonging to the great institutions were divided under the authority of officials (na-gada / nāqidum) who themselves entrusted smaller numbers to individual shepherds (Greco 2021.2: 8-9). For example, an administrative document records 12,867 sheep entrusted to more than twenty shepherds (BM 024954).

Ur III 'shepherd texts' from Lagash present a fixed list of sheep and goats entries used to record the different categories of animals composing the flocks (Mander 2008: 99-108). The head ram bucks (maš2 sag) are few and usually placed at the head of the flock because of their intelligence, so they could lead it.

During the early II nd millennium, shepherds are often responsible for some 200 animals or less, the number of sheep being five to six times more important than the number of goats, perhaps because sheep by-products were much more valuable. The herds often mix different sheep breeds and age as in the following example from Mari (ARM 23, 219: 1-7): 65 amorite sheep, 16 sheep grass fattened, 4 long-fleece Amorite sheep, 2 long-fleece grass fattened sheep, 6 spring male lamb, 6 spring female ? lamb, total: 99 sheep.

The male sheep were normally used for divination or for consumption. Those selected for reproduction were kept in the herd and defined as geš-du3. All types of sheep could be used for that purpose and defined as such, whether coming for different regions or with tail or fleece characteristics. The text PDT 2, 1052: iii 18-25) mentions several of these breeders. Selection for breeding applies also to young animals as it the case with lambs in BIN 3,359: 4 In Ur III texts, the hybrid sheep is referred to as literally 'sheep (of) mountain sheep sperm' (udu a-udu hur-sag). By contrast to the sheep breeder (udu geš-du3), there are no other types of sheep used for hybridization besides the 'moutain sheep'. The female, never specified, is probably a domestic ewe. This use of mountain sheep could be an evidence for the practice of 're-breeding' which involves occasional use of a wild male (here a mountain sheep) with domestic females in order to avoid inbreeding, or to introduce a particular wild character to the domestic.

An interesting point emerging from Ur III documentation is that the resulting hybrids are fertile. As a matter of fact, some sheep are recorded as 'hybrid breeder sheep' (udu a-udu hur-sag geš-du3, see for instance PDT 2, 1052: iii 24).

Hybrid sheep are considered as a specific breed among others (long-tail, long fleece, fat-tail, and sheep from different geographical origins), as shown through the descriptive categories which composed some lists of flocks (PDT 1, 345).

Sheep were grazed in pasture (udu u2) or fattened (udu niga) with various foodstuffs, starting with barley. According to their age, sheep received different daily rations of barley, quantities varying around a litre or more. They were sometimes given complements such as bran (duh) (Nisaba 31/2, 86), and even flour (zi3-kum) (CUSAS 15, 28). The fattener (kurušda) was in charge of feeding the sheep.

Royal archives from Mari provide interesting data in relation to movement of herds toward pasture. This could concern a very wide spatial range, stretching from the Hamrim Moutains (northeast Iraq, ancient Mount Ebih) to a western branch of the Habur, as noted in the following letter: (Fig. 6 Range of the grazing flocks according to FM 6, 9). "News of the grazing flocks (nawû) does not reach me regularly!" This is what my Lord wrote to me. The grazing flocks are scattered throughout the land. On the right of the pasture is Mount Ebih and on its left is Talhayum' (FM 6,9). The sheep were entrusted to a professional, the shepherd (sipa / rē'û), who was responsible for taking care and ensuring the day-to-day management of the herd, most probably during a year's tenure which ended at the plucking time in the spring (Snell 1986: 138).

The shepherd duties and their wages in kind, consisting for example of wool and his own share of the growth of the flocks, were sometimes recorded in contracts documented at different periods. The exact composition entrusted to a shepherd is usually indicated (Meek AJSL 33, 3) and the expected yield of the growth of the herd that the shepherd has to reach as well as his wage paid in wool. Losses, with a bearable rate accepted of 10%, are placed under the responsibility of the shepherd who has to replace the lost animals.

'1 ewe, brand (šimtum) of Amurrum-tukulti, 1 ewe, brand of Abba, 1 ewe, brand of Tawirum, 1 ram. Total: 4 skins of dead (litt. fallen) (animals)' (UET 5, 611: 1-7), (Béranger 2020: 254).

Branding the animal is used to identify who was in charge of it or to whom it belongs, as is also indicated in the (Hammurabi Codex, Law 265).

Third millennium and early second millennium BC cuneiform sources provide an impressive amount of data on sheep. This presentation focused on breeds, physical characteristics, provenance as well as sheep husbandry including for example feeding, reproduction and hybridization. In this ongoing study, other aspects were left aside such as the by-product that characterize the animal as a resource fully utilised in all its various forms, i. e. wool, skin, sinews, carcase, dairy and meat products. As for the latter, it is interesting that the animal parts such as the head, neck, chest, shoulder, scapula, thigh and leg are also documented as pieces of meat for consumption not to mention the terminology of the innards in relation to oracular consultation which is primarily based on the sacrifice of lambs.

As philologists, it is extremely important for us to be part of the multidisciplinary Evosheep project which enables the necessary discussion between specialists of different fields, thus improving our understanding of the written documentation on such issues as the process of hybridization or the reasons behind the different proportions of animals in herds.

Fig. 1

 1 Fig. 1 Old Babylonian Nippur ura3 lexical list CDLI P228711©Cuneiform Digital Library Initiative

Fig. 2

 2 Fig. 2 Fat tailed Awassi breed Lebanon ©EvoSheep

Fig. 3

 3 Fig. 3 Long tailed Gumz breed Ethiopia ©EvoSheep A. Amane

  Fig. 4 Brown-red sheep-Menz Breed Ethiopia ©EvoSheep

Fig. 5

 5 Fig. 5 Spotted sheep Jacob Sheep ©Livestock of the World

Fig. 6 :

 6 Fig. 6 : Range of the grazing flocks according to FM 6, 9. Shaded in orange: extension of the Mari kingdom under Zimri-Lim (1774 B.C -1762 B.C), M. Sauvage (dir.), Atlas historique du Proche-Orient ancien, Les Belles Lettres, Institut français du Proche-Orient 2020, p. 82

Acknowledgements: EvoSheep was funded by the French National Research Agency (ANR).


domain: stat sde

text:
Camille Champion

Raphaelle Momal

Emmanuelle Le Chatelier

Mahendra Mariadassou
email: mahendra.mariadassou@inrae.fr

Magali Berland
email: magali.berland@inrae.fr

OneNet -One network to rule them all: consensus network inference from microbiome data

Keywords: Network inference, stability selection, microbial ecology, microbial guild, Gaussian Graphical Models

Introduction

The human gut microbiota is a complex ecosystem consisting of trillions of microorganisms, mainly viruses, bacteria, archeae and microbial eucaryotes, that play critical roles in host physiology including digestion, immune function and metabolism [START_REF] Belkaid | Role of the microbiota in immunity and inflammation[END_REF]Hand, 2014, Chatelier et al., 2013]. Recent advances in sequencing technologies have enabled the characterization of gut microbiota composition and function at a fine scale, providing opportunities to understand the microbial communities that reside within the human gastrointestinal tract. However, despite these technological advancements, understanding the interactions within the bacteria of the gut microbiota remains a major challenge. These interactions are complex as microorganisms can interact with each other in a multitude of ways: through mutualism, parasitism, commensalism and competition to only cite a few [Weiss et al., 2017, Faust and[START_REF] Faust | Conet app: inference of biological association networks using cytoscape[END_REF]. To address this challenge, network-based approaches have been developed to infer microbial interactions and construct microbial interaction networks. The resulting networks can reveal potential interactions between microbial taxa and support the identification of microbial guilds. Those guilds are defined as groups of microorganisms that co-occur and may interact with each other. Identifying microbial guilds is crucial for understanding the ecological dynamics of the gut microbiota and can provide insights into the role of the microbiota in health and disease [START_REF] Wu | Guild-based analysis for understanding gut microbiome in human health and diseases[END_REF][START_REF] Xiao | Large-scale microbiome data integration enables robust biomarker identification[END_REF].

Formally, microbial interaction networks consist of nodes, which correspond to microbial species, and edges, which correspond to interactions between those species. Positive and negative interactions are rarely observed directly. They are instead often reconstructed from abundance data, using either longitudinal data (see the generalized Lotka-Volterra model in [START_REF] Bucci | Mdsine: Microbial dynamical systems inference engine for microbiome time-series analyses[END_REF]) or co-occurrence data. We focus here on the latter suite of methods.

The simplest way to identify microbial interactions is to perform a correlation analysis. However, correlation-based methods model total dependencies and are therefore prone to confusion by environmental factors (e.g. shared habitat preferences or susceptibility to the same abiotic factors) and do not lend themselves to a clear separation between indirect and direct effects [START_REF] Friedman | Inferring correlation networks from genomic survey data[END_REF]. By contrast, conditional dependency-based method eliminate indirect correlations from direct interactions and lead to sparser and easier to interpret networks, at the cost of increased computational burden and more sophisticated models. The problem of network inference is complicated by the adverse characteristics of microbial abundance data, which are sparse, heterogeneous, heteroscedastic and show extreme variability. These data are thus tricky to model, which leads to poorly reproducible and/or sparse networks, with many missing edges [START_REF] Peschel | Netcomi: network construction and comparison for microbiome data in r[END_REF].

The most common framework for the estimation of the conditional dependencies is Gaussian Graphical Models (GGM) [START_REF] Lauritzen | Graphical Models[END_REF], which describe the conditional dependency structure of multivariate Gaussian distributions. As microbiome abundance data don't directly fit within the gaussian framework, three main workarounds are commonly used: data transformation, models based on alternative distributions and models based on latent variables ; the whole strategy is also illustrated on Fig 1, and each step is detailed in the Supplementary Methods (Section 6).

The complexity of reconstructing networks from co-occurrence data has spawned a rich literature with many methods relying on the solutions exposed above, including (i) approaches based on correlation as SparCC [START_REF] Friedman | Inferring correlation networks from genomic survey data[END_REF], CoNet [START_REF] Faust | Conet app: inference of biological association networks using cytoscape[END_REF], (ii) approaches based on probabilistic graphical models as SpiecEasi [START_REF] Kurtz | Sparse and compositionally robust inference of microbial ecological networks[END_REF], gCoda [START_REF] Fang | gcoda: Conditional dependence network inference for compositional data[END_REF], SPRING [START_REF] Yoon | Sparse semiparametric canonical correlation analysis for data of mixed types[END_REF], PLNetwork [START_REF] Chiquet | Variational inference for sparse network reconstruction from count data[END_REF], ZiLN [START_REF] Prost | A zero inflated log-normal model for inference of sparse microbial association networks[END_REF], cozine [START_REF] Ha | Compositional zero-inflated network estimation for microbiome data[END_REF], Magma [START_REF] Cougoul | Magma: Inference of sparse microbial association networks[END_REF], EMtree [START_REF] Momal | Tree-based inference of species interaction networks from abundance data[END_REF] and (iii) approaches based on the inference of the latent correlation structure as CCLasso [START_REF] Fang | Cclasso: correlation inference for compositional data through lasso[END_REF], SparCC [Friedman and Alm,2 Figure 2: High level summary of the OneNet pipeline: (i) bootstrap subsamples are constructed from the original abundance matrice, (ii) each inference method is applied on the bootstrap subsamples to compute edge selection frequencies using a fixed λ grid, (iii) a different λ is selected for each method to achieve the same density in all methods, (iv) edge selection frequencies are summarized and (v) thresholded to compute the consensus graph. 2012]. The most recent methods include: mixPLN and ms-mixPLN [Tavakoli and[START_REF] Tavakoli | 5: Compared recall (TPR) of inference methods and OneNet-* variants after removing glassobased methods from the set of methods[END_REF][START_REF] Yooseph | Variational approximation-based model selection for microbial network inference[END_REF] which consider the problem of inferring multiple microbial networks (one per host condition) from a given sample-taxa abundance matrix when microbial associations are impacted by host factors. [START_REF] Jiang | Harmonies: A hybrid approach for microbiome networks inference via exploiting sparsity[END_REF] proposed Hybrid Approach foR MicrobiOme Network Inferences via Exploiting Sparsity (HARMONIES) which addresses some critical aspects of abundance data (compositionality due to fixed sampling depth, over-dispersion and zero-inflation of the abundances) while maintaining computational scalability and sparsity of the interaction network, in contrast to mixPLN and ms-mixPLN. Finally, Network Construction and comparison for Microbiome data (NetCoMi) [START_REF] Peschel | Netcomi: network construction and comparison for microbiome data in r[END_REF], provides a one-step platform for inference and comparison of microbial networks, by implementing many existing methods for abundance data preprocessing, network inference and edge selection in single package.

All these methods have been designed to infer networks based on different mathematical hypotheses and thus have different strengths and weaknesses when modeling microbiome data. Each microbial network inference algorithm usually returns distinct edges to connect the taxa together, as many facets of the same reality. In this article, we present OneNet, an ensemble method that generate robust and reliable consensus network that will facilitate the identification of microbial guilds and generation of new hypotheses.

Methods

Overview of OneNet

We developed OneNet, a three-step procedure for robust consensus network reconstruction, illustrated on Figure 2.

We included seven inference methods in OneNet, all of which rely on Gaussian Graphical Models (GGM) to estimate conditional dependencies networks: Magma, SpiecEasi, gCoda, PLNetwork, EMtree, SPRING and ZiLN. We excluded mixPLN and ms-mixPLN as they do not reconstruct a single network but rather a collection of networks and NetCoMi as it collects already existing methods rather than introducing a new one. We left out HARMONIES from the comparison as its implementation doesn't allow the user to specify the regularization grid, a crucial step in the ensemble method, and achieved worse performance than included methods in preliminary tests. We also excluded cozine from OneNet because its implementation doesn't rely on resampling and prevents it from being integrated, but we nonetheless 1: Characteristics of the network inference methods: general methods, abundances normalization (Centred Log Ratio (CLR), Geometric Mean of Pairwise Ratios (GMPR), Relative Log Expression (RLE)), distribution transformations, inference approaches (Meinshausen-Bühlmann (MB), glasso, tree averaging), covariates integration and references. OneNet is based on all methods but cozine.

included it in the benchmark as it compared favorably to others methods in preliminary tests. Table 1 summarizes the inference strategies adopted by each method and the potential integration of covariables in the model.

2.2

Step 1: Assign to each edge a sequence of selection frequency values for each inference method.

Each inference method assigns a score to the edges: either a probability (for the tree averaging method) or the maximal penalty level λ below which the edge is selected in the network. An optimal penalty λ * on these scores is then needed for an edge to be selected in the final network. Several approaches exist but the concept of stability selection [START_REF] Liu | Stability approach to regularization selection (stars) for high dimensional graphical models[END_REF] is the most widely used and the one considered in this work as it yields a compromise between precision and recall, while fostering reproducibility. The associated method, called Stability Approach to Regularization Selection (StARS) uses a resampling strategy to select the value of λ * leading to the most stable graph. We describe briefly the StARS algorithm formally in the following and the modification we propose in this work.

StARS algorithm.

The original data X is subsampled B times and the network inference is conducted on each sub-sample for each value of λ in a grid (λ 1 , . . . , λ K ) to obtain a graph G b,k , with k ∈ {1, ...K} and b ∈ {1, . . . , B}. The selection frequency of edge e for parameter λ k , is computed as its selection frequency across the subgraphs:

f k e = B b=1 1 {e∈G b,k } .
The selection frequency over resamples gives an idea of edge reproducibility: frequency and robustness of the edges are clearly related. StARS aggregate those frequencies to construct a network-level measure of edge variability defined as:

S k = 1 -4 1 q e f k e (1 -f k e )
where q = p(p -1)/2 is the total number of possible edges and S k can be thought of as the mean of (edge-level) Bernoulli variances. Each value λ k is associated to a single selection frequency vector, and a resulting stability value. Finding the right edge frequency is therefore equivalent to finding the right stability level. Classical choices for stability are stab = 80% or stab = 90% to have a good compromise between recall and precision and the optimal level λ * is chosen as λ * = min λ S(λ) ≥ stab. Once the optimal level λ * is fixed, one solution is to return the refit coefficient (1 if the edge is selected and 0 otherwise) computed by running inference methods on the complete dataset, to get the final network.

2.2.2

Using edge-level selection frequencies rather than network-level stability.

Instead of computing the network-level stability S k and doing a final refit step, we propose instead to select directly the edges with high selection frequency for penalty λ as E λ (c) to create the set E λ (c) = {e, f λ e > c}, where c is a constant close to 1. In this way, we guarantee both high precision and high reproducibility for edges in E λ (c) as they are selected many times in the resampling. Smaller values of λ give rise to larger sets E λ (c) and higher recalls. Two advantages of using frequencies rather than refitting the network are (i) filtering out edges with low support that could be included in the refit graph and (ii) making it easier to combine the edges inferred by the different methods.

2.3

Step 2: Equalize of the densities of the networks.

In order to include the best edges in the consensus network, we must choose one λ per method. A natural choice would be the value λ ⋆ m selected by StARS for the method m. However, we observed that StARS computes a very different precision/recall for each method. We select instead the smallest λ m such that (i) the sets E λm m (c) are roughly of equal sizes and (ii) the mean stability is above a given threshold:

1 M M m=1 S m (λ m ) ≥ stab.
It forces all methods to contribute with a similar number of edges to the consensus while ensuring that each edge set is reproducible. In practice, to match edge set sizes, we worked with edge density rather than with λ values as the two are monotonically related.

2.4

Step 3: Summarize the sets of selection frequencies across methods.

Build a consensus network from the sets of edges E λ (c) produced by the different methods is akin to an ensemble procedure where many methods are combined together.

In order to produce a stable and accurate consensus network, we define several summary metrics in the objective to mitigate the drawback of each method while benefiting from their strength. The consensus is obtained by summarizing edges frequencies across the methods. Denoting by f m the selection frequency of a given edge with method m ∈ {1, ..., M }, we define:

• mean: average selection frequency m f m M , • norm2: euclidean norm (2-norm) ( m f 2 m ) 1/2 M 1/2 , • IVW: inverse-variance weighted average m f m × 1 fm(1-fm) m 1 fm(1-fm) , where f m follows a Bernoulli with V ar(f m ) = f m (1 -f m ),
• min3: high frequency for at least 3 methods 1 { m (fm>c)≥3} .

Note that min3 is the only one that returns a binary summary, all the other ones take value in [0, 1], like the original selection frequencies.

Datasets

Simulated dataset

In this work we simulate data using the methodology described in [START_REF] Yoon | Sparse semiparametric canonical correlation analysis for data of mixed types[END_REF] which is based on gaussian copula to control the network structure followed by sampling from the species marginal distributions to preserve the peculiarities of abundance data. This method yields synthetic data with marginal distributions that are closer to the original empirical dataset, while enforcing a given correlation structure between the species.

To generate the simulated dataset, we use in input the empirical dataset described below restricted to diseased individuals. The dataset is simulated in the framework of an unknown undirected graph G(V, E), with no retroactive loop, consisting of p vertices V = {1, . . . , p} and a set of edges E ⊆ V × V connecting each pair of vertices. The graph G is represented by its adjacency matrix A = (A ij ) (i,j)∈E of size p × p, defined as:

∀(i, j) ∈ 1, p 2 , A ij = 1 if (i, j) ∈ E, 0 otherwise.
The package 'EMtree' v.1.1.0 [START_REF]Emtree: Infers direct species association networks using tree averaging[END_REF] is used to generate a precision matrix Ω defined as the graphical Laplacian A of a cluster graph. Ω is inverted to create the correlation matrix Σ and the idea was then to simulate variables with arbitrary marginal distributions from multivariate normal variables with correlation structure given by Σ using gaussian copula. Specifically, we generate a n × p matrix Z with independent normal rows Z i ∼ N (0, Σ). We then get uniform random vectors by applying standard normal cdf transformation to each column of Z, u j = ϕ(Z j / Σ jj ) element-wise, and finally apply the quantile functions of the empirical data marginal distributions to each u j . The function synthData_from_ecdf from the 'SPRING' package v.1.0.4 [START_REF] Yoon | Semi-parametric rank-based approach for inference in graphical model (spring)[END_REF] is used for these simulation steps. To assess the effect of sample size, we simulate datasets of size n ∈ {50, 100, 500, 1000}.

Evaluation criteria

Each methods is evaluated by comparing the inferred network structure to the known simulated network structure using the following metrics:

• Precision (positive predictive value): PPV= TP/(TP+FP),

• Recall (true positive rate): TPR= TP/(TP+FN),

where TP stands for True Positive (a correctly detected edge), FP for False Positive (an edge detected where none should be) and FN for False negative (an undetected edge). The precision measures the proportion of real edges among what the detected ones, whereas the recall measures the proportion of real edges which are detected.

Empirical dataset

The empirical dataset, studied in [START_REF] Qin | Alterations of the human gut microbiome in liver cirrhosis[END_REF], corresponds to stool samples from 216 Chinese individuals sequenced using whole-metagenome sequencing techniques. The raw sequences are available as BioProject PRJEB6337 in the European Nucleotide Archive (ENA). Among this population, 102 individuals are healthy and 114 suffer from liver cirrhosis. Abundances of all microbial species (metagenomic species or MSP) detected using 10.4 million IGC2 gut gene catalogue [START_REF] Wen | Quantitative metagenomics reveals unique gut microbiome biomarkers in ankylosing spondylitis[END_REF] are extracted using the Meteor software suite that create a gene abundance table by mapping high quality reads onto the gene catalogue, using Bowtie2. Abundance of each MSP is computed as the mean abundance of 100 marker genes selected for that MSP, where the gene abundance is the read abundance normalized by the gene length. The final table of size 1990 MSP by 216 individuals records all the normalized abundances [START_REF] Champion | Microbial species abundances from public project PRJEB6337 on liver cirrhosis[END_REF]. In the Application section, we used only the 114 cirrhotic patients.

Results

In this section we evaluated the performances of both OneNet and the network inference methods on the simulated dataset.

Influence of the stability level on the inferred graphs

We first evaluated the effect of stability level on the performance of the inference methods. Instead of fixing a target stability at 0.8 or 0.9, we studied the relationship between the precision and recall of the inferred edges by each method for different stability levels. Because interactions between highly prevalent species are easier to reconstruct, we only kept metagenomic species with a prevalence greater than 50% (159 species). We let the sample size vary in {50, 100, 500, 1000} and we considered B = 40 resamples each time.

Methods have distinct precisions for a given stability level. Figure S1 shows the relationship between the performance obtained with the edge set E λ * (0.90) (precision PPV90 and recall TPR90), and the corresponding stability. The difference in patterns grows with sample size n, revealing peculiarities inherent to each method. Clearly, methods have distinct performances for the same stability level. We observed that methods cluster in groups (glasso, neighorhood selection, tree aggregation) with different precision/recall tradeoffs. As a result, they produce edge sets that greatly differ both in size and quality. This suggests that the stability value is not a good indicator of the precision level achieved by each method.

Methods have comparable precision and recall for a given density level. Unlike precision, which is unavailable when dealing with empirical datasets, the density, or number of detected edges, can always be computed. Figure S2 shows the relationship between precision (resp. recall) and density for all methods at increasing sample sizes. The curves are almost superimposed for values of n up to 100, after which different behaviors appeared. However, the gap in performance between methods stayed small when imposing the same density, rather than the same stability. This also meant that, whatever the method used, the m first edges included in the graph achieve similar graph reconstruction quality, although they correspond to different stabilities. We thus selected individual graphs based on density, rather than stability, to include only graphs with similar precision and recall in the consensus phase.

Mean stability as a proxy of the density level. The previous observation prompted us to explore the link between density and stability for different values of λ. Figure S3 shows how stability decreases with increasing density. We set a target mean stability value (e.g. 90%) for each value of n (here between 50 and 1000). As n grows, we observed that the density increased from 90 to 170, as well as the spread between methods. For n = 1000, stabilities ranged from 0.8 The TPR/PPV compromise achieved for λ * corresponding to a stability of 0.9 is shown with a circle, the one achieved by a mean stability across methods of 0.9 is shown with a square. Whenever the selected λ is the same, the circle and the square are replaced with a diamond. Finally, note that cozine relies on minimization of a BIC criteria (shown with a triangle) rather than on the resampling-based stability selection to choose the regularization parameter.

for EMtree to 1 for SPRING. We can see how targeting the mean stability rather than the same stability for all methods allows to adapt the precision level of each method through density to make them more similar.

Mean stability increases the consistency between the performances of the network inference methods. We compared, for different sample sizes, the precision -recall tradeoff achieved by the mean stability to the ones achieved by a fixed stability (e.g. 0.9). Figure 3 shows that the sample size has a major impact on the precision and stability. The ROC curves stabilize to near-perfection starting at n = 500 (Figure 3 c). It is also noticeable that the adapted stabilities reduced the range of the method's precision. Furthermore, for glasso-based methods (gcoda, PLNnetwork), the new target led to a 20 points improvement in TPR at almost no cost in PPV, for large sample sizes.

OneNet versus the classical network inference methods

We computed, for a frequency threshold of 90%, the precision and recall values obtained by the classical network inference methods (cozine, gCoda, PLNnetwork, SPRING, Magma, SpeicEasi, ZiLM and EMtree) and we compared them to the OneNet networks (with the summay metrics mean, norm2, IVW ans min3). Note that because of the similar density, each method provided roughly the same number of edges to OneNet. Figure 4 shows that OneNet with the mean and norm2 consensus methods, achieved the best precision levels but the worse recall values. OneNet with the min3 summary has comparatively lower precision but higher precision and OneNet with the IVW summary has both worse recall and worse precision. This illustrates how OneNet generally leads to sparser networks with higher precision.

We also observed that the MB-based ones tend to outperform glasso-based ones. However, OneNet still demonstrated precision levels equivalent to that of the best methods (Figures S4 andS5). This reflects the inherent robustness of consensus measures to methods with outlying performance.

On top of that, sample size n has a dramatic effect on all criteria and methods. For large sample sizes (n ≥ 500), most methods exhibit precision above 99% and recall between 30% and 60% (Fig. S6 and S7, panels (c) and (d)). By contrast (Fig. S6a drops below 60% (resp. 80%) for all methods except OneNet -Mean and OneNet -norm2 which both remain above 75% (resp. 85%). Likewise, the recall drops below 40% for n = 100 (Fig. S7b) and below 20% for n = 50 (Fig. S7a).

We observed a discrepancy in terms of precision and recall for the cozine method between figures S6 and S7. We hypothesized that it's due to the original cozine procedure (BIC criteria) used to select the optimal network, which selects to a dense graph (very high recall, very low precision, and therefore many spurious edges, see Figure 3). By contrast, the Figures S6 andS7 show the precision and recall values obtained with the resampling approach. This is an extreme example of lack of robustness, where the network reconstructed from the full dataset differs drastically from the ones reconstructed on random subsets of the data and illustrates the benefits of combining the resamples rather than doing a refit.

Application to liver cirrhosis

To investigate the potentiality of OneNet relative to the other methods, we inferred all the networks from the microbiome dataset of cirrhotic patients presented in the Methods [START_REF] Qin | Alterations of the human gut microbiome in liver cirrhosis[END_REF]). Because of the small size of the dataset (114 samples), only metagenomic species with a prevalence greater than 50% were kept (155 species). The networks have been clustered using the CORE-clustering algorithm to reconstruct microbial guilds [START_REF] Champion | Detection of representative variables in complex systems with interpretable rules using core-clusters[END_REF]).

Following the guidelines of this paper, we fixed the number of clusters between 10 and 19.

Figure 5 illustrates the inferred and clustered networks. We note that the OneNet network, with the mean summary, is the sparsest one and is closely related to the SpiecEasi and the SPRING's one. This result is in line with the results on the simulated data and with Figures 4 andS8.

The second objective of this application was to characterize microbial species responsible for liver cirrhosis using OneNet. One guild -the "cirrhotic guild", represented in Figure 6, is notable as it contains species of the genera Anaeroglobus, Campylobacter, Haemophilus, Prevotella, Streptococcus known in the literature to be associated to the following diseases: liver cirrhosis [START_REF] Qin | Alterations of the human gut microbiome in liver cirrhosis[END_REF], obesity after weight-loss intervention [START_REF] Liu | Gut microbiome and serum metabolome alterations in obesity and after weight-loss intervention[END_REF], schizophrenia [START_REF] Zhu | Metagenome-wide association of gut microbiome features for schizophrenia[END_REF], atherosclerotic cardiovascular disease [START_REF] Jie | The gut microbiome in atherosclerotic cardiovascular disease[END_REF] and Crohn disease [START_REF] He | Two distinct metacommunities characterize the gut microbiota in crohn's disease patients[END_REF].

8

Figure 5: Networks inferred on the liver cirrhosis data set with EMtree, gCoda, PLNnetwork, SPRING, ZiLN, Magma, SpiecEasi and OneNet -mean followed by Core-clustering algorithm to identify the microbial guilds. All graphs are shown using the same node layout for ease of comparison and the nodes are coloured by cluster. Methods are grouped based on the underlying inference technique (tree aggregation, glasso, neighborhood selection).

Figure 6: Detailed view of the cirrhotic guild identified in the OneNet -mean network with taxonomic information on the nodes. All species known to be associated with diseases are annotated with a star.

Discussion

The proposed framework, with a microbial consensus network inference method, offers new insights about inferring robust and sparse microbial networks. OneNet is robust in the sense that i) it uses GGM adapted to deal with the peculiarities of microbial abundance properties (inclusion of environmental effects as covariates, stabilization of data variability, adaptation to abundances with high proportion of zeros, etc), ii) it depends on seven network inference methods aiming for sparse and reproducible microbial network using either glasso, neighborhood selection or tree averaging approaches, iii) it relies on a three-steps procedure to improve the precision and reproducibility of both inference methods and OneNet. Indeed, the selection of edges with high inclusion frequencies and harmonization of stability selection achieve similar precision levels accross methods. The resulting consensus network uses a summary of the edge inclusion frequencies.

Results from the studies on synthetic and real data illustrated the first major and reassuring fact, that the inference methods overall agree with one another and with the truth (Figures S9 andS10). It then showed the effectiveness of OneNet compared to the inference methods. Among the different summaries considered, the mean or norm2 are preferred since they lead to slightly sparser networks but achieved much higher precision than any inference method, especially for sample sizes around n = 100, which is typical in microbiome studies. By contrast, min3 and IVW summaries gave a significant additional quantity of edges compared to the other summary metrics, yielding TPR levels that are comparable to those obtained with glasso-based methods without increasing the PPV, especially when the number of samples is small (n ≤ 100) (Figures S9).

In all numerical experiments, we showed that a minimal sample size to maintain high robustness was n = 100. In this scenario we suggest to use the mean summary in OneNet, as it proved to be more robust to small sample sizes.

Obviously, the precision is affected by both the number of samples and microbial species in the system, the latter being controlled by the prevalence threshold imposed at the very beginning of the analysis. As illustrated on Figures S11, S12 and S13, the prevalence threshold can be adjusted to increase the precision of the method depending on the number of samples (prev = 0.50 for n = 100 and n = 500, prev = 0.20 for n = 1000). From n = 1000, when considered individually, the neighborhood selection and the tree averaging approaches showed performances that were similar to OneNet. In this context, it could be possible to select one of these three approaches.

An advantage of OneNet is it's ability to easily incorporate new inference methods as soon as they are amenable to the modified stability selection framework used here. This is the case for all the methods considered in this work but cozine. Indeed, cozine relies on the BIC criteria to tune the regularization parameter λ: it doesn't allow the user to provide a fixed λ grid for comparisons with other methods and doesn't produce the table of edge selection frequencies required to compute summaries. This is however due to implementation choices (using BIC instead of StARS for selecting λ) rather than to fundamental incompatibilities with the OneNet framework.

Supplementary Methods

We detail here each step of the whole network inference strategy, illustrated on Fig 1 . As microbiome abundance data don't directly fit within the gaussian framework, three main workarounds are commonly used: data transformation, models based on alternative distributions and models based on latent variables:

Transformations. A small constant is added to each abundance before log-transforming them. However, this transformation does not stabilize data variability because the log-transformed abundances scale with sequencing depths and covary with it, making dependencies modeling tricky. On the contrary, the Centered Log Ratio (CLR) transformation [START_REF] Aitchison | The statistical analysis of compositional data[END_REF] guarantees the study of dependencies. It is however highly criticized when data contain a high proportion of zeros: proportions higher than 90% are typical in whole-metagenome or amplicon sequencing data. To circumvent this problem, [START_REF] Yoon | Sparse semiparametric canonical correlation analysis for data of mixed types[END_REF] introduced a modified version of the CLR transformation (mCLR) that respects the original ordering of the data but doesn't account anymore for the compositional nature of the data, the primary motivation of the CLR transformation. An alternative to compositional transformations is to use a normalization factor, such as Geometric Mean of Pairwise Ratios (GMPR) [START_REF] Chen | Gmpr: A robust normalization method for zero-inflated count data with application to microbiome sequencing data[END_REF], Relative Log Expression (RLE) [START_REF] Anders | Differential expression analysis for sequence count data[END_REF] and others like Wrench normalization factors [START_REF] Kumar | Analysis and correction of compositional bias in sparse sequencing count data[END_REF] or Cumulative Sum Scaling [START_REF] Paulson | Differential abundance analysis for microbial marker-gene surveys[END_REF]. The GMPR normalization is designed for abundances with a high proportion of zero values. It compares pairs of samples based only on the species they share, and considers the geometric mean of the median ratio. This makes this technique robust to both differentially abundant species and extreme values. When two samples do not share any species, the computation of GMPR fails (this happens when samples come from very contrasted conditions with no or limited species overlap). In this case, Relative Log Expression (RLE) normalization method can be used. This method is based on the assumption that most of the species are not differentially abundant.

However, this normalization factor fails when no single species is shared across all samples, which is frequently the case in microbiome data. A modified version of RLE only considers positive abundances to avoid this drawback.

Distributions and models. The second workaround is to use models adapted to abundance data characteristics : overdispersion (excess of variability in the data) and zero-inflation (excess of zeros). The Poisson-log normal (PLN) model [START_REF] Chiquet | The poisson-lognormal model as a versatile framework for the joint analysis of species abundances[END_REF] is designed for the analysis of abundance tables. It accounts for both structuring factors and potential interactions between the species. In the presence of overdispersion, the Poisson regression model is not adequate and can lead to biased parameter estimates and unreliable standard errors estimates. The Negative Binomial (NB) model is then often used [START_REF] Forbes | 9: Upset plot of the edges identified by the inference methods, the OneNet-* variants and the ground truth on the synthetic dataset for n = 100[END_REF]. Both models can be seen as compound Poisson model (with a lognormal for the PLN distribution and Gamma for the NB) that are overdispersed compared to a vanilla Poisson distribution but the PLN is multivariate and can account for correlations between abundances. Contrary to the NB model, the zero-inflated model [START_REF] Greene | Accounting for excess zeros and sample selection in poisson and negative binomial regression models[END_REF] is often motivated by an excess of zeros in the data, but less flexible than the zero outcome model. An intuitive approach to analyzing zero-inflated abundance data is to view the data as arising from a mixture distribution of a point mass distribution at zero and an abundance distribution. Hurdle models [START_REF] Cragg | Some statistical models for limited dependent variables with application to the demand for durable goods[END_REF] are a class of models for abundance data that help handle excess zeros and overdispersion. In contrast to Zero inflated-models, hurdle models capture both an excess or a lack of zeros in the dataset. The zero-inflated negative binomial (ZINB) model [START_REF] Cheung | Zero-inflated models for regression analysis of count data: a study of growth and development[END_REF], obtained by applying ZI to NB model, takes into account both overdispersion and excess of zeros. Finaly, copulas are a multivariate cumulative distribution functions for which the marginal probability distribution of each variable is uniform on the interval [0, 1]. As they fully describe the dependency structure, models with copulas allow to separate the modeling of marginal distributions (e.g. overdispersed, with excess zeros, etc) from the modeling of dependencies. Recent developments used gaussian copula coupled with arbitrary discrete marginal distributions to study multivariate abundance data [START_REF] Anderson | A pathway for multivariate analysis of ecological communities using copulas[END_REF]. [START_REF] Popovic | Untangling direct species associations from indirect mediator species effects with graphical models[END_REF] showed that Gaussian copulas are a relevant and promising approach to the problem of network inference from abundance data, even if the computational cost is higher than for other methods. One way of taking advantage of the copula theory without having to actually estimate the joint distribution is to use copulas as a sophisticated data transformation technique to transform abundances into pseudo-Gaussian data.

Latent variables. The third popular idea is to model multivariate discrete data using latent variables and push the depdendency back to the latent layer. Latent variables models have recently received increasing attention as they provide a convenient way to model the dependence structure between species. Two specifications of latent variables stand out in community ecology [START_REF] Warton | So many variables: joint modeling in community ecology[END_REF]: the Multivariate Generalized Linear Mixed Model (GLMM) [START_REF] Ovaskainen | Modeling species co-occurrence by multivariate logistic regression generates new hypotheses on fungal interactions[END_REF][START_REF] Pollockand | Understanding co-occurrence by modelling species simultaneously with a joint species distribution model (jsdm)[END_REF], and the Latent Variable Model (LVM) [START_REF] Ovaskainen | Using latent variable models to identify large networks of species-to-species associations at different spatial scales[END_REF][START_REF] Ovaskainen | How to make more out of community data? a conceptual framework and its implementation as models and software[END_REF]. The difference between these models lies in the dimension of their respective random effects: there are as many latent variables as there are species in the GLMM, whereas in the LVM their number is a parameter of the model.

Most methodologies to infer networks from abundance data first use a rationale (data transformation, latent variable modeling, etc.) to solve the problem of network inference in the Gaussian setting. There, they take advantage of the GGM framework to perform network inference using penalized likelihood or tree-based approaches to estimate the precision matrix, from which is finally derived the network.

Penalized likelihood approaches. There exist two main penalized approaches for the estimation of GGM: the graphical LASSO (glasso) [Friedman et al., 2008], and the neighborhood selection, also called the Meinshausen-Bühlmann approach (MB) [Meinshausen and Bühlmann, 2006]. Both are penalized likelihood approaches which perform a sparse estimation of the precision matrix, either all at once for the glasso or row by row for MB.

Tree averaging approach. Another GGM inference method considered in this article is the tree averaging approach [Meila and Jordan, 2000], which leverages specific algebraic properties to perform a complete and efficient exploration of the space of spanning tree structures. Note that this approach does not require the GGM Markov faithful property to hold. Each edge is given a posterior probability of being present in the network and those probabilities are thresholded to build the network.

Figure 1 :

 1 Figure 1: The classical network inference pipeline.

Figure 3 :

 3 Figure 3: Precision -recall curves of each inference method for different sample sizes: (a) 50 (b) 100 (c) 500 (d) 1000.The TPR/PPV compromise achieved for λ * corresponding to a stability of 0.9 is shown with a circle, the one achieved by a mean stability across methods of 0.9 is shown with a square. Whenever the selected λ is the same, the circle and the square are replaced with a diamond. Finally, note that cozine relies on minimization of a BIC criteria (shown with a triangle) rather than on the resampling-based stability selection to choose the regularization parameter.

Figure 4 :

 4 Figure 4: Quality of both individual and aggregated frequencies in terms of PPV (top row, values above 0.6) and TPR (bottom row, values below 0.5) for n = 100


David Roger
email: droger@inrialpes.fr

Ulf Assarsson

Nicolas Holzschuch
email: holzschu@cs.cornell.edu

Efficient Stream Reduction on the GPU

Stream reduction is the process of removing unwanted elements from a stream of outputs. It is a key component of many GPGPU algorithms, especially in multi-pass algorithms: the stream reduction is used to remove unwanted elements from the output of a previous pass before sending it as input for the next pass.

In this paper, we present a new efficient algorithm for stream reduction on the GPU. Our algorithm works by splitting the input stream into smaller components of a fixed size, on which we run a standard stream reduction pass. We then concatenate the results of these stream reduction pass with line drawing. This last pass is very efficient because we already know the size of the data.

I. INTRODUCTION

Stream reduction [START_REF] Horn | Stream reduction operations for GPGPU applications[END_REF], sometimes called stream compaction [START_REF] Sengupta | A work-efficient stepefficient prefix sum algorithm[END_REF] is the process of removing unwanted elements from a data stream. It is especially useful for multi-pass GPGU algorithms, where the stream of data output from a pass is used as input in the next pass: in the first pass, the algorithm finds that some data is not required, and flags them for deletion. Since fragment shaders lack the ability to write data at specific places in memory (scatter), these elements are still present in the output stream, and the stream reduction pass must remove them before the next pass.

The stream reduction pass, while being essential in many GPGPU algorithms, is a difficult problem, because it implies writing data at specific locations. Given the specificity of the GPU architecture, it is often more efficient to do a stream reduction pass on the entire data set after a processing pass than to try to do scattered writing on the fly during the processing pass. The consequence is that the stream reduction pass must work efficiently on very large datasets.

Because of the size of the datasets and the difficulty to specify the address for writing memory access, stream reduction passes are the bottleneck in many GPGPU algorithms [START_REF] Horn | Stream reduction operations for GPGPU applications[END_REF], [START_REF] Roger | Whitted ray-tracing for dynamic scenes using a ray-space hierarchy on the GPU[END_REF].

In this paper, we present a new approach to stream reduction. Our algorithm is hierarchical, and splits the input stream into smaller components of a fixed size, s. A stream reduction pass is run on each of the components, and the results of each pass are concatenated using line drawing. The advantages of our hierarchical approach are numerous: stream reduction algorithms are more efficient on smaller data sets, and the final pass of concatenation has a very low cost, because we already know the size of the data sets.

In our study, we have found that the theoretical complexity of our hierarchical stream reduction algorithm is O(n + n log s), where s is a constant, compared to O(n log n) for previous work. Experimental studies confirm that our stream reduction algorithm is substantially faster than existing algorithms.

1: i ← 0 2: for j = 0 to n -1 do 3: if x[j] must be kept then 4: x[i] ← x[j] 5: i ← i + 1
Recent GPUs (such as NVidia G80) offer the ability to easily remove elements from a stream of input data, through the geometry engine. Our experiments show that our stream reduction algorithm is also faster than this approach, especially for large datasets.

However, our algorithm is mostly destinated toward frameworks that lack scattering abilities, such as OpenGL. Other systems that are able to perform scatter operations, such as CUDA, can perform reduction through a more straightforward process, and would benefit less from our algorithm.

Our paper is organized as follows: in the next section we review previous work on stream reduction on the GPU. In section III, we then describe our algorithm in detail, and conduct a study of the theoretical complexity of our algorithm (section III-D). Section IV is a review of the different available prefix sum scan techniques, that we need as a subroutine. In section V, we study the experimental behavior of our algorithm as a function of input size and parameters, and compare it with state-of-the-art algorithms. In sectionVI, we present our conclusions and directions for future works.

II. PREVIOUS WORKS

On sequential processors, stream reduction is trivially achieved through a single loop over the stream elements (see code in Figure 1). However, this sequential algorithm requires the ability to specify the position for writing memory access. For GPUs, stream reduction is a more complex task.

Although it is a fundamental element in many GPGPU applications, surprisingly little research has been published on stream reduction techniques. Horn [START_REF] Horn | Stream reduction operations for GPGPU applications[END_REF] was the first to propose an efficient stream reduction algorithm for the GPU. He performs stream reduction in two steps (see Figure 2): first, he computes the offset (displacement) for each non-empty element in the stream, using a parallel prefix sum scan [START_REF] Hillis | Data parallel algorithms[END_REF] Fig. 2. Main steps of Horn's stream reduction illustrated with an example. The unwanted elements are the crosses. The prefix sum scan computes the displacement, and the dichotomic search performs it (i.e. moves each valid element to the left at a distance given by the prefix sum). Fig. 3. The prefix sum scan used as the first part of Horn'algorithm has log n render passes. The first pass applies a boolean operation to each element of the input stream: if the element has to be discarded 1 is output, and 0 otherwise. Each subsequent pass performs the sum of two elements of the output of the previous pass. The prefix sum is the number of empty elements before each element, which is also the displacement that will be applied in the second part of the algorithm.

that counts the number of unwanted elements that are located before each element. Figure 3 details the different steps of the prefix sum scan. Second, each element is moved by the computed offset, using a logarithmic search. Both steps have O(n log n) complexity.

Sengupta et al. [START_REF] Sengupta | A work-efficient stepefficient prefix sum algorithm[END_REF] improved this algorithm by using a different method for the prefix sum scan [START_REF] Blelloch | Prefix sums and their applications[END_REF] of complexity O(n).The second pass (moving each element by the computed offset) is left unchanged, so the overall complexity of their algorithm is still O(n log n).

Recent GPUs, such as NVidia's G80, support geometry shaders, which can trivially be used to remove unwanted elements from a stream of input. All the elements of the stream are stored as vertices in a vertex buffer object, which is given as input to a geometry program. This geometry program tests each element and discards the unwanted ones. The fragment program is skipped and the resulting Fig. 4. We add a hierarchical step to the standard stream reduction algorithm: the input stream is divided in blocks, then each blocks is reduced, and finally the results are concatenated.

stream is directly output as a vertex buffer object, using the NV_transform_feedback OpenGL extension.

Sorting algorithms can trivially be used to perform stream reduction. The fastest GPU sorting algorithms have been based on bitonic sort [START_REF] Govindaraju | A cache-efficient sorting algorithm for database and data mining computations using graphics processors[END_REF]- [START_REF] Greß | GPU-ABiSort: Optimal parallel sorting on stream architectures[END_REF].

Our theoretical and experimental studies show that our stream reduction method is much faster than Horn [START_REF] Horn | Stream reduction operations for GPGPU applications[END_REF] and Sengupta et al. [START_REF] Sengupta | A work-efficient stepefficient prefix sum algorithm[END_REF]. It is also significantly faster than the naive approach using geometry shaders, especially on large streams. Finally, because they address a simpler problem, all dedicated stream reduction methods outperform sorting algorithms by an order of magnitude.

In CUDA, fast stream reduction is possible, by adding a scatter step to an efficient prefix sum scan implementation, such as Sengupta et al. [START_REF] Sengupta | Scan primitives for GPU computing[END_REF] or Harris et al. [START_REF] Harris | Parallel prefix sum (scan) with CUDA[END_REF]. However, it requires scatter abilities. We will show how our algorithm could benefit from scattering too.

III. ALGORITHM

In this section, we present our stream-reduction algorithm. We begin by a presentation of the algorithm (section III-A), then we discuss specific details for each pass, in sections III-B and III-C. The overall complexity of our algorithm is summarized in section III-D.

A. Overview

Our algorithm builds upon existing stream reduction algorithms, using a divide-and-conquer approach. First, we divide the input stream into smaller blocks of size s, and we perform a stream reduction pass on each of these blocks. Then, we concatenate the results of the stream reduction pass, using line drawing (see Figure 4). This hierarchical approach reduce the total number of operation and yields to a better asymptotic complexity. The blocks can be concatenated by openGL line drawing, and thus the algorithm is efficient and easy to implement.

B. Stream-reduction pass on each block

Once we have split the input stream into blocks of constant size s, we run a stream reduction pass on each of these blocks. We use recent stream reduction algorithms, running on the GPU [START_REF] Horn | Stream reduction operations for GPGPU applications[END_REF], [START_REF] Sengupta | A work-efficient stepefficient prefix sum algorithm[END_REF], in two passes: a prefix sum scan followed by a dichotomic search.

Here are the notations used in this section:

• x is one of the blocks.

• s is the size of the block.

• pSum is the prefix sum, an array of size s.

We will now detail the two passes. 1) Prefix sum scan: First, we run a parallel prefix sum scan (e.g [START_REF] Hillis | Data parallel algorithms[END_REF], [START_REF] Blelloch | Prefix sums and their applications[END_REF]) on each block. For every element in the block, this computes the number of empty elements before it, which is also the length of the displacement to apply to this element. The prefix sum is stored in pSum. This step has a complexity of O(s) using [START_REF] Blelloch | Prefix sums and their applications[END_REF].

2) Dichotomic search: Each element is moved by the computed offset. As GPUs lack scattering abilities from the fast fragment shaders, this is achieved through gathering: for each position i of the output (i.e. the reduced stream), a dichotomic search is performed in order to find the element x[j i ] of the input that has to be moved to i. This element is the only valid element of the input stream that satisfies

j i -pSum[j i ] = i.
As j -pSum[j] is monotonously increasing with j, it is possible to find j i by dichotomic search. The starting bounds of the research interval are given by:

i + pSum[i] ≤ j i ≤ i + pSum[s -1]
On older GPUs, all the loops are constrained to have the same length, thus the algorithmic complexity for each element is log s (complexity of the worst case), leading to s log s for the whole block.

With recent GPUs such as the GeForce 8800, it is possible to slightly improve this algorithm. First, the fixed length loop (for) can be replaced by a conditional loop (while), as branching is more efficient. Second, the bounds of the research interval can be refined at each step of the loop in order to reduce the number of iterations. j -pSum[j] is contracting with respect to j:

|j -j i | ≥ |j -pSum[j] -i|
Thus, at each step, the research interval is cut in half at a position j, and then further tightened by |j -pSum[j] -i|. The pseudocode of these improvements is given in Figure 5.

C. Concatenation of the intermediate results

After the reduction inside the blocks, we have n/s blocks, each of them with at most s non-empty elements grouped at its beginning.

For each block, the number of empty elements is known, as a side result of the first pass of stream compaction: it is the last element of the prefix sum. Those n/s last elements (one per block) form a sub-stream of the running sum. By running

1: lowerBound ← i + pSum[i] 2: upperBound ← i + pSum[s -1] 3: if upperBound > s -1 then 4:
There is no element at position i. Stop here. 5: j ← (lowerBound + upperBound)/2 6: f ound ← j -pSum[j] -i 7: while (f ound = 0 or x[j] is unwanted) do f ound ← j -pSum[j] -i 14: return x[j] Fig. 5. Improved dichotomic search at the i th position in a block x of size s: the program returns x[j i ] such as j i -pSum[j i ] = i. A conditionnal loop is used (line 7), and the bounds of the research interval are improved at each step of the loop (lines 9 and 11). The runing sum, previously computed, is stored in the pSum array. Fig. 6. In this example the block size s is 16. After the first part of the algorithm (reduction in the blocks), the last terms of the partial prefix sums, represented in grey color, are the number of empty elements in the blocks. These terms (except the last) form a stream of n/s elements. By applying a new prefix sum scan, the displacement of each block is computed.

a second pass of sum scan on this sub-stream, we get a new stream describing the number of empty elements before each block, which is also the length of the displacement that has to be applied to the blocks, as seen on Figure 6. The complexity of this process is O( n s ) using [START_REF] Blelloch | Prefix sums and their applications[END_REF]. Each block is then moved: the positions of the two ends of each segments are computed with the vertex shader using scatter, and the positions of the intermediate elements are linearly interpolated by rasterization. The vertex shader processes 2 n/s elements, and the fragment shader processes all the elements. Thus the number of operations required decreases with block size: using too small blocks stresses the vertex engine and implies more memory transfers due to the large number of segments. The complexity of these displacements is O(n).

As we store the streams in 2D textures, we have to convert the 1D stream into a 2D array, and thus we have to perform the wrapping of the segments, as illustrated on Figure 7. This can be done by sending all the segments twice: the first time draws the left part (or the full segment if it does not need wrapping), and the second draws the right part if needed or discards it otherwise. As most of the segments do not need such wrapping, it is better to use the geometry engine (if available) to split only when necessary.

D. Overall complexity

A summary of the algorithmic complexities of the different steps is given here:

1) Reduction of the blocks For one block:

• Prefix sum scan using [START_REF] Blelloch | Prefix sums and their applications[END_REF]: O(s)

• dichotomic search in the block: O(s log s) Total for the n s blocks: O(n + n log s) 2) Concatenation of the line segments

• Prefix sum scan using [START_REF] Blelloch | Prefix sums and their applications[END_REF]:

O( n s ) • Line drawing: O(n)
The asymptotic complexity of our algorithm is O(n + n log s) (where s is a constant), and is to be compared with the O(n log n) complexity of the previous methods.

IV. PREFIX SUM SCAN TECHNIQUES REVIEW

Our algorithm performs prefix sums computations at two different scales: inside the blocks and then across the blocks. Any technique can be used to compute those prefix sums.

We implemented both Hillis and Steele [START_REF] Hillis | Data parallel algorithms[END_REF] and Blelloch [START_REF] Blelloch | Prefix sums and their applications[END_REF]. In our experiments, despite its better asymptotic complexity, [START_REF] Blelloch | Prefix sums and their applications[END_REF] has a higher overhead. [START_REF] Hillis | Data parallel algorithms[END_REF] yielded better results for streams of 1 M elements and smaller, whereas [START_REF] Blelloch | Prefix sums and their applications[END_REF] is slightly faster for 4 M and 16 M elements. However, the relative performance difference was less than 10%.

Sengupta et al. hybrid algorithm [START_REF] Sengupta | A work-efficient stepefficient prefix sum algorithm[END_REF], which is a combination of Blelloch [START_REF] Blelloch | Prefix sums and their applications[END_REF] and Hillis and Steele [START_REF] Hillis | Data parallel algorithms[END_REF] method, could also be used for theoritically even faster results, although we have not tested this. However we would not expect a spectacular speed up, as we keep working on relatively small streams (s = 64) whereas [START_REF] Sengupta | A work-efficient stepefficient prefix sum algorithm[END_REF] is tailored for large ones.

V. RESULTS

A. Behavior

All the results of these section have been measured using a GeForce 8800 GTS, and an Intel Pentium IV 3 GHz with 2 Gb RAM. The algorithm was implemented using openGL. As shown in section III-D, the asymptotic complexity of our algorithm is linear with respect to the number of elements in the stream. This is confirmed by our experiments (Figure 8).

The execution time also depends on the ratio of valid elements: it is faster when either most of the elements are valid or most are unwanted. The maximum time is reached around 70% of valid elements (Figure 9). However these variations are relatively small. If not stated otherwise, all the timings in this document are given for a 60% ratio of valid elements, which is a difficult case for our algorithm.

The size of the blocks is an obvious parameter for our algorithm. Small blocks are reduced faster (complexity O(n + n log s)) but are concatened slower because of the complexity O(n + n s ) and the stress on the vertex engine. On the other hand, large blocks are reduced slower but concatened faster. Despite the lesser efficiency of the vertex engine, line drawing is not a bottleneck of the algorithm, except for very small block sizes (8 or less). Figure 10 shows the time repartition between the different steps of the algorithm for several block sizes, and Figure 11 shows the influence of the block size. We found experimentally that the best block size is 64 for streams up to 16 M elements. GeForce 8800 specific functionalities allows the two optimizations mentioned previously: single line drawing (section III-C) and improved dichotomic search (section III-B2). These optimizations save approximately 20% of execution time, as shown in Figure 12.

B. Comparison with other stream-reduction methods

We implemented the standard Horn's algorithm [START_REF] Horn | Stream reduction operations for GPGPU applications[END_REF], and a modified version that uses Blelloch's prefix sum scan [START_REF] Blelloch | Prefix sums and their applications[END_REF].

We also compared our algorithm to straightforward packing using GeForce 8800 advanced functionalities: elements are sent as vertices in a vertex buffer object, and the geometry shader discards the ones.

The comparison of these methods with ours summarized in Figure 13 and Table I. Our algorithm is 9 times faster than Horn's on 4M elements stream and 2 times faster than using geometry shaders. Notice how the speed-up ratio increases with the stream size.

In CUDA, scatter abilities are available. Sengupta et al. sum scan-algorithm, based on [START_REF] Blelloch | Prefix sums and their applications[END_REF]. Thus stream reduction can be efficiently implemented in this framework by adding a fast scattering step to [START_REF] Sengupta | Scan primitives for GPU computing[END_REF] or [START_REF] Harris | Parallel prefix sum (scan) with CUDA[END_REF]. We have reported the timings of [START_REF] Harris | Parallel prefix sum (scan) with CUDA[END_REF] in Table I: our algorithm is ≈2.5 X slower, however this CUDA algorithm requires scatter abilities whereas ours does not. We did not implement Sengupta et al. [START_REF] Sengupta | A work-efficient stepefficient prefix sum algorithm[END_REF], which achieves a 4 times speed up compared to Horn on streams of 1 M elements. As shown in Table I, our speed up is 7.8.

The largest stream on which we could run our algorithm had 16 M elements. We encountered the same limit with Blelloch [START_REF] Blelloch | Prefix sums and their applications[END_REF]. We could not make Horn [START_REF] Horn | Stream reduction operations for GPGPU applications[END_REF] and Geometry shaders work on streams with more than 4 M elements. In most cases, the limit is related to memory issues: such very large streams are taking a significant portion of the memory available on the card.

VI. CONCLUSION

In this paper, we have presented a new algorithm for stream reduction, which is an essential step of many GPGPU applications [START_REF] Horn | Stream reduction operations for GPGPU applications[END_REF], [START_REF] Roger | Whitted ray-tracing for dynamic scenes using a ray-space hierarchy on the GPU[END_REF]. OUR ALGORITHM IS ≈2.5 X SLOWER THAN CUDA, WHICH USES SCATTER ABILITIES WHEREAS WE DO NOT.

Our algorithm works by a hierarchical approach: we divide the input stream into smaller blocks, perform a fast stream reduction pass on these smaller blocks, and concatenate the results. We thus achieve a new order of asymptotic complexity (O(n) whereas previous methods where O(n log n)) and an impressive speed up (Table I). Our algorithm even outperforms doing stream reduction using the geometry shaders.

As our algorithm can use any method to perform the prefix sum, it can benefit from research in this domain.

Furthermore, our work can be seen as a meta-algorithm that improves any stream reduction algorithm by adding a hierarchical step to it.

Our divide-and-conquer approach results in impressive speed-ups for stream-reduction. In the future, we expect that this approach can be used for improving other GPGPU algorithms.

As a future work, we would like to implement our algorithm in CUDA. This would allow to speed up several steps: the reduction of the blocks could be computed sequentially using the algorithm Figure 1 (instead of a scan and a dichotomic search), line drawing and the line wrapping would disappear. The complexity would be then O(n), and could be comprared more accurately to existing methods in CUDA.

Fig. 1 .

 1 Fig.1. On sequential processors, the stream reduction of an array x of n elements is trivially done with a single loop.

Fig. 7 .

 7 Fig.7. Since we store the stream in a 2D texture, we have to wrap the line segments as we concatenate them. Segments 5, 16 and 20 have to be split, either by drawing them two times or in the geometry engine (if available).

Fig. 8 .Fig. 9 .

 89 Fig.8. Experimentally, our algorithm has a linear complexity with a slope depending on block size.

Fig. 10 .

 10 Fig. 10. Time spent in the different steps of the algorithm, for 4 M elements and various block sizes. Time spent in line drawing and global sum decreases with block size, whereas dichotomic search time and sum inside blocks time increases.

Fig. 11 .

 11 Fig. 11. Behavior of the algorithm with respect to block size. The best block size was 64 in all our experiments.

Fig. 12 .Fig. 13 .

 1213 Fig. 12. Effect of two optimizations requiring a recent GPU such as the GeForce 8800: single line drawing and improved dichotomic search.

TABLE I SPEED

 I UP OF OUR ALGORITHM COMPARED TO OTHER STREAM-REDUCTION METHODS (i.e. RATIO OF THE EXECUTION TIMES).

	#elements	Horn	Blelloch	G. Shaders	CUDA (sum only)
	256k	5.5	3.9	1.5	0.32
	1M	7.8	3.7	1.9	0.37
	4M	9.1	3.6	2	0.36
	16M	N.A.	4.8	N.A.	0.35