Gemma C Garriga

Esa Junttila

Heikki Mannila

Banded Structure in Binary Matrices ⋆

A binary matrix has a banded structure if both rows and columns can be permuted so that the non-zero entries exhibit a staircase pattern of overlapping rows. The concept of banded matrices has its origins in numerical analysis, where entries can be viewed as descriptions between the problem variables; the bandedness corresponds to variables that are coupled over short distances. Banded data occurs also in other applications, for example in the physical mapping problem of the human genome, in paleontological data, in network data and in the discovery of overlapping communities without cycles. We study the banded structure of binary matrices, give a formal definition of the concept and discuss its theoretical properties. We consider the algorithmic problems of computing how far a matrix is from being banded, and of finding a good submatrix of the original data that exhibits approximate bandedness. Finally, we show by experiments on real data from ecology and other applications the usefulness of the concept. Our results reveal that bands exist in real datasets and that the final obtained orderings of rows and columns have natural interpretations.

⋆ A preliminary version of this paper appeared in the proceedings of SigKDD 2008 (see [16]). This current submission corresponds to a significant extension which (1) includes proofs of formal statements; (2) presents new algorithms that outperform the ones presented in the previous paper and lift previous fixed-column permutation

Introduction

Matrices with binary values occur in many different applications. A typical example is market basket data gathered by retail companies [START_REF] Agrawal | Mining association rules between sets of items in large databases[END_REF]. Further than this, binary matrices abound in a large variety of fields ranging from information retrieval (documents and words occurrences) [START_REF] Baeza-Yates | Modern Information Retrieval[END_REF], to bioinformatics and computational biology (genes and probes mappings) [START_REF] Alizadeh | Physical mapping of chromosomes: A combinatorial problem in molecular biology[END_REF][START_REF] Myllykangas | DNA copy number amplification profiling of human neoplasms[END_REF], or ecology and paleontology (sites and species occurrences) [START_REF] Atkins | A spectral algorithm for seriation and the consecutive ones problem[END_REF][START_REF] Puolamäki | Seriation in paleontological data using Markov chain Monte Carlo methods[END_REF]. Understanding the properties of such matrices is therefore important for many applications. A fundamental problem is to uncover structures that will reveal the nature of the relations between the rows and the columns of the binary dataset.

In this paper we study the banded structure of binary matrices. A binary matrix is fully banded if both rows and columns can be permuted so that the non-zero entries exhibit a staircase pattern of overlapping rows. See Figure 1 for an illustration of a fully banded matrix. The concept of banded matrices has its origins in numerical analysis, where the matrix entries indicate connections between variables. From the computational point of view, working with banded matrices is always preferable: The work load involved in performing certain operations, such as multiplication, falls significantly for banded matrices [START_REF] Cuthill | Reducing the bandwidth of sparse symmetric matrices[END_REF], often leading to huge savings in terms of calculation time and complexity. There has been much research focused on minimizing the bandwidth of a matrix (broadly, the distance of non-zero entries from the main diagonal of the matrix) by applying permutations on the original matrix [START_REF] Aykanat | Permuting sparse rectangular matrices into block-diagonal form[END_REF][START_REF] Cuthill | Reducing the bandwidth of sparse symmetric matrices[END_REF][START_REF] Rosen | Matrix bandwidth minimization[END_REF].

From the data analysis perspective banded structures can occur in many applications. Consider for example the physical mapping problem of the human genome. Genome biologists break the genome into pieces (clones) which by recursive breaking can be eventually sequenced together. Unfortunately, the information about the relative positions of clones is lost during the breaking process. The physical mapping process starts with the experimental data from which information about the clone overlaps can be derived. The biological community has invested considerable efforts in the analysis of clone-probe matrices, in order to determine useful properties of both clone and probe orderings [START_REF] Alizadeh | Physical mapping of chromosomes: A combinatorial problem in molecular biology[END_REF].

For another application consider the presence/absence data from paleontology. Rows represent sites and columns represent species. A banded structure signifies an overlapping pattern between a set of species occurring in a spatially correlated set of sites. Or similarly, consider dialect word data described by a binary matrix of words used in several locations or municipalities. For this linguistic application a band provides a comprehensive visualization of the spatial distribution of dialects across the different municipalities of a country.

Another application where bands are potentially visible is in the discovery of overlapping communities without cycles in network data [START_REF] Banerjee | Model-based overlapping clustering[END_REF]. After finding a suitable permutation of rows and columns, a band of 1s from the adjacency matrix of a network should reveal communities of nodes that are strongly connected in a overlapping fashion. For example, consider the Football network dataset [START_REF] Girvan | Community structure in social and biological networks[END_REF] containing a match graph of football teams in US colleges in year 2000: a band corresponds to teams that frequently play together and are geographically close to each other. Or the network of characters in the novel Les Misérables [START_REF] Knuth | The Stanford GraphBase: a platform for combinatorial computing[END_REF] by Victor Hugo: a band shows the coappearence patterns by groups in the novel. Notice however, that communities in networks define smaller groups of highly interacting components and they cannot be always mapped in the shape of a band. The idea of bandedness is that the locality structure of the communities can be mapped in a overlapping fashion if this exists.

Fully banded matrices are not expected to arise in a real noisy environment. Therefore we study the problem of determining the minimum number of transformations one needs to do on the original binary matrix to uncover a banded structure. The number of such transformations will measure how far a matrix is from being banded. A simple example is shown in Figure 2. In its original form the matrix seems to be random, yet when permuted suitably, it exhibits a high concentration of 1s confined close to the main diagonal band. Fig. 2. An example of a binary matrix in its original form (a), and after permuting rows and columns to uncover its banded structure (b).

To better illustrate this optimization problem let us examine in more detail the network data applications mentioned above. Both Football and Les Misérables datasets are small matrices of sizes 115 × 115 and 77 × 77 respectively. The approximate bands exhibited by the datasets are shown in Figure 3. The band shown by the Football dataset is only 534 bit flips away (that is, number of 0-to-1 or 1-to-0 transformations) from being fully banded, while the band exhibited by Les Misérables only requires 201 flips. In both cases the ordering of rows and columns reveals clusters that are linked in an overlapping fashion. For the football dataset we observe a clear organization of teams in small clusters; indeed, each of those non-overlapping communities represent teams that are geographically close to each other and thus, they usually play one against each other. For Les Misérables dataset we can clearly see which group of characters co-occur together in the novel, and on the other hand, which characters appear more independently across the chapters. Not all network datasets might contain bands, but still they might contain communities. The notion of a community expresses a group of highly interacting nodes, which is more general. The concept of bandedness we study here forces the mapping of the discovered communities in a overlapping fashion if this is possible. These simple examples already reflect the complexity involved in this problem. Another optimization problem we consider to cope with the unavoidable noise is finding a maximum subset of rows and columns from the original matrix that will exhibit an almost banded structure. Here we explore the combinatorial properties related to bandedness and propose algorithmic solutions to solve the associated optimization problems. Finally, we demonstrate the usefulness of the band-concept by applying our methods to an extensive collection of synthetic datasets and several real life-sciences applications: paleontological dataset, mammal dataset, DNA amplification data and Finnish dialect word data.

The rest of this paper is organized as follows. In Section 2 we give a formal definition of the concept of bandedness and study its theoretical properties. Section 3 introduces the definition of the banded optimization problems. In Section 4 and Section 5 we propose algorithmic solutions for the banded augmentation problem; and in Section 6 for the banded submatrix problem. Finally Section 7 shows empirical results and Section 8 discusses related work.

Bandedness and its combinatorial properties

This section studies the combinatorial nature of banded matrices. We introduce an incompatibility graph between the columns of the matrix that will characterize the proper permutations for obtaining a fully banded structure. We also define an intuitive binary relation between the rows which, under those proper column permutations, will lead to a banded structure. The combination of these two results drive to a polynomial test for bandedness.

Bandedness

Consider a n × m binary matrix M . In the genetic fingerprinting application, rows would correspond to probes and columns to clones; in the ecological application, rows would be sites and columns would be species. We denote the i-th row of M by M i and the j-th column of M by M j . Given a permutation π of natural numbers {1, . . . , m} and a row vector M i , we denote with M π i the permutation of the vector M i induced by π. In general M π represents the matrix resulting from applying the column permutation π to each one of the rows of M . Similarly, M π κ represents the permutation of rows and columns according to κ and π respectively. Sometimes it is useful to interpret the row and column vectors of M as sets of indices, that is a row M i is as well a set of column indices that appear in the row.

Roughly speaking, a matrix is fully banded if both rows and columns can be permuted in a way that the non-zero entries exhibit a staircase pattern of overlapping rows (see Figure 1). A formal definition follows.

Definition 1.

A binary matrix M is fully banded if there exists a permutation of rows κ and a permutation of columns π such that (1) for every row i in M π κ the entries with 1s occur in consecutive column indices {a i , a i + 1, . . . , b i } and (2) these indices satisfy a i ≤ a i+1 and b i ≤ b i+1 .

From the definition we observe that bandedness is a hereditary property: submatrices of a fully banded matrix preserve that property. Also, if M is fully banded also its transpose, M T , is fully banded. Another observation is that for a matrix to be banded it has to satisfy at least the consecutive-ones property on the rows and the columns. Formally, a binary matrix M is a consecutive-ones matrix if it is possible to order the columns so that, in every row, the non-zero entries occur in consecutive positions [START_REF] Tucker | A structure theorem for the consecutive 1's property[END_REF]. A consecutive interval of 1s from a row will be denoted as [a, b]: all the entries between columns a and b (inclusive) are 1s and the rest 0s. Unless explicitly stated, when a binary matrix M is said to satisfy the consecutive-ones property we assume that it comes permuted to exhibit this consecutive arrangement of ones in all rows. Thus, we can omit permutation π from the notation. Testing the consecutive-ones property of a binary matrix and representing permutations π (if they exist) can be done in polynomial time via PQ-trees [START_REF] Booth | PQ-tree algorithms[END_REF], and more recently, in linear time with a certifying algorithm [START_REF] Mcconnell | A certifying algorithm for the consecutive-ones property[END_REF]; unfortunately, these algorithms do not handle noise.

Not all matrices with the consecutive-ones property are banded, i.e., part (2) in the definition is not superfluous. For example, the following is a consecutiveones matrix for the rows and the columns yet none of the row-column permutations leads to a staircase pattern of ones.

  1 1 1 0 0 1 1 1 0 0 1 0   (1)
On the other hand, even if the matrix is fully banded, we cannot expect that all column permutations from the consecutive-ones property will lead to an overlapping sequence of rows. An example is the following matrix, which as shown in Equation 2 has consecutive ones for the rows but it is only banded when the last column is placed first.

  0 0 0 1 1 1 1 1 0 1 1 0   (2)
Next we characterize exactly the relation of banded matrices with consecutiveones matrices. This will be possible via the following binary relation between rows.

Definition 2. Let M i = [a, b] and M j = [a ′ , b ′
] be two rows with the consecutiveones property. We say that M j is properly included in M i , denoted M j ≺ M i , if and only if a < a ′ and b ′ < b. We denote by M i ⌢ M j whenever M i ⊀ M j and M j ⊀ M i .

As an example, the third row in the small matrix of Equation 2 is properly included in the second row. The rows of a consecutive-ones matrix M are a Sperner family of intervals if for any two rows M i , M j we have M i ⌢ M j . This Sperner property on the family of row intervals can be seen as a restricted version of the Sperner property on the family of row sets: two rows M i and M j with M i ⊂ M j might, or might not satisfy M π i ⊀ M π j . This depends on the column permutation π of the consecutive-ones arrangement for M , which, at the end, determines the starting and ending points of the row intervals. For example, the third row in Equation 2would not be properly included in the second row if the last column was placed first. The following statement characterizes exactly this relation.

Lemma 1.

A binary matrix M is fully banded if and only if M is a consecutiveones matrix for a permutation of columns π where every two rows i and j satisfy M π j ⌢ M π i .

Proof. If for every two rows M π i = [a i , b i ] and M π j = [a j , b j ] with consecutive ones we have M π i ⌢ M π j , then: either a i ≤ a j and b i ≤ b j , or a j ≤ a i and b j ≤ b i . Under these conditions we can easily establish a preorder of rows by sorting them in ascending value of a i 's, while resolving ties with the ascending value of their b i 's; if two rows are exactly the same, then ties between them can be resolved arbitrarily. This preorder defines directly the banded structure of Definition 1. The other direction of the statement is trivially implied from the definition of fully banded matrices. ⊓ ⊔

An incompatibility graph for bandedness

The set of orderings that satisfies bandedness is a subset of those satisfying consecutive-ones, as argued with Lemma 1. Unfortunately, checking whether a matrix is fully banded would be unfeasible if we had to go through all the permutations of the consecutive-ones property of the same matrix: there might be exponentially many such permutations. The key lies in identifying the difference between bandedness and consecutive-ones. In the remaining of this section we will study the incompatibilities that arise under certain permutations of the columns to derive a polynomial test for bandedness.

A very convenient data structure to deal with the valid column orderings of a consecutive-ones relation is the incompatibility graph, as introduced in [START_REF] Mcconnell | A certifying algorithm for the consecutive-ones property[END_REF]. An incompatibility graph is simply an undirected graph whose vertices are pairs of elements (a, b), for all column indices a = b, and whose edges reflect the incompatibility relation derived from the following observations. Suppose we want to discover the consecutive-ones ordering of a set of columns. For each pair of columns a, b of M , we have that the order (a, b) (column a before b) is incompatible with order (b, a). Moreover, for a triplet of columns a, b, c and a row i of M with a, c ∈ M i and b / ∈ M i , the pairs (a, b) and (b, c) cannot appear in the same consecutive-ones relation. Namely, they would imply that b goes between a and c in the column ordering, which would imply that M fails to be consecutively ordered in the i-th row. Therefore, (a, b) and (b, c) are incompatible in this case. It follows from definition that only these types can violate consecutive-ones property. The incompatibility graph G(M ) of M is an undirected graph formed by adding all edges of these types. A toy example of an incompatibility graph is shown in Figure 4. It is then clear that a consecutive-ones relation must not contain incompatible pairs, thus giving the following result.

Proposition 1 ([24]

). A binary matrix M has the consecutive-ones property if and only if its incompatibility graph G(M ) is bipartite.

Next we study further the observation given by Lemma 1 and use the incompatibility graph to test the bandedness property in polynomial time. The idea is to augment the graph G(M ) with more incompatibility edges; basically, those incompatibility ordering pairs that restrict the consecutive-ones orderings to those preserving a Sperner family of row intervals. We denote the augmented graph with Ĝ(M ).

Suppose we have two rows i and j of M and M i ⊂ M j . Let a, b, c ∈ M j be column indices such that a ∈ M i but b, c / ∈ M i . Then if the matrix is fully banded, the orderings (b, a) and (a, c) cannot both hold. Namely, if b is before a and a is before c, we have entries as follows. row b a c M i : 0 1 0

M j : 1 1 1 (3)
But in this case, the ordering cannot satisfy the conditions of a full band as one interval is included in the other (as seen in Lemma 1 this is a necessary condition for bandedness). We call these new edges added to G(M ) to form Ĝ(M ) the set of Sperner conflicting pairs. As an example consider the matrix in Equation 1: it has the consecutive-ones property, as shown in the incompatibility graph of Figure 4. However, when adding to that graph the incompatible pairs from M 1 \M 3 = {a, b} and M 2 \M 3 = {b, d}, the final graph would not be bipartite anymore. This implies that the matrix in Equation 1 is not banded in any of the consecutiveones permutation.

Adding Sperner conflicting pairs to graph Ĝ(M ) restricts the set of consecutiveones orderings to fully banded orderings. To see this, recall that part (2) in Definition 1 separates bandedness from consecutive-ones. The only way for matrix M to violate part (2) is that the matrix in Equation 3 is a submatrix of column-ordered M . The new edges prevent this. Notice that by construction the incompatibility graph is complete: for each edge there is an incompatiblity pair (due to part (1) or part (2) of Definition 1) and for every incompatibility pair there will be an edge reflecting such conflict in the graph.

The number of Sperner conflicting pairs that can be potentially added to the graph Ĝ(M ) is at most quadratic in n. A way of testing bandedness in linear time in the number of edges of the incompatibility graph is then the following. Proposition 2. A matrix M is fully banded if and only if its incompatibility graph augmented with the set of Sperner conflicting pairs is bipartite.

Proof. It follows directly from the construction of the incompatibility graph and the property stated by Lemma 1. In the incompatibility graph an ordering is defined by a set of nodes. An ordering of columns satisfying the fully banded property must have no incompatible pairs in the graph augmented with the Sperner conflicting pairs. Thus, a proper column ordering for the fully banded property is an independent set in this graph that consists of half of the vertices.

The reverse of a fully-banded ordering is also a fully-banded ordering, so the remaining vertices of the incompatibility graph must be also an independent set. Therefore, the incompatibility graph must be bipartite whenever the matrix is fully-banded. This is also a sufficient condition: the incompatibility graph fails to be bipartite whenever the matrix fails to be fully-banded. It follows that if the graph is not bipartite then it must have an odd cycle, and therefore there is no independent set that contains half of the vertices in the graph, implying that there is no column ordering that would make the matrix fully banded.

⊓ ⊔

We derive next a generalization from Proposition 2. The advantage is that it will summarize the conflicting Sperner pairs directly into an augmentation of M . This will be particularly important for the algorithmic solutions we later present to our optimization problems.

Lemma 2. Let M be the binary matrix M augmented with a set of new rows M ij = M j \M i for every two rows M i ⊂ M j . Then we have that M is fully banded if and only if M has the consecutive-ones property.

Proof. If any two rows i and j of M such that M i ⊂ M j satisfy M i ⌢ M j under a certain permutation π of the consecutive-ones, then M j \M i will be consecutiveones for that permutation π as well. On the other direction: any permutation π of the consecutive-ones on M , preserves both the consecutive-ones property of M and the Sperner family of intervals of their rows.

⊓ ⊔

The maximum number of rows that can be potentially added to the original

M is at most n • (n -1)/2.

Problem definitions

Real world matrices are not expected to be fully banded. Therefore we introduce a measure that computes how far a matrix is from being banded. We will do so by looking at the minimum number of 0s that need to be transformed into 1s to make a matrix banded. We name this optimization problem as the Minimum Banded Augmentation problem. Problem 1. (Minimum Banded Augmentation, Mba) Given a binary matrix M , find the minimum number of 0s that need to be transformed into 1s so that M becomes fully banded.

The formulation of the problem above offers a way to estimate how strong the underlying banded pattern is; this will be closely related to the permutations of rows and columns needed to find the full band. In addition, if a dataset contains a banded pattern, we can identify entries that are likely errors by checking which entries were transformed. In the best case, the fully banded matrix induced by transformations is closer to reality than original erroneous data is. From now on, the transformations are also called flips, and the solutions to the Mba problem are called minimum-flips values.

Example 1. Finding the minimum number of 0s that need to be transformed into 1s in order to make the matrix band is not easy. Consider the following input matrix and the transformation.

  1 0 1 0 1 1 1 1 0 1 1 0     1 0 1 0 1 1 1 1 0 1 1 1  
Despite the appearances, only one flip (bottom right) is enough to transform the matrix into a full band. Notice that after this flip, the third of the columns in the matrix can be placed in the second position to display the full band.

We denote with β(M ) the optimal number of flips to make the matrix fully banded. Matrices that are fully banded satisfy β(M ) = 0. A trivial observation following from the symmetry of our definition of fully banded matrices is the following.

Proposition 3. For a binary matrix M , we have β(M ) = β(M T ), where M T is the transpose of M .

For ecological and paleontological data we often have the situation where the 1s are reasonably certain, but the 0s can be missing values. Thus the Mba problem is quite natural. However, for other types of applications it is useful to allow flips in both directions, from 1 to 0 and from 0 to 1. We refer to the version of Mba as the Bidirectional Mba. We will denote with β B (M ) the optimal number of such bidirectional flips to obtain a fully banded matrix. Note that β B (M ) ≤ β(M ).

A variation of the Bidirectional Mba is to consider its weighted version: suppose we have a cost associated to 0-to-1 flips and a cost associated 1-to-0 flips; the weighted problem becomes that of finding a transformation with the minimum cost. This weighted variation of the Bidirectional Mba can be especially useful in ecological applications, as it will be justified in the experiments.

Example 2. To illustrate the importance of bidirectional flips, consider the following input matrix as in Example 1 with one extra row.

    1 0 1 0 1 1 1 1 0 1 1 0 0 1 1 0         1 0 1 0 1 1 1 0 0 1 1 0 0 1 1 0    
Even if we have one more row, only one flip is necessary to transform the band into a full band: in this case, the last 1 on the second row is flipped. Again, move the third column to the second position to see the full band.

The algorithmic complexity of the Mba problem is expected to be hard. The basis of this assumption arises from Lemma 2 and the known fact that the Consecutive-ones Matrix Augmentation problem (with 0-to-1 flips) is NP-complete [START_REF] Papadimitriou | The NP-completeness of the bandwidth minimization problem[END_REF]. Also symmetric banded matrices have an interpretation as proper interval graphs [START_REF] Roberts | Indifference graphs[END_REF], and general banded matrices as proper interval digraphs [START_REF] Sen | Indifference digraphs: A generalization of indifference graphs and semiorders[END_REF]. Both the edition and deletion of edges to transform a general graph into an interval graph is NP-complete [START_REF] Burzyn | NP-completeness results for edge modification problems[END_REF].

A second problem we study here is the following.

Problem 2. (Maximum Banded Submatrix, Mbs) Given a binary matrix M and an integer k, find the maximum submatrix M ′ of M such that it is banded after at most k flips.

The Mbs problem will be useful in datasets exhibiting several independent band structures, or also, when noise is too high to identify a clear band from the complete dataset. In essence, a solution extracts the most relevant banded pattern in a dataset. Since Mbs is only a generalization of the Mba and Bidirectional Mba problem, the algorithmic complexity is also expected to be hard.

Algorithms with fixed column permutation

Since the values of β(M ) and β B (M ) cannot be computed exactly for larger instances, we consider algorithms for their estimation in both the Mba and the Mbs problems. We denote those upper bounds with β(M ) and βB (M ).

We first investigate a special form of the problems where we assume that a column permutation that establishes the banded structure is already known. Later we will show ways to lift this requirement.

Algorithms for Mba

By Lemma 1, when having a fixed column permutation the Mba problem decomposes into enforcing the consecutive-ones property for single rows, and after that, making some final flips to ensure that rows form a Sperner family of intervals. Algorithm 1 outlines this Fixed Permutation solution.

Assuming that the column permutation π is given and that Mba only allows 0-to-1 flips, there is only one exact way of producing a consecutive-ones relation: by flipping all possible 0s falling between 1s for each row of M π (lines 4-6). A second phase has to ensure that all row intervals will be pairwise overlapping, that is M π i ⌢ M π j for all rows i = j (lines 7-13). Since only 0-to-1 flips are allowed, the solution is simply to extend the row intervals. An extension of

M i = [a, b] means to update the endpoints of the interval for a new [a ′ , b ′ ] such that a ′ ≤ a ≤ b ≤ b ′ .
At every step the algorithm takes a row M π i and calculates its optimal extension (if this is needed). To do so it will select all the superintervals M π j ≻ M π i and check all the potential extensions for M π i that would resolve the Sperner conflicts that the i-th row has with those j-th rows.

An extension of M π i that will always resolve all Sperner conflicts for that row can either be a left-hand side extension to the leftmost M π j ≻ M π i (line 10 (a)); a Algorithm 1 The Fixed Permutation algorithm for Mba 1: Input: An n × m binary matrix M and a permutation π 2: Output: A permutation κ of rows 3: Fix the column permutation of M to be π 4: for each row i in M π do 5: Flip all those 0s falling between 1s 6: end for 7:

for each row i in M π s.t. M π i = [a, b] do 8: Let C = {M π j = [aj , bj ] | M π i ≺ M π j } 9:
Extend right-hand side extension to the rightmost M π j ≻ M π i (line 11 (b)); or, extending M π i to both left and right-hand sides with a combination of two superintervals (line 12 (c)). This only requires checking the starting point of each M π j ≻ M π i and combining it with the right-most ending point from all other superintervals M π k ≻ M π i whose starting point comes before M π j . Eventually, for a row i the algorithm takes the extension that represents fewest transformations.

M π i = [x,
Resolving the Sperner conflicts for a row i does not change the optimal extension of other rows j = i (there are no cascade effects). To see this suppose that i includes a new interval k after its extension, and let j be the interval that previously contained i and that originated the best optimal extension of i (in lines [START_REF] Burzyn | NP-completeness results for edge modification problems[END_REF][START_REF] Chen | Best k": critical clustering structures in categorical datasets[END_REF][START_REF] Cormen | Introduction to algorithms[END_REF]. This interval j contains k by transitivity, and indeed, it also contained k even before i was extended. This implies that when processing k, its optimal extension cannot be changed because of a previous extension of interval i.

Example 3. To illustrate the Fixed Permutation algorithm, consider again the small matrix in Example 1.

  1 0 1 0 1 1 1 1 0 1 1 0     1 1 1 0 1 1 1 1 0 1 1 0     1 1 1 0 1 1 1 1 0 1 1 1  
Under this given fixed column permutation, the algorithm would transform first all 0s falling between 1s: this corresponds to flipping one 0 in the first row of the matrix. Then, the algorithm would resolve the Sperner conflicts between the rows: there is a conflict between the second and the third row, which can be resolved by flipping the bottom right 0 into a 1. The matrix is fully banded for this permutation after doing two 0-to-1 flips. The number of flips done in this illustration is 2, while we showed in Example 1 that only one flip would be Algorithm 2 The Bidirectional Fixed Permutation algorithm for Bidirectional Mba 1: Input: An n × m binary matrix M and a permutation π 2: Output: A permutation κ of rows 3: Fix the column permutation of M to be π 4: for each row i in M π do 5:

Let W π i be the weight vector for row i of M 6:

Let if

M π i ⊂ M π j then 11: Let A = M π j \M π i 12:
Let WA be the weight vector for A 13:

Let [a, b] be the solution of the maximum consecutive subarray of WA 14:

Update

M π i preserving M π j \M π i = [a, b] 15:
end if 16: end for 17: Sort rows [a, b] of M π in an ascending order of as, while deciding ties with the ascending order of their bs.

necessary to make the matrix fully banded. Therefore, an important component of this fixed-permutation algorithm is the initial permutation that is given to the algorithm.

For a fixed permutation π of the columns, the Fixed Permutation algorithm computes exactly β(M π ) in polynomial time. This will be always an upper bound of β(M ). The number of comparisons between rows made by the algorithm will be at most n • (n -1)/2. Since every comparison requires at most n comparisons (line 12 (c)), the final complexity of the algorithm is of the order O(n 3 ).

Finally we have the following property derived from the above algorithm.

Proposition 4. If a permutation π of the columns is fixed, Mba for M π is solvable in polynomial time.

Algorithms for bidirectional Mba

Next we study the properties of the algorithms for the Bidirectional Mba problem. Again we will use the same principle of assuming a fixed column permutation π. The basic idea here is also simple: first solving optimally the consecutive-ones property for bidirectional flips on the M π and after that, resolving Sperner conflicts between rows. Algorithm 2 gives the outline of the Bidirectional Fixed Permutation algorithm.

The first simple observation we use is the following.

Proposition 5. If a permutation π of the columns is fixed, the minimum number of bidirectional flips that lead to a consecutive-ones relation on M π can be solved exactly in linear time.

To see this consider a n × m matrix of weights W where each non-zero entry of M will be assigned a +1 in W , and each zero entry of M will be assigned a weight of -1. This kind of weight matrix assumes equal cost between the two types of flips.

W j i = +1 if M j i = 1 -1 if M j i = 0
In this form, the problem of finding the optimal consecutive-ones solution with bidirectional flips for M π i corresponds to solving the maximum subarray problem on W π i . The aim of the maximum subarray problem is, for a given array of numbers, to find a consecutive subarray such that the sum of the numbers in the subarray is maximum. This can be done in linear time in the size of the array by making use of a scan-line algorithm [START_REF] Cormen | Introduction to algorithms[END_REF]. In essence, for a solution [a, b] of the maximum subarray algorithm over W π i , we should update M π i with consecutive ones between the columns a and b, while setting to 0s the rest of its entries. This corresponds to lines 4-8 in Algorithm 2. Note that weights of the matrix W can be tuned according to the apriori knowledge of the application. In this way we would be solving the weighted Bidirectional Mba. As an example, after substituting +2 for +1, 0-to-1 flips have half the cost of 1-to-0 flips.

In the second phase (lines 9-16) the algorithm proceeds by removing the Sperner conflicts between the rows of M π , which at this stage already exhibit a consecutive-ones property. The useful technical observation comes from Lemma 2: A matrix M will be fully banded as long as the augmented M has consecutive ones. Therefore, to eliminate all possible Sperner conflicts between the row intervals of M π , the algorithm simply has to go through all the extra rows described in M and make them consecutive ones. When propagating the changes back to M π we would end up with a banded matrix. As above, we can make use of the maximum subarray problem on the extra rows of M to solve the problem exactly for bidirectional flips. It only remains to update rows in M π so that they are kept consistent with the changes made over M . The final obtained solution on M π will be always banded. Basically, this corresponds exactly to deciding the best number of bidirectional flips that will eliminate Sperner conflicts in a pairwise comparison of rows in M π . The final complexity of the algorithm is O(n 2 m).

The Bidirectional Fixed Permutation algorithm, however, is not optimal for the Bidirectional Mba problem on M π . The problem resides in the second phase: rows are compared in a pairwise fashion and globally beneficial updates may be missed. In general, it seems that the Bidirectional Mba problem would be still hard for a fixed permutation of the columns.

Example 4. To illustrate the Bidirectional Fixed Permutation algorithm consider the matrix in Example 2.

    1 0 1 0 1 1 1 1 0 1 1 0 0 1 1 0         1 1 1 0 1 1 1 1 0 1 1 0 0 1 1 0         1 1 1 0 1 1 1 0 0 1 1 0 0 1 1 0    
Under this given fix column permutation, the algorithm would transform first the matrix into a consecutive-ones relation; as shown before, this can be solved exactly by finding the consecutive maximum subarray problem for each row. In the example, this corresponds to flipping second 0 in the first row to a 1 (notice that we could have also flipped any of the 1s of the first row into a 0 with the same final cost). In the second phase, the algorithm resolves the conflicts between the rows: there are two conflicts between the second row and the third and fourth row respectively, which will be resolved after flipping the last 1 of the second row into a 0.

As a final step in both Algorithm 1 and Algorithm 2, it only remains to sort the rows with the preorder given by the row intervals to visually exhibit the fully banded structure. The banded pattern in original dataset is described by this row-ordering.

Finding good column permutations

An essential component of Algorithm 1 and Algorithm 2 is to identify a good column permutation that will be fixed since the beginning; new algorithms that lift this restriction will be introduced in the next section. Intuition says that a good permutation will tend to put similar columns close enough to each other in the order. If the matrix exhibits banded structure, an order based on the similarity between columns should also preserve the band. We define next several similarity measures to compare columns and use two different methods to find a proper order.

One measure of similarity between two columns is the correlation similarity. Given two columns M a and M b the correlation similarity is:

CorrS(M a , M b ) = (1 + ρ ab )/2,
where ρ ab represents the Pearson coefficient between the columns. Values range from 1 (identical columns) to 0 (anticorrelated columns).

As an alternative we will use the overlapping measure computed by the Jaccard coefficient:

J(M a , M b ) = |M a ∩ M b | |M a ∪ M b | .
This similarity measure captures the particularities of our problem where columns should be increasingly overlapping one with the other. It takes the highest value of 1 when columns M a and M b are exactly the same. The overlapping similarity for non-intersecting columns is 0. Other straightforward measures we will use are the dot product and the Hamming distance metric between columns.

In order to find a good permutation of the columns preserving the similarity or distance relationship, we consider a complete undirected graph whose nodes are the columns of the input matrix M . The weight of an edge {a, b} is defined by the similarity between a and b. This similarity can be defined by any of the measures mentioned above. We will use two methods to find a good order between Algorithm 3 The Alternating algorithm for Mba or Bidirectional Mba 1: Input: An n × m binary matrix M , a number of iterations t 2: Output: A permutation κ of rows and a permutation π of columns 3: Initialize π with a random permutation of the columns 4: Let A = M π 5: repeat 6:

Apply a fixed column permutation algorithm on A 7:

Transpose the current matrix A = A T 8: until reaching t iterations 9: Return the pair of row-column permutations, κ and π, with the best band found.

columns: spectral ordering [START_REF] Atkins | A spectral algorithm for seriation and the consecutive ones problem[END_REF] and approximation of the minimum Hamiltonian path (by means of constructing a minimum spanning tree) [START_REF] Cormen | Introduction to algorithms[END_REF]. Note that similarity measures will be transformed into distances for the Hamiltonian path.

Lifting the fixed column permutation requirement

Algorithm 1 and Algorithm 2 are dependent on fixing a column permutation beforehand. In this section we propose algorithmic strategies that avoid this dependency: two new methods (Alternating and Simulated Annealing) and a pre-existing method (Barycentric), now adapted to solve band-problems.

Alternating method

The first of these strategies follows from Proposition 3. It means that in practice we can solve the Mba problem or Bidirectional Mba problem on either M or M T . This suggests the following alternating approach: solving the problem of finding a good row permutation given the current column permutation (with any of the fixed permutation algorithms 1 or 2) and then, transposing the matrix to iterate with the current configuration until convergence, or until a certain number of iterations is reached. The pseudo-code of this Alternating algorithm is shown in Algorithm 3.

Notice that the alternating strategy just described does not necessarily converge for all matrices. Indeed, it is not possible to guarantee that the value of β(M ) or βB (M )-the approximated value for β(M ) and β B (M ) respectivelywill decrease after each iteration. To stop the alternating process, we bound the number of iterations by an input parameter t. The output of this algorithm is the pair of row-column permutations that achieved the best value of β(M ) among the t iterations.

Barycentric method

Another strategy that transposes the matrix is presented by the Barycentric algorithm, originally proposed in [START_REF] Sugiyama | Methods for visual understanding of hierarchical system structures[END_REF]. In essence, the Barycentric algorithm finds good permutations of both rows and columns based on a barycenter measure (average position of 1s in a row). Let the barycenter of row i be defined as follows:

Barycenter(i) = j=1..m j • M j i j=1..m M j i .
The Barycentric algorithm first computes the barycenter for all rows, then it orders (stable ordering) the rows from smallest to largest barycenter, and finally, it transposes the matrix M to iterate again following the same strategy until convergence. Notice that this sorting process does not use flips or compute in any way the borders of the band at any iteration. Indeed, Barycentric only orders rows and columns according to the barycenter measure in an iterative fashion.

Simulated annealing method

Finally, we present a simulated annealing metaheuristic to find a good approximation of the global minimum for the Mba and Bidirectional Mba problems. We do not expect this stochastic method to be as competitive in terms of running time as the other algorithms presented this far. In general, the simulated annealing strategy comes with a high computational cost, but it yields nearoptimal results in most cases. For this reason we use this strategy as a reference method for comparing the different bandedness solutions. Simulated annealing [START_REF] Kirkpatrick | Optimization by simulated annealing[END_REF] is a general stochastic optimization method, used successfully to solve various combinatorial problems. In band-related problems, different row-column permutations of the matrix are interpreted as "states" and minimum-flips costs as "energy values". We seek an order of rows and columns for the matrix M that produces the smallest energy, that is β B (M ). The search space includes all permutations of rows and columns.

The simulated annealing algorithm starts with a random state and a large temperature parameter T , for example 10.0. The algorithm keeps track of the current state S of the system all the time. The temperature is decreased at each iteration of the algorithm until a pre-determined number of iterations has passed. At each iteration, a new candidate state S ′ is generated stochastically by function Neighbor(S). The new state S ′ is then accepted as a new current state with probability min{1, e ((E(S)-E(S ′ ))/T ) }, where E is the energy function whose minimum we seek. In our case, this energy function is the minimum number of flips needed to make the matrix fully banded, with the restriction that both current row and column permutations remain fixed. An optimal polynomial-time algorithm (Optimal Visual Band) for the problem is given shortly. Algorithm 4 shows the scheme of the simulated annealing strategy.

We next define the technical details for the auxiliary functions Neighbor and Optimal Visual Band for the simulated annealing algorithm.

Picking a neighbor

The most important requirement for the function Neighbor is that all states (in our case, permutations) in the search space should be ac- Pick a neighbor S ′ = Neighbor(S) 7:

Compute energy

E ′ = Optimal_Visual_Band(S ′ ) 8: if E ′ < E then 9:
Update best solution S best = S ′ and E best = E ′ 10:

end if 11:

Update S = S ′ and E = E ′ with probability min{1, e ((E-E ′ )/T ) } 12:

Decrease temperature T = T • α 13: until reaching t iterations 14: Return the row-column permutations of the best solution S best .

cessible by repeating the candidate generation process sequentially. We consider several candidate generation methods [START_REF] Tian | Application of the simulated annealing algorithm to the combinatorial optimisation problem with permutation property: An investigation of generation mechanism[END_REF] for finding good permutations and we apply them in our 2-dimensional permutation problem. For all these methods we treat an ordering as a cycle: after the last row the first row follows. We summarize the methods here.

-Swap-k: Choose two random rows and swap them; do the same for two random columns. Repeat k times. -Adj-swap-k: Randomly choose a row index r, then swap (adjacent) rows r and r + 1; do the same for a randomly chosen column. Repeat k times. -Reverse: Choose two random row indices r and r ′ . Starting from r, move forward in row-order until r ′ is found. Reverse the order of all encountered rows, including r and r ′ . Do the same for columns. -Relocate: Choose two random row indices r, r ′ , and a random integer k ∈ {1, . . . , n}. Starting from r, move forward in row-order until r ′ is found. Shift all encountered rows k steps forward in row-order. Do the same for columns. -Reverse+Relocate: Apply both reverse and relocate methods simultaneously: first choose a set of rows, then relocate the rows and reverse their order. Do the same for columns.

Finding the optimal visual band For Simulated Annealing to produce good results, we must compute energy values accurately and efficiently. To this end, we propose a polynomial-time exact algorithm Optimal Visual Band when both row and column permutations are fixed.

A fully banded matrix is visually banded, if the row-column permutations of the matrix conform to bandedness. The problem to be solved is as follows: given a binary matrix M with fixed row-column permutations and a non-negative weight matrix W , find the minimum cost of transformations needed to make the matrix visually banded.

Let a tuple of column indices s, e represent the interval [s, e -1] including columns positions s, . . . , e -1; the tuple of indices s, s represents an empty interval. The cost of an interval on a row is the total cost for transforming the row so that the values inside the interval are all 1s and other values are 0s. The goal is to define intervals on each row r in M so that the intervals of 1s across the rows have visual bandedness property and the total cost is minimum.

In order to simplify the algorithm, assume that the weight matrix W has been split in two separate matrices: W1 for the costs of transforming 1s and W0 for the costs of transforming 0s. If M c r = 1, then W1 c r = W c r and W0 c r = 0, otherwise W1 c r = 0 and W0 c r = W c r . The algorithm, which is based on dynamic programming, goes through all rows one by one. It ensures that the 3-dimensional array C contains optimal cumulative costs for intervals on the previous rows. Using array C, it evaluates the cost of each interval on current row with respect to M and W and stores the cumulative interval cost in C. Visual bandedness of intervals is preserved at all times. For each triplet of indices r, s, e, varying along the row indices for r and along the column indices for s and e, we update array C(r, e, s ) via a classical recurrence: The recurrence is only applied to triplets of indices r, s, e that satisfy 1 ≤ r ≤ n and 1 ≤ s ≤ e ≤ m + 1, where n is the number of rows and m is the number of columns. As seen in the initialization of the recurrence, the interval of the first row starts as empty; the cost of this interval is the sum of costs of all errors: all 1s on that row. We next describe the three recurrence operators r 1 , r 2 , r 3 used in the recurrence.

C(1, 1, 1 ) ← j=1...m W1 j 1 C(r, s, e ) ← min{r 1 , r 2 , r 3 }, where
-The value of r 1 is available if and only if s < e. It represents moving the end of the interval one position ahead on row r, essentially adding one entry to the interval. This either introduces an error with cost W0 e-1 r , or eliminates an error with cost W1 e-1 r .

-The value of r 2 is available if and only if s > 1. It represents moving the start of the interval one position ahead on row r, thus removing the leftmost entry from the interval. This either introduces an error with cost W1 s-1 r , or eliminates an error with cost W0 s-1 r .

-Value r 3 is available if and only if r > 1. It represents fixing the interval on previous row r -1 and proceeding to the next row r of M . Initially, the interval for row r is identical to that of row r -1, updating the cumulative cost as follows: all 1s outside the current interval are errors, as well as all 0s inside the interval. Each error contributes to the total cumulative cost.

Once the recurrence has completed, the optimal solution is the minimum cumulative cost among the last-row entries of C: min 1≤s≤e≤m+1 {C(n, s, e )}. An exact sequence of optimal intervals can be traced by storing the operator used on each entry.

The time complexity for the Optimal Visual Band is O(nm 2 ), caused by the size of array C. Each sum-operation in the recurrence can be evaluated in constant time: instead of using the costs in W1 and W0 directly, we construct cumulative weight matrices of both W1 and W0 on rows beforehand and evaluate the sums using these matrices. Memory consumption can be lowered to O(m 2 ), if cumulative costs for row r -1 are replaced by costs for row r as they are computed; the time consumption remains unchanged.

Algorithms for Mbs

Next we turn to the banded submatrix problem Mbs. Consider a solution matrix of a band-algorithm. As is common in information retrieval, we can compute, for each row and column the number of true positive (TP) entries as the number of 1s in the row (or column) that were not transformed into a 0. That is, they correspond to relevant 1s contained in the row (column) for bandedness. Similarly, the false positive (FP) entries of a row (or column) correspond to original 0s that were transformed into 1s by the algorithm. Finally, false negatives (FN) are represented by those original 1s that were switched into 0 by the banded algorithm.

From here many different measures for evaluating the performance of the banded retrieval system can be computed. For example we can consider precision =TP/(FP+TP), recall=TP/(FP+FN) or accuracy=(TP+ TN)/n of each one of the row and column vectors. The algorithm proposed is called Iterative row-column elimination. As the name indicates, it will eliminate rows/columns one at a time, by selecting the one with worst performance so far. Pseudocode of the algorithm is shown in Algorithm 5.

Experiments

In this section we give experimental ground of the bandedness concept. We demonstrate that our algorithms are computationally effective and produce good results in synthetic datasets. On top of that, we show that the algorithms are able to find banded patterns in a variety of real datasets from life sciences. The banded structures in Figure 2 were found by Alternating algorithm. Rank columns/rows with a performance measure 6:

Select the column/row i with worst performance 7:

Remove i from M 8: until ǫ ≤ k

Synthetic Data

Data generation Briefly, we generated synthetic data as follows. For given n and m and a width parameter 2w, we generate a fully banded matrix by means of a random walk in the matrix grid. Starting at coordinate (0, 0) of an n × m matrix initially set to all zeros, the random walk chooses to move either one step down (i.e. from (i, j) to (i + 1, j)) or one step to the right (i.e. from (i, j) to (i, j +1)) with equal probabilities. Whenever a step to the right is chosen, we will set to 1 all the w entries above the current position and all the w entries below the current position (or less than w, if i < w). The random walk is designed so that it always reaches the final position (n, m). When it does, the matrix has a clear solid band of maximum thickness 2w like shown in Figure 1.

Additionally we can introduce noise by flipping the original values, 0 to 1 or 1 to 0, according to given probabilities. We use Pr(0 → 1) as the noise probability for 0-to-1 flips and Pr(1 → 0) for 1-to-0 flips.

We generated samples of 50 × 55 binary matrices, which are used in all synthetic data experiments. The width parameter was 2w = 30, which yields approximately 50% fill of 1s. Before running the experiments on these matrices, noise was added and row-column permutations were randomized. We prepared a variety of test settings, including both balanced noise (when Pr(1 → 0) is the same as Pr(0 → 1)) and unbalanced noise (when they are different), and different weighting schemes.

In all the synthetic experiments we generated 30 sample matrices for each one of our parameter settings. The values in the figures of this section are the averages of the values obtained from the samples. We will next summarize the performance exhibited by our algorithms on synthetic data.

Methods The algorithms used in experiments include Alternating, Barycentric, and Simulated Annealing, as well as Fixed Permutation and Bidirectional Fixed Permutation coupled with several similarity measures and column-ordering methods. Eventually, we used Optimal Visual Band method to evaluate minimumflips values βB (M ) for each algorithm after row-column permutations were fixed by each one of the methods.

In order to make the evaluation of algorithms more reliable, we introduce two new methods, Original and Random. In Random method, row-column permutations are selected randomly. Another competitor is Original method, which has access to the original "correct" permutations. By correct we mean that the original fully banded structure was generated under these permutations. Once noise has been added, the matrix is no longer fully banded, but Original gives good permutations nonetheless.

As expected, the Bidirectional Fixed Permutation strategy outperformed the Fixed Permutation in all cases where noise exists, so we decided to leave the Fixed Permutation algorithm out of the results presented next. For the Bidirectional Fixed Permutation algorithm, we sorted columns via the spectral ordering with three similarity measures: Pearson correlation, dot product, and Jaccard similarity. We also experimented with finding good columnsorderings via the Hamiltonian path approximation (specifically, minimum spanning tree approximation) with both Hamming distance and Jaccard distance. For visual clarity, only the best results are shown: spectral ordering using dot product and Hamiltonian ordering with Jaccard distance.

For the Simulated Annealing method, we first compare the different neighbor schemes: Swap-k and Adj-swap-k for k ∈ {1, 2, 4}; Reverse, Relocate, and Reverse+Relocate. Results are shown in Figure 5(a). The number of iterations was one million and the starting temperature was 10. These parameters produced results that we consider near-optimal. In terms of convergence speed, Swap-1 proved to be the fastest, as seen in Figure 5(a). This scheme was chosen to be the neighbor scheme for all experiments that include Simulated Annealing.

Because of the 2-dimensional nature of band-problems, the number of iterations that Simulated Annealing needs is larger than other ordinary permutation problems of same size. We studied the convergence of Simulated Annealing with Swap-1 scheme extensively by using different number of iterations and temperature-multipliers. Using a large multiplier requires significantly more iterations, which makes the largest multiplier values impractical for this experimental setting. We chose a running-time limit of 2 minutes and settled for multiplier α = 0.9999 with 100000 iterations. Other parameter choices that satisfy the time limit were tested but they failed to produce superior results consistently. The chosen parameter combination is conservative: temperature is allowed to drop low enough so that further improvement is unlikely, thus reducing variance. These parameters yield solutions that have, on average, 3% higher cost than those from one million iterations.

In most cases, both Alternating algorithm and Barycentric algorithm converge very fast, taking less than 30 iterations in average to find good permutations, regardless of matrix size and initial permutations. After the first iterations, the Alternating algorithm does not strictly converge towards an optimal solution, but starts oscillating between good and very good solutions, as seen in Figure 5(b). Because of this, we return the best solution among the first 100 iterations. Results A summary of the main results follows; first, we use the minimum-flips as quality measure and then a notion of rank correlation measure. We added balanced noise to the sample of 30 synthetically generated fully banded matrices and then run each of the presented algorithms to obtain the βB (M ) values for the matrices. The average results from the 30 runs were then computed for each method, as shown in Figure 6(a). We can see that both Alternating algorithm and Simulated Annealing algorithm consistently perform as well as the Original method, and at higher noise levels, they perform even better. The reason is that the original banded structure starts to disappear with more noise-the original permutation is no longer the best one with respect to minimum-flips value. The other methods do not beat Original until extreme noise levels. We discovered that the Hamiltonian path method (in a fixed permutation algorithm) has large variance, caused by non-optimal cycle-to-path transformations.

A A A A A A A A A A A A A A A A A A A A A A A A A A A R R R R R R R R R R L L L L L L L L L L B B B B B B B B B B A S A B A L S R
In the unbalanced noise case we chose to fix one of the noise probabilities Pr(0 → 1) = 0.1 and added different noise levels of Pr(1 → 0) to the input matrices; again, with a sample of 30 matrices in total. To make the comparison reliable, we use an equal cost weighting scheme. As in the balanced noise case, the average result of 30 runs was collected for each method; the results are shown in Figure 6(b). Again, the Alternating algorithm and the Simulated Annealing algorithm perform very well, beating the Original method at noise levels Pr(1 → 0) > 0.25. Overall, Alternating seems to produce best average results, although we noticed its variance to be a bit larger than that of the Simulated Annealing algorithm. Curiously, the Barycentric method performs poorly here: at high noise levels, the results hardly beat a random permutation. Spectral methods perform a bit better, whereas the results for Hamiltonian methods are two-fold: mediocre results at low noise levels; good results at high noise levels.

Finally, we would like to know how well our methods can recover the original permutations after we have added noise and randomly reordered the rows and columns. For this we use Spearman rank correlation, which compares two sets of rankings and computes a correlation coefficient that describes the similarity of the rankings. Here we understand a ranking as an ordering of rows and columns. We compare two rankings, namely original and recovered: the original ranking matches the row order in the generation process of noiseless data; the recovered ranking comes from the new row ordering retrieved by any of our algorithms. Correlation 1.0 means perfect agreement between rankings; values close to 0.0 indicate that there is no correlation between the two.

The rank correlation results for balanced noise are shown in Figure 7(a). We see that none of the methods dominates others: at low noise levels under 0.15, Spectral methods are able to recover the original order almost perfectly; from 0.15 to 0.30, the Alternating algorithm seems to be the best choice; at levels of noise over 0.30, the Barycentric method is the best. Overall, these three methods are very good in recovering the original row order. Hamiltonian methods, however, are able to find the original order only at lowest noise levels.

Rank correlation with ties

Pr(0<->1) rank correlation mean

Occasionally the Alternating method converges towards local minimum, as demonstrated by correlation value below 1.0 with noiseless data.

Figure 7(b) shows the rank correlation results for unbalanced noise. To deal with unbalanced noise we analyze the input data and pass a weight parameter to the methods (except for Barycentric). We use a simple but effective weighting scheme: compute the relative proportion of 1s and 0s in input matrix, namely p 1 and p 0 , where p 1 + p 0 = 1. If the weight for a 0-to-1 flip is 1, then the weight of a 1-to-0 flip is log(p 1 )/ log(p 0 ).

Unlike other methods, Barycentric is unaffected by the nature of noise (balanced or unbalanced), and does not need a weight parameter. Once the appropriate weight parameter is analyzed from input data, Alternating and Spectral methods perform very well. Until noise level 0.6 these methods are the best at recovering original band-order. As it was the case with balanced noise, Hamiltonian methods have larger variance and difficulties at recovering the original order.

The following table shows the running times (ms) for the methods on a 3.1 GHz Intel machine. All synthetic matrices have 50% fill and 10% balanced noise. The Barycentric and Alternating are the fastest methods, even though both perform 100 iterations. Their running times increase steadily along the matrix size, while the increase is more pronounced for Hamiltonian and Spectral methods. Overall, all these methods are fast and able to handle even larger matrices in reasonable time.

In the next sections we will present a series of real data applications. We will use the algorithms that are most appropriate for the data. Given the analysis in synthetic data these algorithms will be either Alternating, Barycentric or Bidirectional Fixed Permutations that uses the Spectral method and dot product for fixing the column permutation. Notice that even if Simulated Annealing could give good results, it is too time consuming for practical real data.

Mammals Data

The Mammal dataset consists of presence/absence records of European mammals: for 2179 cells of size 40 × 40 km 2 we have binary information about 124 different species. The full version of the Mammal dataset is available for research purposes upon request from the Societas Europaea Mammacologica (www.european-mammals.org).

Figure 8(a) shows the band obtained by the Alternating method that allows bidirectional flips with equal costs. The discovered structure of sites and mammals is interesting in itself: it shows a highly nested pattern from sites and species so far unknown in the dataset. Nested patterns have been widely studied in ecological applications because they give insight into the processes that govern the distribution of species. The correlation between the order of the sites and the temperature variable is 0.64, which is high for life-science data. Intuition is that the number of mammal species increases to some extent with temperature and few species live in cold conditions. Furthermore, of the 54155 total 1s in the dataset, 70% of them are accumulated within the borders of the band retrieved by Alternating algorithm (see Figure 8(b)). This means we found a dense structure. The matrix needs 23997 bidirectional flips to become a fully banded matrix.

As a comparison, Figure 8(c) depicts the band obtained by the Barycentric method. The bidirectional flips needed to make this matrix fully banded is 42041 flips. Nonetheless, the row ordering has high correlation with certain variables associated with the location, such as latitude coordinate (0.92) and average temperature (-0.85). It seems that both Alternating and Barycentric produce reasonable, yet visually different results here: Alternating prefers large dense clusters of 1s, whereas Barycentric assumes that all rows and columns belong to the banded pattern. Which shape is preferable depends on application. Note that assigning different weights for flips causes the Alternating method to change the shape of the band.

Dialect Data

The Dialect dataset [START_REF] Embleton | Computerized dialect atlas of Finnish: Dealing with ambiguity[END_REF], originally published in 1940, contains data about the usage of dialectical features in spoken Finnish language. The data is in binary form and represents 1334 phonological features and their usage in 506 municipalities.

The basic division of Finnish dialects has long remained static among linguistics. The division into two dialects results in Western and Eastern dialects; further divisions bring out more detail inside these two main dialects. The known division into eight dialect areas is shown in Figure 9(a).

We used both the Alternating algorithm and the Barycentric algorithm on the Dialect data to uncover a band structure. The best ordering found by the Alternating method after 1000 iterations is shown in Figure 10 the band of the Barycentric method. The two bands are very different: while Alternating aims at creating dense bands of 1s, the Barycentric method aims at creating wider banded patterns.

Interestingly, the order obtained for municipalities in the band from the Alternating algorithm can be interpreted as a crude estimate of the language spoken in a municipality: the order captures the variation from Western dialects to Eastern dialects. Indeed, the municipalities interpolated from the ordering of the band are plotted in Figure 9(b). We can see that this interpolation is very similar to the known division shown aside in Figure 9(a), despite some outliers. The outliers are explained by the small number of data points (1s) for those municipalities. Of course the dialectical variation cannot be fully described by a single ordering: some municipalities do not really fit anywhere in the linear order, for example those in Karelia. Note that the results show many geographically coherent dialect areas, although the algorithm did not have access to spatial information about the municipalities.

Paleontological Data

The Paleontological data contains information of fossils genera in Europe [START_REF] Fortelius | Neogene of the old world database of fossil mammals (NOW)[END_REF] 1 . Columns correspond to genera and rows correspond to sites. There are 124 sites and 139 genera. The banded structure found in the Paleontological dataset is shown in Figure 11. The result is given by the Alternating method.

As is common with data gathering in real life, the paleontological data has more missing 1s than false 1s. Therefore 0-to-1 flips should be preferred over 1-to-0 flips. In general it is difficult to decide for a real life application which are the best weights, and many criteria for "correct" band-shape can be given. On this data, we decided to give to 1-to-0 flips a weight that is four times larger that that of a 0-to-1 flip, as it matches with our estimation of error rates.

DNA Amplification Data

The DNA amplification data is available upon request from the authors in [START_REF] Myllykangas | DNA copy number amplification profiling of human neoplasms[END_REF]. It contains information on the DNA copy number amplifications recorded in 4590 cases (rows) and 393 band specific chromosomal locations (columns). A 0 denotes no DNA copy number amplification in the corresponding chromosomal location, and a 1 denotes a finding in the DNA copy number amplification in the location. More than this we also have available the neoplasm labels (cancer label) coupled to the cases of the matrix. The goal is to investigate DNA amplifications in different neoplasms types. This justifies a solution where different subsets of columns and rows can be evaluated separately, so we will run our Iterative row-column elimination algorithm. We use the performance measures of accuracy and recall as a deletion criteria of rows and columns, while allowing a maximum of 100 flips in the retrieved submatrices.

Figures 12(a) and 12(b) show the two retrieved submatrices (from the original data) found by accuracy and recall. For both submatrices we observe immediately that they exhibit a very clean banded structure (almost zero noise outside the band). They both preserve nicely the column permutation of the chromosomes, and moreover, each one of these submatrices identifies a different reduced collection of cancer types.

Related Work

The property of consecutive ones in binary matrices is particularly important in presence/absence data of paleontological applications [START_REF] Puolamäki | Seriation in paleontological data using Markov chain Monte Carlo methods[END_REF]. The question of determining how far away a given matrix is from satisfying the consecutiveones property corresponds to counting Lazarus events. This is also known as the seriation problem, to which Spectral ordering has been shown to produce good approximations [START_REF] Atkins | A spectral algorithm for seriation and the consecutive ones problem[END_REF].

Another important ecological concept related to the structure of binary matrices is nestedness [START_REF] Mannila | Nestedness and segmented nestedness[END_REF]. A dataset is nested if for all pairs of rows one row is either a superset or subset of the other. From this definition it is direct to see that all nested matrices are banded, although the reverse does not hold. Bandedness is therefore a generalization of these ecological concepts.

A fully banded binary matrix is also known in some contexts as the monotone consecutive arrangement of ones [START_REF] Sen | Indifference digraphs: A generalization of indifference graphs and semiorders[END_REF]. Then, we can establish a connection between banded matrices and zero-partitionable matrices from [START_REF] Lin | Interval digraphs that are indifference digraphs[END_REF]. We say that Two submatrices from the DNA data. On the left (a), the submatrix with 1402 rows and 282 columns, using accuracy as the elimination criteria. On the right (b), the submatrix of 291 rows and 244 columns, using recall as the elimination criteria. Permutations on the columns were fixed by the spectral ordering with dot product similarity. The number of flips to obtain a full band is at most 100. a binary matrix is zero-partitionable if the rows and columns can be independently permuted such that each 0 can be labeled with R or C in such a way that every position to the right of an R is an R and every position below a C is a C. Therefore, any matrix satisfying the consecutive ones property will be always zero-partitionable, we only need to list rows in nondecreasing order of the position of their left most ones. The inverse statement is not true. This establishes the following hierarchy between classes of matrices,

N ⊂ B ⊂ SC1P ⊂ C1P ⊂ Z
where N are nested matrices, B are banded matrices, SC1P are simultaneously consecutive-ones on row and columns, C1P are consecutive-ones matrices on rows, and Z are zero-partitionable matrices.

We cannot ignore the relation of banded matrices with the field of numerical analysis [START_REF] Aykanat | Permuting sparse rectangular matrices into block-diagonal form[END_REF][START_REF] Cuthill | Reducing the bandwidth of sparse symmetric matrices[END_REF][START_REF] Rosen | Matrix bandwidth minimization[END_REF]. The solutions there also find good permutations to confine non-zero entries to the main diagonal band. Noisy data is, however, a problem for them, as the algorithms require that all 1s must be close to the main diagonal. The approach in this paper allows banded structures outside the main diagonal.

In data mining, bands are closely related to reorderable matrices, see e.g. [START_REF] Bertin | Graphics and graphic information processing[END_REF][START_REF] Mäkinen | Reordering the reorderable matrix as an algorithmic problem[END_REF]: this is a simple visualization method to explore multivariate or multidimensional data; the basic principle is to transform a multidimensional data set into a 2D interactive graphic. The graphical presentation of a data set closely resembles the underlying matrix in that it contains rows and columns which can be permuted, allowing different views of the data set. The problem of reorderable matrices can be applied to real-valued data and to binary data. The property of bandedness can be seen as a subcase of reorderable matrices where we only have binary values and the visualization must satisfy the band constraints. Therefore, our contribution offers new algorithms and theoretical properties for a special case of reorderable matrices.

Finding a band structure is closely related to biclustering problems [START_REF] Madeira | Biclustering algorithms for biological data analysis: A survey[END_REF] and subspace clustering [START_REF] Moise | Subspace and projected clustering: experimental evaluation and analysis[END_REF][START_REF] Parsons | Subspace clustering for high dimensional data: a review[END_REF]: the goal is to identify groups of columns and rows that exhibit similar value patterns (for example, similar expression patterns in the case of microarray data). Many combinatorial problems have been studied around biclustering and the associated bipartite graph. The most related one to the one we study here is the edge editing problem, consisting in adding or removing the fewest edges from a bipartite graph so that it becomes a vertex disjoint union of complete bipartite graphs. However, the definitions are slightly different, as in our case the different biclusters (groups of columns and rows with ones) do not require to be vertex disjoint.

Finally, our approach can be see as a cluster ranking approach, as in [START_REF] Bar-Yossef | Cluster ranking with an application to mining mailbox networks[END_REF]. Rather than simply partitioning a network into clusters, a cluster ranking algorithm also orders the clusters by their strength. In our problem, the different "clusters" or groups of rows are mapped into a one-dimensional rank provided by the columns of the same matrix. An important difference with those approaches is that we do not depend on parameters such as the number of clusters, which is a well studied problem in clustering categorical data [START_REF] Chen | Best k": critical clustering structures in categorical datasets[END_REF].

Conclusions

We introduce the new concept of banded structures for binary data and characterize exactly its combinatorial properties. We first establish how row and column permutations affect the structure of the band and derive a polynomial test for bandedness. Based on the derived properties, we propose principled algorithms to approach two important combinatorial problems related to bands: minimum banded augmentation and maximum banded submatrix.

A first set of algorithms relies on fixing a column permutation beforehand. We then lift this column permutation requirement and propose two new algorithms that search for both column and row permutations at the same time. We show in synthetically generated data that our proposals are able to find banded structures efficiently.

Experiments on real data show that bands potentially occur in many applications, and here we concentrate on datasets from life sciences and a word dialect data. For two of the datasets, we discover previously unknown bands that have natural interpretations in the final order of rows and columns. Our results suggest that a hierarchy of mammals exists in Mammals occurrence data. For the word dialect data we have discovered that the band order has captured the main variation between Eastern and Western dialects of the Finnish language. In bandedness real life binary data is given a combinatorial treatment so as to uncover hidden structure in the data.

Fig. 1 .

 1 Fig. 1. An example of a fully banded matrix.

Fig. 3 .

 3 Fig. 3. Band exhibited in two network datasets after proper row and column permutation: Football (a) and Les Misérables (b).

Fig. 4 .

 4 Fig. 4. The incompatibility graph corresponding to the matrix from Equation 1, with row sets M1 = {a, b, c}, M2 = {b, c, d}, and M3 = {c}. The upper figure shows the incompatibilities introduced by only the first row M1 = {a, b, c}, together with the rule (x, y) and (y, x) are incompatible. The lower figure shows the whole graph obtained by adding the incompatibilities from all rows. M has the consecutive-ones property as this graph is bipartite.

  y] with the best from the following options: 10: (a) x = min{aj | [aj, bj] ∈ C} and y = b 11: (b) x = a and y = max{bj | [aj , bj] ∈ C} 12: (c) x = aj and y = max{b k | [a k , b k ] ∈ C, a k < aj}, for every M π j = [aj , bj] ∈ C 13: end for 14: Sort the rows [a, b] of M π in ascending order of as, resolving ties with the ascending order of their bs.

7 :

 7 [a, b] be the maximum consecutive subarray on W π i Update M π i = [a, b] 8: end for 9: for each pair of rows i, j in M π do 10:

Algorithm 4 1 :

 41 The Simulated Annealing algorithm for Mba or Bidirectional Mba Input: An n × m binary matrix M , a number of iterations t, a large temperature T and a temperature multiplier α ∈ [0.95, 1) 2: Output: A permutation κ of rows and a permutation π of columns 3: Initialize state S = M and energy E = Optimal_Visual_Band(S) 4: Initialize best solution: S best = S and E best = E 5: repeat 6:

r 1 = 1 r-W1 e- 1 r r 2 = 1 r+ W1 s- 1 r r 3 =

 1112113 C(r, s, e -1 ) + W0 e-C(r, s -1, e ) -W0 s-C(r -1, s, e ) +

Algorithm 5 The

 5 Iterative row-column elimination algorithm for Mbs 1: Input: An n × m binary matrix M and an integer k 2: Output: A submatrix M ′ of M 3: repeat 4: Compute ǫ, the number of flips estimated for Mba or Bidirectional Mba problem on M 5:

Fig. 5 .

 5 Fig. 5. Left part (a) shows the convergence analysis and comparison of neighbor schemes for Simulated Annealing algorithm. The fastest converging schemes are swap-1 and reverse. Right part (b) shows the convergence analysis for the Alternating algorithm. On the y-axis: minimum number of flips necessary to make the matrix fully banded; on the x-axis: number of iterations used.

Fig. 6 .

 6 Fig.6. Left part (a) shows the average minimum-flips comparison of the different methods on synthetic data when balanced noise has been added to the original matrices; the x-axis corresponds to the variation of noise for both Pr(0 → 1) and Pr(1 → 0). Right part (b) shows the average minimum-flips comparison of methods on synthetic data with unbalanced noise; the x-axis corresponds to the variation of noise for Pr(1 → 0) (where Pr(0 → 1) = 0.1 is fixed all the time). The average value of 30 samples is shown for both plots. Alternating and Simulated Annealing are the best performers.

Fig. 7 .

 7 Fig.7. Left part (a) shows the average rank correlation comparison of all the methods on synthetic data with balanced noise. On the y-axis: Spearman rank correlation measure; on the x-axis: variation of noise for both Pr(0 → 1) and Pr(1 → 0) (balanced noise scenario). Spectral methods, Alternating, and Barycentric all produce good results. Right part (b) shows the average rank correlation comparison of all the methods on synthetic data with unbalanced noise. On the y-axis: Spearman rank correlation measure; on the x-axis: variation of noise for Pr(1 → 0) (unbalanced noise scenario, where Pr(0 → 1) = 0.1 is fixed). Basic Barycentric method performs well by itself; Spectral methods and Alternating are equally good, when the weight parameter is employed.

Fig. 8 .

 8 Fig.8. Reordered Mammals data. Left part (a) shows the data reordered by the Alternating algorithm, using equal-cost as a weighting scheme. Middle part (b) shows the borders of the band, that is all the 1s in original Mammals data that belong to the structure. In total the band includes 70% of all the 1s. Right part (c) shows the data as reordered by the Barycentric algorithm. The order of locations has strong correlations with latitude coordinates (0.92) and average temperatures (-0.85) of the locations.

Fig. 9 .

 9 Fig. 9. Left part (a) shows the traditional division of Finnish municipalities into eight areas by their dialectical features. Right part (b) visualizes the ordering of municipalities interpolated from the band discovered by the Alternating algorithm. The position of a municipality in band-order is depicted by its color, ranging from black to white. The ordering of the band has captured the main variation between Western and Eastern dialects.

Fig. 10 .Fig. 11 .

 1011 Fig. 10. Left part (a) shows the best band found by the Alternating algorithm on the Dialect data. This band is 54461 flips away from being fully banded. The right part (b) shows the band found by the Barycentric algorithm. This band is 105410 flips away from being fully banded.

  Fig.[START_REF] Cormen | Introduction to algorithms[END_REF]. Two submatrices from the DNA data. On the left (a), the submatrix with 1402 rows and 282 columns, using accuracy as the elimination criteria. On the right (b), the submatrix of 291 rows and 244 columns, using recall as the elimination criteria. Permutations on the columns were fixed by the spectral ordering with dot product similarity. The number of flips to obtain a full band is at most 100.

requirements; (3) compares to previous related work; and (4) includes an extended experimental evaluation with a view to new applications.

FromAugust 31,

2007, the version of the dataset selects all columns with at least 10 occurrences and all rows with at least 10 genera present.


Yaset Oliva
email: yaset.oliva@insa-rennes.fr

Maxime Pelcat
email: maxime.pelcat@insa-rennes.fr

Jean-Francois Nezan
email: jean-francois.nezan@insa-rennes.fr

Jean-Christophe Prevotet
email: jean-christophe.prevotet@insa-rennes.fr

Building a RTOS for MPSoC Dataflow Programming

Multiprocessor Systems-on-Chip (MPSoC) are becoming the standard high performance Digital Signal Processing (DSP) systems. Hardware complexity abstraction is needed to enable efficient MPSoC programming. A major challenge of MPSoC programming is efficiently handling the combination of new features necessary in a MPSoC operating system: load balancing and efficient use of the parallel resources, with the more traditional features of Real-Time Operating Systems (RTOS): resource sharing between applications, task priorities and reactivity to events. This paper presents a method to combine dataflow methods and RTOS features. The resulting system prototypes an RTOS for symmetric multiprocessing MPSoCs whose inputs are dataflow graphs of applications. The prototype is built on the µC/OS-II RTOS. Experimental results are given on a 3GPP Long Term Evolution algorithm executed on a 4-core MPSoC.

I. INTRODUCTION

In [START_REF] Lee | The problem with threads[END_REF], Edward Lee shows that programming with threads is an error-prone operation and proposes several alternatives, including using process networks, to make software behavior more predictable. The dataflow process network [START_REF] Lee | Dataflow process networks[END_REF] Model of Computation (MoC) models an algorithm by concurrent and independent modules known as actors which communicate ordered tokens (data quanta) through First-In First-Out channels. A set of firing rules defines when an actor executes. Dataflow models have been shown to favor parallel algorithm description as they favor data locality and reduce multi-core scheduling constraints to data dependencies [START_REF] Ackerman | Data flow languages[END_REF]. Thus, dataflow models are well suited for use with signal processing algorithms MPSoC systems used for signal processing are increasingly complex to program. Tools that ease MPSoC programming are thus more and more needed. This paper outlines a method of using dataflow graphs instead of thread declarations as inputs to an RTOS. In this case, the RTOS is then able to dispatch the actors which compose the dataflow graph to the cores of an MPSoC. These actors are directly executed by the MPSoC RTOS, instead of requiring the programmer to program primitives for task migration and synchronization. This method ensures the synchronization of actors and the management of actor input and output data.

This paper presents experiments performed with a MPSoC RTOS prototype incorporating dataflow model management into the µC/OS-II kernel. The first target application is the uplink data decoding algorithm executed in base stations supporting the 3GPP Long Term Evolution (LTE) telecommunication standard.

Section II presents related works on multi-core runtime management and dataflow MoCs. Section III introduces the MPSoC RTOS structure. Section IV explains the dataflow management part and Section V considers the multi-core scheduling part. Experimental results of the prototype are shown in Section VI.

II. RELATED WORKS

A. RTOS and Runtime Management Systems

In [START_REF] Nollet | A safari through the MPSoC Run-Time management jungle[END_REF], Nollet, et al. present an overview of runtime management systems for MPSoCs. The overview covers both industrial and academic systems. Commonalities between systems are identified: their structure can be divided into two parts, the quality manager and the resource manager. The quality manager tries to optimize the Quality of Service (QoS) of the system, i.e. to find the best application configuration; the resource manager offers mechanisms to allocate cores, communication media and memory.

The runtime systems reviewed in [START_REF] Nollet | A safari through the MPSoC Run-Time management jungle[END_REF] consider only sequential applications which share cores of an MPSoC. The MPSoC RTOS presented in this paper differs, as it partitions each application between available cores. Parallelism of each application is explicitly stated using a dataflow graph. Such partitioning aims at ensuring good load balancing between cores and at reducing the computation latency of each application. In this way, the MPSoC RTOS may be considered to be equivalent to a quality manager of low granularity. The StreamIt [START_REF] Thies | StreamIt: a compiler for streaming applications[END_REF] runtime system has similar goals than the MPSoC RTOS. However, it processes a specific streaming language while the MPSoC RTOS reuses C/C++ legacy code for the actors and combines it with a parameterized dataflow coordination language.

The MPSoC RTOS prototype is built on an existing kernel: µC/OS-II [START_REF] Labrosse | MicroC/OS-II: the real-time kernel[END_REF]. This choice was made through consideration of the small footprint, the simplicity and the available source code for this kernel. Since the original µC/OS-II doesn't implement any multiprocessor mechanism, its sources are modified to manage MPSoCs.

B. Modeling Applications with a Parameterized Dataflow Model

Many dataflow MoCs have been introduced in the literature, each offering a tradeoff between compile-time predictability and capacity to model dynamic run-time variability. It may be seen that the most obvious difference is their firing rules. Signal processing applications, such as telecommunication or video processing, are based on a loop which repeats a pattern of execution on a sequence of input data. Describing the repeated pattern with a dataflow MoC consists of dividing the computation into actors that exchange data only through input and output data queues without sharing any state.

In this paper, applications are described using parameterized dataflow modeling [START_REF] Bhattacharya | Consistency analysis of reconfigurable dataflow specifications[END_REF]. The chosen model of Parameterized Cyclo-Static Directed Acyclic Graph (PCSDAG) states that the graph can be totally reconfigured once before starting a new execution. This model is a parameterized and acyclic version of the Cyclo-Static Dataflow (CSDF) model [START_REF] Bilsen | Cyclostatic data flow[END_REF]. The reconfiguration enables dynamic behavior of the algorithm, i.e. computation strongly depending on the input data. The PCSDAG model has been shown to be suitable for describing uplink and downlink data processing algorithms of 3GPP LTE base stations [START_REF] Pelcat | Adaptive multicore scheduling for the LTE uplink[END_REF].

III. STRUCTURE OF THE MPSOC MODEL-BASED RTOS

The RTOS obtained is divided into two modules: dataflow management and RTOS scheduling (Figure 1). The dataflow graph is described with C++ objects and compiled within the dataflow management module. At each OS clock tick, the dataflow graph is parameterized, i.e. new parameter values are retrieved, and a temporary graph of execution is obtained. For each actor within the temporary graph, an RTOS task is created. The task corresponding to an actor must manage actor data and call actor code. Flags are automatically generated to synchronize actors based on their dependencies. RTOS scheduling is completed by assigning each actor to a core.

IV. DATAFLOW GRAPH MANAGEMENT

The test application used is the uplink decoding of 3GPP LTE base stations. The PCSDAG model was specifically designed for this application but is likely to be also suitable for video encoding and decoding algorithms. Managing the graph consists of performing an expansion. This is the process where the dataflow graph is transformed into the temporary graph in which each actor is executed only once in a graph iteration. The temporary graph depends strongly on parameter values that fluctuate between consecutive iterations.

Figure 2 illustrates the application and the expansion. Uplink decoding consists of retrieving data sent simultaneously by nb user users. This data is divided into Code Blocks (CB) of variable size. Both the number of users nb user and the number of code blocks nb CB for each user vary every millisecond. nb CB is different for each user and may be represented by a cyclic pattern nb CB = nb CBuser0 , nb CBuser1 , .... The computation is divided into four phases: Multiple Input Multiple Output (MIMO) decoding, symbol to bit conversion, bit processing, and Cyclic Redundancy Check (CRC) [START_REF] Pelcat | Adaptive multicore scheduling for the LTE uplink[END_REF]. Each millisecond, the graph is parameterized and transformed into a single rate directed acyclic temporary graph (srDAG) in which production and consumption rates of each FIFO queue are equal. Figure 2 shows 4 possible configurations. The expansion phase has been demonstrated to be executable in real-time while respecting the 3GPP LTE constraints in [START_REF] Pelcat | Adaptive multicore scheduling for the LTE uplink[END_REF]. After the expansion phase, an RTOS task is created for each actor in the srDAG graph. These tasks are then scheduled on the multiple cores of the MPSoC architecture. It may be noted that an alternate algorithm model to PCSDAG can be chosen provided that the srDAG graph can still be produced.

V. MULTI-CORE SCHEDULING

RTOS scheduling is divided into two phases: the master/slave phase and the symmetric phase. In the master/slave phase, the core designated as master executes the dataflow graph management and posts actors to the slave cores and to itself. When the scheduling enters the symmetric phase , each core is part of a pair and calls successively the common µC/OS-II scheduler to execute the highest priority task. The objective of this division is to limit the master/slave phase that naturally puts the master in the position of the bottleneck. The symmetric phase also enables preemptions and passive wait of events.

A. Master/Slave Phase

The master/slave phase is illustrated in Figure 3. Each core has a 1-place queue to receive an order of task execution. After an OS timer tick, the previously described dataflow graph management is performed and one task per actor is created. The master processor then dispatches the ready tasks with highest priority to the available slave processors. The dispatching process is performed by adding a message into the 1-place queue of the selected slave. The message contains the address of the shared memory location where the code of the task has been placed as well as input/output buffer management information. The dispatching process finishes when all ready tasks have been mapped or when there are no more available slave processors. When all slave processors are busy and there are still ready tasks, then the master processor executes the next task itself.

At the same time as the master executing, each slave processor waits for a message to be placed into its 1-place queue. When the message arrives, the slave processor places its program counter to the designed address and executes the task. After this first task dispatching is done, the symmetric phase starts and all cores become pairs.

B. Symmetric Phase

During the symmetric phase (Figure 4), each core can access the schedule function. This common access requires the use of mutex semaphores provided within the architecture support library. As in a typical RTOS, the scheduler function consists of identifying the highest-priority task that is ready and running it. The scheduler is called if either the task execution is preempted by a higher execution task or if the task execution is completed.

After the execution of the schedule function, and if no task is ready to be executed by the core, the scheduler saves the current task's context state and gets returned to its private memory. The slave is now ready to receive another message from the master processor.

Using RTOS scheduling, several independent applications can share the MPSoC with each actor with its own priority. Moreover, passive wait of events is allowed. This is very important in the case of systems with co-processors. In the cases where a programmer wishes to offload a costly operation such as turbo-coding on a co-processor, an actor can wait for turbo-coding completion passively on a core and can be preempted by another ready actor while turbo-coding is running.

The next section gives experimental results of the prototype.

VI. EXPERIMENTAL RESULTS

The experimental prototype is executed on an Altera Cyclone II FPGA integrating 4 NIOS cores [START_REF] Altera | Literature on NIOSII software embedded processor[END_REF] and a shared memory. The memory footprint of the system is shown in Figure 5. As can be seen, the operating system data occupies called from all cores Fig. 4. Symmetric Scheduling Phase: Self-Organizing Execution more than 50% of the 510kb of total system memory. This expensive cost is largely due to the size and number of task stacks statically allocated by the kernel (64 in this case). The prototype has been run for six iterations. At the end of each iteration, the application parameter values are modified to generate a different number of actors, as shown in Table I 6 shows the system performance in terms of execution time, against the number of cores. The iterations correspond to the ones presented in Table I. These curves confirm that the system performances improve as the number of cores increases. However, the concurrent shared memory accesses also increase and become a bottleneck. This is a limitation of the current system. The performance improvement naturally depends on the shape of the srDAG because the exposed parallelism depends on the graph parameters.

VII. FUTURE WORKS

The current prototype may be improved both in terms of latency and in terms of memory. Speed can be gained by exploiting quiescent points [START_REF] Neuendorffer | Hierarchical reconfiguration of dataflow models[END_REF] of the single rate DAG. Quiescent points are points between two actors' execution when no context needs to be maintained during rescheduling. This case corresponds to the scheduler call at the end of the "Task execution code" in Figure 4. This call does not require context switch and the task stack can even be flushed.

From a memory viewpoint, the task stacks can be significantly reduced. µC/OS-II associates independent stacks to each task. Having a pool of tasks and reusing tasks associated to finished actors will automatically reduce task creation time and stack memory. It is planned to test the code on a Texas Instruments TMS320TCI6486, a 6-core SMP DSP with a 500MHz clock. Each core of this platform also has a private memory. Thus, application data and code can be removed from the shared memory and managed automatically using dataflow graph information. The TMS320TCI6486 enables research towards a RTOS for non-SMP MPSoC architectures.

VIII. CONCLUSION

In this paper, we detailed an RTOS for Symmetric Multiprocessing MPSoC combining ideas from dataflow models and from commonly-used RTOS. The advantages of such a system are numerous. Dataflow graphs bring automatic parallelization and ease of application description. The traditional difficulty of manually synchronizing tasks disappears. Moreover, the actors themselves can be written in C or C++ code, so legacy code can be easily reused. The system obtained is portable to SMP architecture with considerable number of cores.

The MPSoC RTOS obtained goes successively through a master/slave phase that parameterizes the dataflow graph and a symmetric phase that permits any core to access the RTOS scheduler. This allows the reactivity to events of an RTOS to be combined with the automatic partitioning and portable execution of an algorithm described in a parameterized dataflow graph. Experimental results are shown on a 3GPP Long Term Evolution algorithm running on a 4-core MPSoC. The MPSoC RTOS prototype is built on the µC/OS-II RTOS.

Based on the present prototype of the MPSoC RTOS, further research will be conducted to identify the dataflow models best suited to signal processing applications. Additional scheduling methods and architectures will also be tested to exploit a prior knowledge on graph execution.

Fig. 1 .

 1 Fig. 1. Structure of the MPSoC Model-Based RTOS Next sections detail the dataflow management and RTOS scheduling modules.

Fig. 2 .

 2 Fig. 2. Expanding the input dataflow model.

  Find highest priority ready task T If T is not the current taskSwitch context Execute the task preempted?

Fig. 5 .

 5 Fig. 5. Memory Usage

Fig. 3 .

 3 Fig. 3. Master/Slave Scheduling Phase: Expanding Dataflow Graph and Launching Slave Cores.

Fig. 6 .

 6 Fig. 6. Execution time performances

  .

	Iterations	1	2	3	4	5	6
	nbuser	10	3	6	8	10	2
	max nb CB per user	5	4	3	1	2	5
	max nb CB	14	10 13	5	8 14
	nbactors	39	31 43 27	30 36