Florence Levé

Richard Groult

Guillaume Arnaud

Cyril Séguin

Rémi Gaymay

Mathieu Giraud

RHYTHM EXTRACTION FROM POLYPHONIC SYMBOLIC MUSIC

In this paper, we focus on the rhythmic component of symbolic music similarity, proposing several ways to extract a monophonic rhythmic signature from a symbolic polyphonic score. To go beyond the simple extraction of all time intervals between onsets (noteson extraction), we select notes according to their length (short and long extractions) or their intensities (intensity +/-extractions). Once the rhythm is extracted, we use dynamic programming to compare several sequences. We report results of analysis on the size of rhythm patterns that are specific to a unique piece, as well as experiments on similarity queries (ragtime music and Bach chorale variations). These results show that long and intensity + extractions are often good choices for rhythm extraction. Our conclusions are that, even from polyphonic symbolic music, rhythm alone can be enough to identify a piece or to perform pertinent music similarity queries, especially when using wise rhythm extractions.

INTRODUCTION

Music is composed from rhythm, pitches, and timbres, and music is played with expression and interpretation. Omitting some of these characteristics may seem unfair. Can the rhythm alone be representative of a song or a genre? Small rhythmic patterns are essential for the balance of the music, and can be a way to identify a song. One may first think of some clichés: start of Beethoven 5th symphony, drum pattern from We will rock you or Ravel's Boléro. More generally, Query By Tapping (QBT) studies, where the user taps on a microphone [START_REF] Hanna | Query by tapping system based on alignment algorithm[END_REF][START_REF] Jang | Query by Tapping: a new paradigm for content-based music retrieval from acoustic input[END_REF], are able in some situations to identify a monophonic song. On a larger scale, musicologists have studied how rhythm, like tonality, can structure a piece at different levels [START_REF] Cooper | The Rhythmic Structure of Music[END_REF][START_REF] Rigaudière | La théorie musicale germanique du XIXe siècle et l'idée de cohérence[END_REF].

This article shows how simple extractions can, starting from a polyphony, build relevant monophonic signatures, being able to be used for the identification of songs or for the comparison of whole pieces.

Permission to make digital or hard copies of all or part of this work for personal or classroom use is granted without fee provided that copies are not made or distributed for profit or commercial advantage and that copies bear this notice and the full citation on the first page. c 2011 International Society for Music Information Retrieval.

In fact, most rhythm-only studies in Music Information Retrieval (MIR) concern audio signal. These techniques often rely in detection of auto-correlations in the signal. Some studies output descriptors [START_REF] Gruhne | Improving rhythmic similarity computation by beat histogram transformations[END_REF][START_REF] Peeters | Rhythm classification using spectral rhythm patterns[END_REF][START_REF] Robine | Meter class profiles for music similarity and retrieval[END_REF] that can be used for further retrieval or classification. Several papers focus on applications of non-Western music [START_REF] Holzapfel | Rhythmic similarity in traditional turkish music[END_REF][START_REF] Jensen | Rhythm-based segmentation of popular chinese music[END_REF][START_REF] Wright | Analyzing afro-cuban rhythms using rotation-aware clave template matching with dynamic programming[END_REF].

There are other tools that mix audio with symbolic data, comparing audio signals against symbolic rhythmic pattern. For example, the QBT wave task of MIREX 2010 proposed the retrieval of monophonic MIDI files from wave input files. Some solutions involve local alignments [START_REF] Hanna | Query by tapping system based on alignment algorithm[END_REF]. Another problem is rhythm quantization, for example when aligning audio from music performances against symbolic data. This can be solved with probabilistic frameworks [START_REF] Cemgil | Rhythm Quantization for Transcription[END_REF]. Tempo and beat detection are other situations where one extracts symbolic information from audio data [START_REF] Dixon | Automatic extraction of tempo and beat from expressive performances[END_REF][START_REF] Seyerlehner | From rhythm patterns to perceived tempo[END_REF].

Some rhythm studies work purely on symbolic MIDI data, but where the input is not quantized [START_REF] Lima | On rhythmic pattern extraction in Bossa Nova music[END_REF], as in the QBT symbolic task in MIREX 2010. Again, challenges can come from quantization, tempo changing and expressive interpretations. Finally, on the side of quantized symbolic music, the Mongeau and Sankoff algorithm takes into account both pitches and rhythms [START_REF] Mongeau | Comparaison of musical sequences[END_REF]. Extensions concerning polyphony have been proposed [START_REF] Allali | Toward a general framework for polyphonic comparison[END_REF]. Other symbolic MIR studies focus on rhythm [START_REF] Chen | Query by rhythm: An approach for song retrieval in music databases[END_REF][START_REF] Christodoulakis | Identifying rhythms in musical texts[END_REF][START_REF] Thul | Rhythm complexity measures: A comparison of mathematical models of human perception and performance[END_REF][START_REF] Toussaint | The geometry of musical rhythm[END_REF][START_REF] Godfried | A comparison of rhythmic similarity measures[END_REF].

However, as far as we know, a framework for rhythmic extraction from polyphonic symbolic music has never been proposed. Starting from a polyphonic symbolic piece, what are the pertinent ways to extract a monophonic rhythmic sequence? Section 2 presents comparison of rhythmic sequences through local alignment, Section 3 proposes different rhythm extractions, and Section 4 details evaluations of these extractions for the identification of musical pieces with exact pattern matching (Section 4.2) and on similarity queries between complete pieces (Sections 4.3 and 4.4).

RHYTHM COMPARISONS

Representation of monophonic rhythm sequences

For tempo-invariance, several studies on tempo or beat tracking on audio signal use relative encoding [START_REF] Hanna | Query by tapping system based on alignment algorithm[END_REF]. As we start from symbolic scores, we suppose here that the rhythms are already quantized on beats, and we will not study tempo and meter parameters. If necessary, multiple queries handle the cases where the tempo is doubled or halved.

Rhythm can be represented in different ways. Here, we model each rhythm as a succession of durations between notes, i.e. inter-onset intervals measured in quarter notes or fractions of them (Figure 1). Thus, in this simple framework, there are no silences, since each note, except the last one, is considered until the beginning of the following note.

Monophonic rhythm comparison

Several rhythm comparisons have been proposed [START_REF] Godfried | A comparison of rhythmic similarity measures[END_REF]. Here, we compare rhythms while aligning durations. Let S(m, n) be the best score to locally align a rhythm sequence x 1 . . . x m to another one y 1 . . . y n . This similarity score can be computed via a dynamic programming equation (Figure 2), by discarding the pitches in the Mongeau-Sankoff equation [START_REF] Mongeau | Comparaison of musical sequences[END_REF]. The alignment can then be retrieved through backtracking in the dynamic programming table. There can be a match or a substitution (s) between two durations, an insertion (i) or a deletion (d) of a duration. The consolidation (c) operation consists in grouping several durations into a unique one, and the fragmentation (f) in splitting a duration into several ones (see Figure 3). (

S(a, b) = max                                          S(a -1, b -1) + δ(x a , y b ) (match, substitution s) S(a -1, b) + δ(x a , ∅) (insertion i) S(a, b -1) + δ(∅, y b ) (deletion d) S(a -k, b -1) + δ({x a-k+1 ...x a }, y b ) (consolidation c) S(a -1, b -k) + δ(x a , {y b-k+1 ...y b }) (fragmentation f ) 0 (local alignment)
) (114) c. ( 143
) (53) d. (7) 43
(

) (10) f. (0) 6
Music engraving by LilyPond 2.12.3-www.lilypond.org Matches, consolidations and fragmentations respect the beats and the strong beats of the measure, whereas substitutions, insertions and deletions may alter the rhythm structure and should be more highly penalized. Scores will be evaluated in Section 4 where it is confirmed that, most of the time, the best results are obtained when taking into account consolidation and fragmentation operations.

RHYTHM EXTRACTION

How can we extract, from a polyphony, a monophonic rhythmic texture? In this section, we propose several rhythmic extractions. Figure 4 presents an example applying these extractions on the beginning of a chorale by J.-S. Bach.

The simplest extraction is to consider all onsets of the song, reducing the polyphony to a simple combined monophonic track. This "noteson extraction" extracts durations from the inter-onset intervals of all consecutive groups of notes. For each note or each group of notes played simultaneously, the considered duration is the time interval between the onset of the current group of notes and the following onset. Each group of notes is taken into account and is represented in the extracted rhythmic pattern. However, such a noteson extraction is not really representative of the polyphony: when several notes of different durations are played at the same time, there may be some notes that are more relevant than others.

In symbolic melody extraction, it has been proposed to select the highest (or the lowest) pitch from each group of notes [START_REF] Uitdenbogerd | Music Information Retrieval Technology[END_REF]. Is it possible to have similar extractions when one considers the rhythms? The following paragraphs introduce several ideas on how to choose onsets and durations that are most representative in a polyphony. We will see in Section 4 that some of these extractions bring a noticeable improvement to the noteson extraction.

Considering length of notes: long, short

Focusing on the rhythm information, the first idea is to take into account the effective lengths of notes. At a given onset, for a note or a group of notes played simultaneously:

• in the long extraction, all events occurring during the length of the longest note are ignored. For example, as there is a quarter on the first onset of Figure 4, the second onset (eighth, tenor voice) is ignored; • similarly, for the short extraction, all events occurring during the length of the shortest note are ignored. This extraction is often very close to the noteson extraction.

In both cases, as some onsets may be skipped, the considered duration is the time interval between the onset of the current group of notes and the following onset that is not ignored. Most of the time, the short extraction is not very different from the noteson, whereas the long extraction brings significant gains in similarity queries (see Section 4).

Considering intensity of onsets: intensity +/-

The second idea is to consider a filter on the number of notes at the same event, keeping only onsets with at least k notes (intensity + ) or strictly less than k notes (intensity -), where the threshold k is chosen relative to the global intensity of the piece. The considered durations are then the time intervals between consecutive filtered groups. Figure 4 shows an example with k = 3. This extraction is the closest to what can be done on audio signals with peak detection.

RESULTS AND EVALUATION

Protocol

Starting from a database of about 7000 MIDI files (including 501 classical, 527 jazz/latin, 5457 pop/rock), we selected the quantized files by a simple heuristic (40 % of onsets on beat, eighth or eighth tuplet). We thus kept 5900 MIDI files from Western music, sorted into different genres (including 204 classical, 419 jazz/latin, 4924 pop/rock). When applicable, we removed the drum track (MIDI channel 10) to avoid our rhythm extractions containing too many sequences of eighth notes, since drums often have a repetitive structure in popular Western music. Then, for each rhythm extraction presented in the previous section, we extracted all database files. For each file, the intensity +/-threshold k was choosen as the median value between all intensities. For this, we used the Python framework music21 [START_REF] Scott | music21: A toolkit for computer-aided musicology and symbolic music data[END_REF].

Our first results are on Exact Song identification (Section 4.2). We tried to identify a song by a pattern of several consecutive durations taken from a rhythm extraction, and looked for the occurrences of this pattern in all the songs of the database.

We then tried to determinate if these rhythm extractions are pertinent to detect similarities. We tested two particular cases, Ragtime (Section 4.3) and Bach chorales variations (Section 4.4). Both are challenging for our extraction methods, because they present difficulties concerning polyphony and rhythm: Ragtime has a very repetitive rhythm on the left hand but a very free right hand, and Bach chorales have rhythmic differences between their different versions.

Exact Song Identification

In this section, we look for patterns of consecutive notes that are exactly matched in only one file among the whole database. For each rhythm extraction and for each length between 5 and 50, we randomly selected 200 distinct patterns appearing in the files of our database. We then searched for each of these patterns in all the 5900 files (Figure 5). We see that as soon as the length grows, the patterns are very specific. For lengths 10, 15 and 20, the number of patterns (over 200) matching one unique file is as follows: We notice that the long and intensity +/-extractions are more specific than noteson. From 12 notes, the median values of Figure 5 are equal to 1 except for noteson and short extractions. In more than 70% of these queries, 15 notes are sufficient to retrieve a unique file.

The results for average values are disturbed by a few patterns that match a high number of files. Figure 6 displays some noteworthy patterns with 10 notes. Most of the time, the patterns appearing very frequently are repetitions of the same note, such as pattern (a). With long extraction, 174 files contain 30 consecutive quarters, and 538 files contain 30 consecutive eighths. As these numbers further increase with noteson (and short) extractions, this explains why the long extraction can be more specific. (

) (114) c. ( 143
) (53) d. (7) 43
(6)

(0) (1) g.

Music engraving by LilyPond 2.12.3-www.lilypond.org The number of occurrences of each pattern is mostly determined by its musical relevance. For example, in a pattern with three durations, (d) appears more often than (g), which is quite a difficult rhythm. In the same way, among patterns with only quarters and eighths, (b) and (c) can be found more often than (f). We also notice that patterns with longer durations, even repetitive ones such as pattern (e), generally appear in general less frequently than those containing shorter durations.

Similarities in Ragtime

In this section and the following, we use the similarity score computation explained in Section 2.2. Ragtime music, one of the precursors of Jazz music, has a strict tempo maintained by the pianist's left hand and a typical swing created by a syncopated melody in the right hand.

For this investigation, we gathered 17 ragtime files. Then we compared some of these ragtime files against a set of files comprising the 17 ragtime files and randomly selected files of the database. We tested several scores functions: always +1 for a match, and -10, -5, -2, -1, -1/2 or -1/3 for an error. We further tested no penalty for consolidation and fragmentation (c/f ). Figure 7 shows ROC Curves for A Ragtime Nightmare. A ROC Curve [START_REF] Fawcett | An introduction to ROC analysis[END_REF] plots sensibility (capacity to find true positives) and specificity (capacity to eliminate false positives) over a range of thresholds, giving a way to ascertain the performance of a classifier that outputs a ranked list of results.

Here one curve represents one rhythm extraction with one score function. For each score function, we computed the true positive and the false positive rates according to all different thresholds. The long extraction, used with scores +1 for a match and -1 for all errors, gives here very good results: for example, the circled point on Figure 7 corresponds to 0.88 sensitivity and 0.84 specificity with a threshold of 45 (i.e. requiring at least 45 matches).

Considering the whole curve, the performance of such a classifier can be measured with the AUC (Area Under ROC Curve). Averaging on 9 different queries, the best set of scores for each extraction is as follows:

Extraction Scores Mean AUC s/i/d c/f match noteson -5 0 +1 0.711 short -1 -1 +1 0.670 long -1 0 +1 0.815 intensity + -1/3 0 +1 0.622 intensity - -1 -1 +1 0.697
Most of the time, the matching sequences are long sequences of eighths, similar to pattern (a) of Figure 6. If such patterns are frequent in noteson database files (see previous section), their presence in long files is more frequent in Ragtime than in other musical styles. For example, pattern (a) is found in 76 % of Ragtime long extractions, compared to only 25 % of the whole database.

Indeed, in ragtime scores, the right hand is very swift and implies a lot of syncopations, while the left hand is bet- ter structured. Here the syncopations are not taken into account in the long extraction, and the left hand (often made of eighths, as in Figure 8) is preserved during long extractions.

Finally, intensity + does not give good results here (unlike Bach Chorales, see next Section). In fact, intensity + extraction keeps the syncopation of the piece, as accents in the melody often involve chords that will pass through the intensity + filter (Figure 8, last note of intensity + ).

Similarities in Bach Chorales Variations

Several Bach chorales are variations of each other, sharing an exact or very similar melody. Such chorales present mainly variations in their four-part harmony, leading to differences in their subsequent rhythm extractions (Figure 9). For this investigation, we considered a collection of 404 Bach chorales transcribed by www.jsbchorales.net and available in the music21 corpus [START_REF] Scott | music21: A toolkit for computer-aided musicology and symbolic music data[END_REF]. We selected 5 chorales that have multiple versions: Christ lag in Todesbanden (5 versions, including a perfect duplicate), Wer nun den lieben Gott (6 versions), Wie nach einer Wasserquelle (6 versions), Herzlich tut mich verlangen (9 versions), and O Welt, ich muss dich lassen (9 versions).

For each chorale, we used one version to query against the set of all other 403 chorales, trying to retrieve the most similar results. A ROC curve with BWV 278 as a query is shown in Figure 10. For example, with intensity + extraction and scores -1 for s/i/d, 0 for c/f , and +1 for a match, the circled point corresponds to a threshold of 26, with 0.80 sensitivity and 0.90 specificity. Averaging on all 5 chorales, the best set of scores for each extraction is as follows: We thus managed to align almost globally one chorale and its variations. We further checked that there is not a bias on total length: for example, BWV 278 has a length of exactly 64 quarters, as do 15% of all the chorales, but the score distribution is about the same in these chorales than in the other ones.

DISCUSSION

In all our experiments, we showed that several methods are more specific than a simple noteson extraction (or than the similar short extraction). The intensity -extraction could provide the most specific patterns used as signature (see Figure 5), but is not appropriate to be used in similarity queries. The long and intensity + extractions give good results in the identification of a song, but also in similarity queries inside a genre or variations of a music.

It remains to measure what is really lost by discarding pitch information: our perspectives include the comparison of our rhythm extractions with others involving melody detection or drum part analysis.

Figure 1 .

 1 Figure 1. The monophonic rhythm sequence (1, 0.5, 0.5, 2).

Figure 2 .

 2 Figure 2. Dynamic programming equation for finding the score of the best local alignment between two monophonic rhythmic sequences x 1 . . . x a and y 1 . . . y b . δ is the score function for each type of mutation. The complexity of computing S(m, n) is O(mnk), where k is the number of allowed consolidations and fragmentations.

Figure 3 .

 3 Figure 3. Alignment between two rhythm sequences.

Figure 4 .

 4 Figure 4. Rhythm extraction on the beginning of the Bach chorale BWV 278.

Figure 5 .

 5 Figure 5. Number of matching files for patterns between length 5 and 35. Curves with points indicate median values, whereas other curves indicate average values.

Figure 6 .

 6 Figure 6. patterns with 10 durations, with the number of matching files in noteson and long extractions.

Figure 7 .

 7 Figure 7. Best ROC Curves, with associated AUC, for retrieving 17 Ragtime pieces from the query A Ragtime Nightmare, by Tom Turpin, in a set of 100 files.

Figure 8 .

 8 Figure 8. Possum Rag (1907), by Geraldine Dobyns.

Figure 9 .

 9 Figure 9. Extraction of long rhythm sequences from different variations of the start of the chorale Christ lag in Todesbanden. The differences between variations are due to differences in the rhythms of the four-part harmonies.

Figure 10 .

 10 Figure 10. Best ROC Curves, with associated AUC, for retrieving all 5 versions of Christ lag in Todesbanden from BWV 278 in a set of 404 chorales.

  Even if the noteson extractions already gives good results, long and intensity + bring noteworthy improvements. Most of the time, the best scores correspond to alignments between 8 and 11 measures, spanning a large part of the chorales.

	Extraction		Scores		Mean AUC
		s/i/d c/f match
	noteson	-1	0	+1	0.769
	short	-1	0	+1	0.781
	long	-5	-5	+1	0.871
	intensity +	-1	0	+1	0.880
	intensity -	-5	0	+1	0.619

Acknowledgements. The authors thank the anonymous referees for their valuable comments. They are also indebted to Dr. Amy Glen who kindly read and corrected this paper.


M Léonard

L Mouchard

M Salson

On the Number of Elements to Reorder When Updating a Suffix Array

Keywords: FM-index, Burrows-Wheeler Transform, Suffix Array, Longest Common Prefix, Dynamic Data Structure, Self-Index, Complexity

Recently new algorithms appeared for updating the Burrows-Wheeler transform or the suffix array, when the text they index is modified. These algorithms proceed by reordering entries and the number of such reordered entries may be as high as the length of the text. However, in practice, these algorithms are faster for updating the Burrows-Wheeler transform or the suffix array than the fastest reconstruction algorithms.

In this article we focus on the number of elements to be reordered for real-life texts. We show that this number is related to LCP values and that, on average, L ave entries are reordered, where L ave denotes the average LCP value, defined as the average length of the longest common prefix between two consecutive sorted suffixes. Since we know little about the LCP distribution for real-life texts, we conduct experiments on a corpus that consists of DNA sequences and natural language texts. The results show that apart from texts containing large repetitions, the average LCP value is close to the one expected on a random text.

Introduction

The Burrows-Wheeler Transform [START_REF] Burrows | A block-sorting lossless data compression algorithm[END_REF][START_REF] Ferragina | The Burrows-Wheeler Transform (special issue)[END_REF], BWT, is a very interesting blocksorting algorithm that reorders the letters of a text T of length n over an alphabet of size σ for easing its compression. It is used as a preprocessor by some of the most popular lossless text compression tools (such as bzip) that chain it to Run-Length Encoding, entropy encoding or Prediction by Partial Matching methods [START_REF] Cleary | Data compression using adaptive coding and partial string matching[END_REF][START_REF] Cleary | Unbounded length contexts for PPM[END_REF].

Conceptually speaking, the suffix array [START_REF] Manber | Suffix arrays: a new method for on-line string searches[END_REF][START_REF] Gonnet | New indices for text: Pat trees and pat arrays[END_REF] is very close to the text produced by the BWT. Due to its intrinsic structure and its similarity with the suffix array, it has also been used for advanced compressed index structures [START_REF] Ferragina | Opportunistic data structures with applications[END_REF][START_REF] Ferragina | Compressed representation of sequences and full-text indexes[END_REF], such as FM-index, that authorize approximate pattern matching [START_REF] Lam | Compressed indexing and local alignment of DNA[END_REF][START_REF] Langmead | Ultrafast and memoryefficient alignment of short DNA sequences to the human genome[END_REF], and therefore can be used by search engines.

The FM-index is based on the close relationship between the Burrows-Wheeler Transform and the suffix array. It is therefore a sort of compressed suffix array that takes advantage of the compressibility of the indexed text in order to reduce space occupancy with respect to the entropy of the text.

Recently algorithms that update BWT [START_REF] Salson | A four-stage algorithm for updating a Burrows-Wheeler Transform[END_REF], corresponding enhanced/extended suffix array [START_REF] Gallé | In-place update of suffix array while recoding words[END_REF] or FM-index [START_REF] Salson | Dynamic extended suffix array[END_REF] appeared when edit operations transform the text T into a text T ′ of length n ′ . At first sight, the overall complexity of such algorithms is bounded by O(n log n log σ) for [START_REF] Salson | A four-stage algorithm for updating a Burrows-Wheeler Transform[END_REF][START_REF] Salson | Dynamic extended suffix array[END_REF], and O(n 2 ) for [START_REF] Gallé | In-place update of suffix array while recoding words[END_REF] making them slower than known linear-time suffix array construction algorithms [START_REF] Kärkkäinen | Simple linear work suffix array construction[END_REF][START_REF] Ko | Space efficient linear time construction of suffix arrays[END_REF].

Nevertheless, practical experiments conducted in the past (see [START_REF] Gallé | In-place update of suffix array while recoding words[END_REF][START_REF] Salson | A four-stage algorithm for updating a Burrows-Wheeler Transform[END_REF][START_REF] Salson | Dynamic extended suffix array[END_REF]) have shown that they globally outperform the quickest known suffix array construction algorithm [START_REF] Puglisi | A taxonomy of suffix array construction algorithms[END_REF] for a reasonable number of edit operations.

The time consuming part when one wants to update a suffix array-related structure is to reorder suffixes. In the worst case a linear number of suffixes can be reordered. In this article, we show that the average time complexity is strongly connected to the average LCP , the longest common prefix. For various texts, such as repetitive genomic sequences or natural language texts, these values are surprisingly small.

In section 2 we briefly sketch properties of the suffix array and the Burrows-Wheeler Transform, we identify four types of suffixes (or cyclic shifts) and, based on that, characterize the number of elements to reorder on average. In section 3 we present the sets of data that have been used for the tests, explain their relevance. In section 4 we present practical results together with discussions and finally we conclude in section 5.

Preliminaries and Algorithms

We consider the text T = T [0 . . n], a word of length n + 1 over Σ, a finite ordered alphabet of size σ. The last letter of T is a sentinel letter $, that has been added to the alphabet Σ and is smaller than any other letter of Σ. A factor starting at position i and ending at position j is denoted by T [i . . j] and a single letter is denoted by T [i] (or T i to facilitate the reading). We add that when i > j, T [i . . j] is the empty word. A factor starting at position 0 is a prefix of T while a factor ending at position n is a suffix of T .

We briefly recall well-known definitions that will be used throughout the document.

The cyclic shift of order i of T is T [i] = T [i . . n]T [0 . . i -1] for 0 ≤ i ≤ n. The Burrows-Wheeler Transform of T , denoted BWT (T ), is the text of length n + 1 corresponding to the last column L of the conceptual matrix whose rows are the lexicographically sorted cyclic shifts T [i] (see Fig. 1(b)). The suffix array of T , denoted by SA(T ), or simply SA, is the list of starting positions of lexicographically sorted suffixes of T (or cyclic shifts of T ). That is SA[j] = i if and only if T [i . . n] is the (j + 1) th suffix of T in ascending order. The inverse of SA, denoted by ISA(T ), indicates for each suffix of T its corresponding position in SA. For each i ∈ [0, n -1], LCP [i] is the length of the Longest Common Prefix between suffixes starting at positions SA[i] and SA[i + 1] in the text. It is clear that BWT , SA and ISA are strongly related, even if BWT is made of letters, with potential multiplicity, while SA and ISA are made of unique consecutive integers. In [START_REF] Salson | A four-stage algorithm for updating a Burrows-Wheeler Transform[END_REF][START_REF] Salson | Dynamic extended suffix array[END_REF] we explained how the FM-index, which is based on these three data structures, can be updated rather than entirely reconstructed when the text T is edited (insertion, deletion or substitution of a letter or a factor). Gallé et al. [START_REF] Gallé | In-place update of suffix array while recoding words[END_REF] explained how a suffix array can be updated when indexing for grammatical inference purposes.

T = C 0 T 1 A 2 G 3 T 4 T 5 A 6 G 7 $ 8 T [0] C T A G T T A G $ T [1] T A G T T A G $ C T [2] A G T T A G $ C T T [3] G T T A G $ C T A T [4] T T A G $ C T A G T [5] T A G $ C T A G T T [6] A G $ C T A G T T T [7] G $ C T A G T T A T [8] $ C T A G T T A G (a) Cyclic shifts F L T [8] $ C T A G T T A G T [6] A G $ C T A G T T T [2] A G T T A G $ C T T [0] C T A G T T A G $ T [7] G $ C T A G T T A T [3] G T T A G $ C T A T [5] T A G $ C T A G T T [1] T A G T T A G $ C T [4] T T A G $ C T A G (b)
In the following we consider indifferently cyclic shifts or suffixes. When a text is updated, in both cases, some suffixes (or cyclic shifts) have to be moved to take into account their new lexicographical order. We will use the word "reordering" whenever a new order has to be assigned (insertion of a new position that shifts everything by one position, rotation, ...).

In this article, we focus on the number of elements that have to be reordered when one wants to update a Burrows-Wheeler Transform or a suffix array. This analysis mainly applies to [START_REF] Salson | A four-stage algorithm for updating a Burrows-Wheeler Transform[END_REF][START_REF] Salson | Dynamic extended suffix array[END_REF] but are of interest for any algorithm that updates dynamically a BWT, a suffix array [START_REF] Gallé | In-place update of suffix array while recoding words[END_REF] or a more complex stucture such as FM-indexes. Moreover, it greatly helps in understanding the behavior of several SACA, as presented in [START_REF] Puglisi | A taxonomy of suffix array construction algorithms[END_REF].

SA 8 T [8] $ C T A G T T A G 6 T [6] A G $ C T A G T T 2 T [2] A G T T A G $ C T 0 T [0] C T A G T T A G $ 7 T [7] G $ C T A G T T A 3 T [3] G T T A G $ C T A 5 T [5] T A G $ C T A G T 1 T [1] T A G T T A G $ C 4 T [4] T T A G $ C T A G F L (a) BWT (CTAGTTAG$) SA Type $ C T A C G T T A G T ′[9] 9 1 A C G T T A G $ C T T ′[2] 2 4 A G $ C T A C G T T T ′[7] 7 1 C G T T A G $ C T A T ′[3] 3 3 C T A C G T T A G $ T ′[0] 0 4 G $ C T A C G T T A T ′[8] 8 1 G T T A G $ C T A C T ′[4] 4 2 T A C G T T A G $ C T ′[1] 1 4 T A G $ C T A C G T T ′[6] 6 1 T T A G $ C T A C G T ′[5] 5 1 F L (b) BWT (CTACGTTAG$)
Fig. 2. T , T ′ and their respective SA and BWT

Without loss of generality we will consider the insertion of a single letter in T at position i. Let T ′ be the final text obtained after modifying T .

For updating the suffix array or the Burrows-Wheeler Transform, we need to differentiate four types of rows in the conceptual matrix. They correspond to the position of the inserted letter in the cyclic shifts.

During the updating phase, the cyclic shifts are processed from Type 1 to Type 4.

Type 1 the inserted letter appears after $ and before L. The number of these rows is variable: the closer to the beginning of the text the edit operation occurs, the larger this number is. We showed in [START_REF] Salson | A four-stage algorithm for updating a Burrows-Wheeler Transform[END_REF][START_REF] Salson | Dynamic extended suffix array[END_REF] that the respective order of these rows is preserved, they are only subject to implicit reordering. See dashed arrows in Fig. 2.

It corresponds to all cyclic shifts T ′[n+1] to T ′[i+2] in decreasing order. Type 2 the inserted letter appears in L. An existing row is modified, the letter appearing in L is stored (in our example A) and replaced by the new letter (in our example C). See dotted arrow in Fig. 2. It corresponds to the unique cyclic shift T ′[i+1] . Type 3 the inserted letter appears in F . A new row is inserted, the new letter is inserted in column F (in our example C) and previously stored letter in column L (in our example A). Since a new row is inserted, indices of all shifted rows are incremented. It corresponds to the unique cyclic shift T ′ [i] . Type 4 the inserted letter appears after F and before $. Each row has a potential impact on the respective order of the other cyclic shifts. It induces a partial reordering of these cyclic shifts that should be performed if needed. See plain arrows in Fig. 2.

It corresponds to all cyclic shifts T ′[i-1] to T ′[0] in decreasing order.

We already proved ( [START_REF] Salson | A four-stage algorithm for updating a Burrows-Wheeler Transform[END_REF]) that only rows of Type 4 are subject to direct reordering, playing a central role in the overall time complexity. The main reason is that the modification appears before $ in those cyclic shifts. Hence that may modify their lexicographical ranking. Since we focus on the number of reorderings, we restrict our study to Type 4 cyclic shifts.

For a given modification, directly reordered cyclic shifts are either moved up or down. Without loss of generality, we consider the case where cyclic shifts are moved up.

We consider the j-th cyclic shift of Type 4, that is T ′[i-j] . We denote by k the position of this cyclic shift in the conceptual matrix and ℓ the corresponding

LCP value, that is ℓ = LCP[k -1].
Lemma 1 When j > ℓ, T ′[i-j] is not directly reordered and the reordering stage ends.

PROOF. When we are processing the j-th cyclic shift of Type 4, T ′[i-j] , all the cyclic shifts T ′ [r] such that r > ij have already been ordered. Suppose that j > ℓ.

Let p = SA[k -1], we also know that SA[k] = ij. For the sake of clarity, we also note

c = T ′[p] [ℓ] and c ′ = T ′[i-j] [ℓ].
The modification appears at position j in T ′[i-j] , hence:

T ′[i-j] [0 . . ℓ] = T [i-j] [0 . . ℓ], since j > ℓ (1)
Following the definitions of LCP and BWT, T [i-j] and T ′[p] share a common prefix of length ℓ and T [i-j] [ℓ] > c. Using (1), we also know that this relation still holds for T ′[i-j] : prefixes of length ℓ of T ′[i-j] and T ′ [p] are equal and c ′ > c (see Fig. 3). -j] does not have to be reordered.

i SA k -1 p c ℓ < T ′[p] k i -j c ′ 0 j n ′ Modification T ′[i-j] Fig. 3. T ′[i
We showed that T ′[i-j] is correctly ordered and has not to be reordered. We proved [21, Lemma 3] that if T ′[i-j] is correctly ordered then all the cyclic shifts T ′[q] , q < ij are correctly ordered as well. Finally, if j > LCP[k -1] then the reordering ends.

It gives us a very interesting upper bound on the number of reorderings that should be performed after a given edit operation.

Lemma 2

The average number of elements to be reordered, after one given edit operation, is at most equal to L ave , the average LCP value.

PROOF. For a modification occurring at position i in T , we denote by r i the number of elements to be reordered and by r the array of all r i , for 0 ≤ i ≤ n.

We also denote by PLCP the array which is the permuted LCP array where values appear in position order rather than lexicographical order. PLCP array has the following property [START_REF] Kasai | Linear-time longestcommon-prefix computation in suffix arrays and its applications[END_REF]:

PLCP[i -1] -1 ≤ PLCP[i], ∀i ∈ [0 . . n -1].
We insert a 0 at the beginning of the array such that PLCP[-1] = 0. By definition PLCP[n -1] = 0. We will show that r is a permutation of PLCP.

From the previous lemma,

PLCP[i] = 0 ⇐⇒ r i+1 = 0. Let us consider PLCP[j . . j ′ ],
where for any k such that j ≤ k ≤ j ′ , we have PLCP[k] > 0, and Finally, we have PLCP

PLCP[j -1] = 0, PLCP[j ′ + 1] = 0. Let PLCP[j] = L 1 , L 1 >
[j . . j ′ ] = [L 1 , L 1 -1, . . , ℓ 1 , L 2 , L 2 -1, . . , ℓ p-1 , L p , . . , 1].
Lemma 1 can be rewritten as follows:

r i =      0 if PLCP[i -1] = 0 d otherwise, with 1 ≤ k ≤ d, PLCP[i -k] ≥ k and PLCP[i -(d + 1)] < d + 1
(2)

j -1 j + L 1 -ℓ 1 L 1 ℓ 3 ℓ 2 L 1 ℓ 1 L 2 ℓ 2 L 3 ℓ 3 L 4 ℓ 4 i d + 1 k Fig. 4. Computation of r i = d using the PLCP array, we have ℓ 3 ≥ k and ℓ 2 < d + 1.
Indeed, a suffix might be reordered whenever the modification belongs to a prefix which appears more than once. When it is not the case anymore (ie. 4), the corresponding suffix does not have to be reordered.

PLCP[i -(d + 1)] < d + 1, see Fig.
Since PLCP[j -1] = 0, r j = 0. By definition r i+1 ≤ r i + 1, hence r j+1 = 1, r j+2 = 2 and so forth until position j + L 1 , where

r j+L 1 = L 1 . The next value r j+L 1 +1 is ℓ 1 (if ℓ 1 = L 1 ) since PLCP[j + L 1 + 1 -(ℓ 1 + 1)] = PLCP[j +L 1 -ℓ 1 ] = ℓ 1 < ℓ 1 +1 (see equation (2)
). Then the sequence increases again until reaching L 2 and so forth.

In the special case where ℓ i = L i (the decreasing sequence is therefore of length 1), with 1 ≤ i < p, let r k = L i , then r k+1 = L i + 1 and the sequence keeps increasing until L i+1 .

Finally

r[j + 1 . . j ′ + 1] = [1, . . , L 1 , ℓ 1 , . . , L 2 , ℓ 2 , . . , L p ] and it is a permutation of PLCP[j . . j ′ ].
We have shown this, considering a 0-free factor of PLCP . Therefore this is true for each such factor. Moreover using the definition of r i , we know that we have as many 0's in r than in PLCP. Hence, r is a permutation of PLCP and the average number of elements to be reordered is equal to L ave .

Hence, when modifying T at position i, we have at most

k reordered elements if PLCP[i-(k +1)] < k +1.
Therefore LCP values can be used for determining an upper bound of the number of elements to directly reorder for updating a suffix array.

3 Materials and methods: experimental data and their relevance

Based on Lemma 2, we conduct experiments on various texts. Our purpose is to determine the distribution of LCP values and give a better insight of the average number of reorderings when one of the studied texts is updated.

In what follows:

• the percentile LCP value, denoted by L perc , is an upper bound of the number of elements to reorder in 99% of the cases; • the average LCP value, denoted by L ave , is an upper bound of the number of elements to reorder on average; • the maximal LCP value, denoted by L max , is an upper bound of the number of elements to reorder, in any case.

Note that these values are upper-bounds since the number of elements to reorder also depends on the inserted (or deleted) letters. For instance, let us

consider T = A 0 A 1 A 2 A 3 $ 4
, we have SA(T ) = 4 3 2 1 0. If one inserts a A at position 4 in T , the resulting suffix array is 5 4 3 2 1 0, no element was reordered. On the opposite if one inserts a G at position 4, the resulting suffix array is 5 0 1 2 3 4, which denotes that many elements were reordered.

We subsequently define R n = L ave /n as the ratio between the average LCP value and the length of the text. In the worst case, for T =A n , L ave =n/2 and R n =0.5.

Remark 3 This analysis is of interest for all algorithms that update BWT or SA by reordering elements. It is also of importance for understanding the behavior of suffix array construction algorithms as described by Puglisi et al. [START_REF] Puglisi | A taxonomy of suffix array construction algorithms[END_REF]. For two texts having the same length, Seward's algorithm [START_REF] Seward | On the performance of BWT sorting algorithms[END_REF] is slower when processing the text with the largest L ave . For solving the problem that arises while considering large LCP values, Manzini and Ferragina's algorithm [START_REF] Manzini | Engineering a lightweight suffix array construction algorithm[END_REF] performs a special processing for LCP values that are larger than a threshold chosen by the user. In order to determine that appropriate threshold, one needs to use the distribution of LCP values for various texts, with respect to their length.

In what follows, we consider two different classes of texts: genomic sequences (unstructured texts) and natural language texts (structured texts). Each studied text possesses specific properties such as number and length of repeats, entropy and alphabet size. These texts can be classified according to a measure of repetitiveness, such as the one introduced by Haubold and Wiehe [START_REF] Haubold | How repetitive are genomes?[END_REF]. This measure is based on the length of the shortest unique substring. They define their index of repetitiveness, denoted by I r , as the logarithm of the ratio of the sum of the shortest unique substring lengths at each position, and the theoretical expectation for a random text having the same composition (as described in [START_REF] Haubold | How repetitive are genomes?[END_REF]). We note that I r can be computed in a simpler way than the one described in their article, based on suffix trees: compute the LCP array for the original text (concatenated to its reverse complement in the case of a DNA sequence); traverse the array top to bottom; for each two consecutive values, pick the largest one and add it to the final result. In a similar manner, compute the sum for the randomly shuffled original text. The logarithm of the ratio of the two sums corresponds to the I r value. The authors mention they used a suffix tree, which is a very space-consuming structure, we instead consider the natural inherent LCP values extending the traditional suffix array. It permits a faster computation of the I r values for genome sequences they were not able to handle, due to memory limitations.

In [START_REF] Fayolle | Analysis of the average depth in a suffix tree under a Markov model[END_REF], Fayolle and Ward showed that the average LCP value, under a Markovian model of order one, is (log n)/H 1 (T ) + C, where H 1 (T ) is the first-order entropy of T and C is a constant. In order to study if our practical values are diverging significantly from the theoretical values, we are defining R FW = L ave × H 1 (T )/ log n. It is clearly a measure of proximity between the number of direct reorderings that should be performed on average and Fayolle and Ward's theoretical value.

We know that L ave ranges from O(log n) to n/2 depending on the texts. On one hand many articles show that on a randomly generated text, the average LCP is O(log n), on the other hand one can easily build a text whose L ave is n/2 (e.g. A n ). Since the texts we chose are neither random nor composed of a single letter, we only know that their L ave is between O(log n) and n/2. In this section, we want to determine the LCP distribution and, among that, the average LCP value on different classes of text containing different types of repetition. This will allow us to evaluate how far typical texts are from the two extremes and how many elements have to be reordered at most for updating a suffix array for these classes of text.

Genomic Sequences

More and more complete genome sequences are available [START_REF] Liolios | The genome on line database (GOLD) in 2007: status of genomic and metagenomic projects and their associated metadata[END_REF], Haubold and Wiehe selected a subset of 336 organisms (330 prokaryotes and 6 eukaryotes) and computed I r values for each genome sequence. They ranked the genome sequences, from the most repeated to the least repeated sequence 1 . Following this study, we focus on the five most repeated among the 330 prokaryote genomes (whose sequence lengths are close to 2 Mbp). The corresponding I r values range from 6.34 to 3.84, and are related to organisms: M. flagellatus, S. agalactiae, D. ethenogenes, F. tularensis, N. meningitidis.

We extended our study to encompass the six eukaryote organisms considered in [START_REF] Haubold | How repetitive are genomes?[END_REF]: mouse-ear cress (A. thaliana), nematode (C. elegans), fruit fly (D. melanogaster ), yeast (S. cerevisiae), human (H. sapiens) and mouse (M. musculus). The genome sequences are longer: the largest chromosome (chr. 1) of H. sapiens contains more than 200 Mbp, and are usually sheltering shorter repeats. These eukaryote organisms have been extensively studied, are well documented and annotated, being part of a larger set of eight model organisms.

In addition we also consider a set of 398 virus genomes (whose sequence length are between 1 kb and 1 Mb). We will study in details two representatives of this set: G. fumiferanae ichnovirus and C. herpesvirus.

To avoid biased results due to unsequenced portions of genomes, we remove all runs of N's from the sequences. It explains why, for H. sapiens, chr1 is shorter than chr2: we had to remove large runs of N's from chr1, while chr2 was not containing that many.

Natural Language Texts

In addition to these large DNA sequences, we retrieve natural language texts from project Gutenberg 2 (plain ASCII texts, no extra tags). We also retrieved some corpora 3 of the Wikipedia encyclopedia (XML document, Wiki tags, plagiarism) in various languages 4 (Afrikaans, Basque, Bosnian, Estonian, Latin, Occitan). From project Gutenberg, we obtained a single file of size 91,070,340 bytes by concatenating 224 files from which we removed specific Gutenberg headers and footers. Due to this processing, the size of the alphabet and the type of text considered, we do not expect to see large LCP values since authors are usually trying to avoid as much as possible repeats and plagiarism. In the following we refer to this text as "etext". Regarding Wikipedia, each corpus for a given language is stored as a XML document and the content of the encyclopedia is written using the wiki language. The underlying XML structure implies more redundancy, compared to plain natural language texts, namely several predefined templates that are specific to the wiki language for displaying information in a standard format. These templates are often large duplicated pieces of codes, explaining a lot of redundancy that directly impacts LCP values. To evaluate the effect of the templates on the LCP values, we considered two versions of the Afrikaans corpus: the original one and the one from which we removed the most repeated template (see Fig. 6). Moreover Wikipedia is known to be a cooperative encyclopedia. Some people creates redundant articles and start by filling them using the content of a previous article (e.g. content of "People's Republic of China" is used in "History of China", in the Basque version) leading to larger expected LCP values.

We compute suffix arrays and LCP tables for all the texts, using Manzini and Ferragina's algorithms [START_REF] Manzini | Engineering a lightweight suffix array construction algorithm[END_REF][START_REF] Manzini | Two space saving tricks for linear time LCP array computation[END_REF] and implementation5 . These algorithms are chosen for their low space-consumption properties. For each text, we calculate the values L max , L ave , L perc as defined previously (see Sect. 3, page 8). We also compute R n , R FW and I r . In the next section we will also plot curves corresponding to L ave and the actual number of reorderings with respect to suffix lengths, for several texts after the insertion of a single letter. We will use these graphs for studying how these values are correlated and are affected by repetitions.

Results

In this section we present experimental results obtained from the sets of sequences we previously described. For each set of texts we compute their characteristics and the actual number of reorderings that are needed for maintaining the associated dynamic suffix array. Furthermore we compare the computed L ave value with the real number of reorderings that had to be performed during our experiments on various texts, with various edit operations. We show that they are strongly correlated for all studied texts.

Genomic Sequences

We are here presenting two analyses which are based on the values we compute from the LCP array or on the curves we plot from the experiments we conducted.

LCP values analysis

The characteristics for the five most repeated prokaryote genomes are presented in Table 1. These short genome sequences have been extensively studied and are well known for containing a lot of short repeats as well as very large repeated regions with a small number of copies. It is therefore unsurprising to observe very large L max relatively to their length. As an example, in M. flagellatus 5% of the genome corresponds to a single repeat. The longest repeated region is 143 Kbp long while the second longest is less than 1 Kbp long. For M. flagellatus, the L ave value is artificially large because of the L max value. The largest I r value is achieved for M. flagellatus although it possesses only one huge repeat. Therefore we can wonder if a text having a unique huge repeat has to be considered as "more repeated" than a text having many mediumsize repeats. Since the average complexity of our algorithm depends on L ave , and we want to test its efficiency in the worst conditions, it is nevertheless pertinent to consider most-repeated texts in the sense of the I r definition. Although we have large LCP values on average, they are still far from the average number of reorderings in the worst case: n/2. Values R n show that on average 0.116% (M. flagellatus) to 0.0115% (N. meningitidis) elements have to be reordered instead of 50% (the worst case one could expect). R FW values are above 20, meaning that they are significantly different from the theoretical values one can expect from [START_REF] Fayolle | Analysis of the average depth in a suffix tree under a Markov model[END_REF] showing that a Markovian model of order one is not suitable for such sequences.

Since we selected the five most repeated prokaryote genome sequences, we do not expect to deal with sequences having larger I r values. Hence, the R FW values which appear to be significantly larger than still costs less than rebuilding the entire new suffix array using a linear time construction algorithm.

The eukaryote genomes like the ones considered in We remark that for the most repeated chromosome sequence from M. musculus, that is chr. Y, only 0.00709% elements might be reordered on average, while in the worst case only 16, 589/2, 702, 555 × 100 = 0.61% might be.

Organism Seq. Len. Similarly to what we observed for the bacteria and the eukaryotes, in Table 5, we remark that for G. fumiferanae ichnovirus, only 0.00769% elements might be reordered on average, while in the worst case only 0.56% might be.

L max L ave L perc R n R FW I r G. f. ichnovirus 291,
On average a very limited number of reorderings have to be performed with respect to the length of the genomic sequences. Even for the most repeated sequence (M. flagellatus) a maximal 5 % of the whole sequence should be reordered if needed.

Curves analysis

In Fig. 5 we trace graphs that study how many reorderings are needed to update the indexes with respect to the length of the texts. For this purpose, we compute L ave for each text length and we count the number of reorderings after the insertion of a single letter for each chosen suffix length. A random letter is inserted at random positions 100,000 times and the average number of reorderings is computed over those insertions.

We choose to plot L ave and the actual number of reorderings for the most repeated prokaryote genome, M. flagellatus, whose sequence contains one huge repeated region. We focus on this specific genome that is a representative of all the other prokaryote genomes we mentioned. They all contain several large repeated regions, their curves are similar to Fig. 5. We also consider the two largest chromosome sequences from M. musculus as well as the two viruses that have been reported in Table 5. In both cases sizes are similar but sequences have totally different LCP distribution.

We remark that, since we consider suffixes, when a repetition appears at the beginning of the text, a peak will appear at the end of the graph. The M. flagellatus genome sequence, that possesses a single large repeat in the first Mbp, illustrates what has been described before (Fig. 5, whole sequence).

Moreover, when a repeat is close to the end, the number of LCP values that have to be considered is small, and therefore the L ave value will be larger compared to the same repeat closer to the start. Hence, two peaks of the same height, one at the beginning another at the end of a graph, represent two different repeat lengths. When necessary, we also provide the graph for the reversed text so that we can appreciate the bias due to a repeat at the end of a text (Fig. 5,M. musculus, chr. 2, reversed). The curve we obtain shows the same behavior for the peak while it reveals a slightly different aspect for the rest of the values.

Due to the length of chr. 1 from M. musculus, the respective values for L max and L ave (see Table 4), we are not expecting steep slopes or noticeable peaks. The curve that has been computed confirms our expectations. Contrarily, for chr. 2, L max is twenty times larger than L max for chr. 1 while L ave of chr. 2 is only six times larger than L ave of chr. 1. We are expecting one or several significant peaks. They are confirmed by the curves: it appears that one huge peak (corresponding to a large repeated region located at the end of the sequence) is masking the remaining LCP values.

Natural Language Texts

Similarly to the study we conducted for genome sequences, we now focus on both LCP values and curves analyses.

LCP values analysis

For etexts, we have very small LCP values as shown in Table 6. L perc is as low as 32 meaning that updating the suffix array is done by reordering at most 32 elements in 99% of the cases.

On the contrary, we have much larger LCP values with Wikipedia corpora. It was also expected, mainly because of inherent duplicates as explained in subsection 3.2. Moreover, removing a single duplicated template in the Afrikaans corpus makes L max , L ave or R FW values drop drastically. The L ave and I r values are significantly smaller for etexts than for Wikipedia texts. The redundant templates inserted in the latter being absent from the former, it reduces automatically these values. Nevertheless, even for the most repeated text, Afrikaans, we observe that less than 0.000 1% elements should be reordered on average and at most 0.05% in the worst case.

Curves analysis

For etexts we are expecting the L ave to quickly grow until it stabilizes. Doubling the size of the text will not have a significant impact on the L ave value. The curve clearly follows [START_REF] Fayolle | Analysis of the average depth in a suffix tree under a Markov model[END_REF], the R FW value indicating that it is close to their theoretical result.

In Fig. 6, the graph for etexts also confirms that we have very short repeats in natural language texts: we do not observe any significant peak. On the other hand, Afrikaans corpus in Wikipedia possesses some large repeats that can be observed with the rapid growth of L ave around 60 MB. In Fig. 6 (1), original Afrikaans is plotted and a peak clearly appears on the right. It corresponds to one single repeat of length 34,205 within a template of length 245,724. We then removed that particular template from Afrikaans, suppressing the corresponding peak in Fig. 6 (2). Similarly to what we did for M. musculus, we reversed the cleaned Afrikaans in order to attenuate the leftmost peak. We observe that Fig. 6 (3) is rather close to the etexts curve and that globally these two curves are following [START_REF] Fayolle | Analysis of the average depth in a suffix tree under a Markov model[END_REF].

We showed that L ave is an upper bound of the number of elements to reorder on average. We can observe that the actual number of reorderings is very close to the L ave for the texts we studied. It confirms that computing L ave gives some precious information about the expected number of reorderings that might be performed. It therefore permits to select a dynamic method or another strategy depending on L ave .

Studying L ave over a collection of texts

In the previous examples, we saw that L ave is always small with respect to n. In order to confirm this assertion, we consider the L ave of a larger set of texts (bacteria, eukaryotes, viruses, etexts and Wikipedia corpora). More precisely for one given type we have a specific graph where each plot corresponds to one text, its length on the x-axis and its L ave on the y-axis. We also plot a logarithmic function computed using non-linear regression, denoted by NLR, that fits best to the data. Results for genome sequences are given in Fig. 7 (bacteria), 8 (eukaryotes), 9 (viruses) and for natural language texts in Fig. 10 (etexts) and 11 (wikipedia).

Apart from few outliers, we observe that in all cases L ave values are at most about 100, and mainly concentrated around NLR. This tends to show that L ave is logarithmic not only for random texts, as shown by Fayolle and Ward [START_REF] Fayolle | Analysis of the average depth in a suffix tree under a Markov model[END_REF], among others, but also for more specific texts such as genome sequences and natural language texts. Now, let us focus on each figure and on the few outliers that can be observed. In Fig. 7, a dozen of bacteria have a L ave larger than 100. Their genome sequences contain very large duplicated regions that create a bias in the computation of the L ave (as presented in Table 1).

In Fig. 8, a majority (62 out of 78) of chromosome sequences have a L ave value below 50. However there is clearly a little number of chromosomes that are not so close to NLR. Indeed the eukaryotic sequences are well known for containing numerous repeats, as an example, more than 80% of the human genome sequence is made of repeats. Moreover there exists a high variability between chromosomes, as shown in Table 2. The maximal L ave value in Table 3 is smaller than 100 while the L ave values in Table 4 are larger, these values are still very modest with respect to the length of the chromosomes. The largest values can be observed for sexual chromosomes in both cases. They are due to the recombination rate that usually affects these chromosomes, mainly due to evolution. In Fig. 9, apart from eight virus sequences whose L ave values are above 30, all the other values are mainly located between 5 and 20. In Fig. 10, only 31 texts have a L ave above 15. These are all relatively small texts and can be classified in two categories. The first one consists in very small texts (around a thousand letters) which contain lyrics of English traditional songs. For instance one of them is "Heave Away", its refrain is 110 letters long and is repeated twice in this song which is 1,020 letters long. This repetition is not very large but it represents an important ratio of the text (about 20%), which explains that L ave is 32.74. The second category consists in larger texts (about a hundred thousand letters) where the content is centered using spaces, leading to large repetition of spaces.

In Fig. 11 limited number of texts have a quite large L ave value. Namely, 26 texts out of 151 have a L ave value larger than 200. We already explained why Wikipedia corpora contain more repetitions. This is mainly due to "copy and paste" of article content and template duplication. However these duplications appear in large corpora and therefore L ave is generally not much affected. Even in the worst case, L ave is 8,238 for a corpus of size 6,826,160. That means L ave represents only 0.12% of the total size.

Conclusions

We studied in this article the average number of reorderings when updating a suffix array-related structure. We showed the correlation between theory and practice, in Lemma 1 and showed in Lemma 2 that L ave is an upper-bound for the number of elements to reorder on average. We conducted experiments on very different types of texts: several DNA sequences of various repetitiveness and natural language texts. These experiments confirmed the correlation between the number of reorderings and L ave .

Fig. 1 .

 1 Fig. 1. BWT (CTAGTTAG$)=GTT$AATCG and corresponding SA, ISA, LCP.

  0 and the next values in the PLCP array decrease one by one until a value ℓ 1 at position j 1 . Then we have PLCP[j 1 + 1] = L 2 , with L 2 ≥ ℓ 1 , and again the values decrease, one by one, until a value ℓ 2 ≤ L 2 . And so forth until PLCP[j ′ ] = 1.

Fig. 5 .

 5 Fig. 5. Theoretical and actual number of reorderings for inserting a single letter in suffixes of DNA sequences.

  Afrikaans cleaned, reversed (3)Average case upper bound (L ave ) Practice

Fig. 6 .

 6 Fig.6. Actual and theoretical number of reorderings for inserting a single letter in suffixes of natural language texts.

3 NLRFig. 7 .

 37 Fig. 7. L ave for bacteria genome sequences

Fig. 8 . 3

 83 Fig. 8. L ave for eukaryote chromosome sequences (A. thaliana, C. elegans, D. melanogaster, H. sapiens, M. musculus, S. cerevisiae)

Fig. 9 .

 9 Fig. 9. L ave for virus genome sequences

  plots are more scattered than in the other figures. However still a

Fig. 10 . 3

 103 Fig. 10. L ave for etexts from Gutenberg project

Fig. 11 .

 11 Fig. 11. L ave for Wikipedia corpora

Table 1

 1 Values for the five most repeated prokaryote genomes.

	Organism Seq. Len.	L max	L ave	L perc	R n	R FW	I r
	M. flagellatus 2,971,519 143,034 3,452 113,320 1.16 • 10 -3 315.49 6.34
	S. agalactiae 2,211,485	47,068	546	24,954 2.47 • 10 -4	50.02 4.84
	D. ethenogenes 1,469,720	21,106	377	11,918 2.56 • 10 -4	36.41 4.03
	F. tularensis 1,892,775	33,912	336	14,984 1.78 • 10 -4	30.49 3.96
	N. meningitidis 2,272,360	32,158	261	9,435 1.15 • 10 -4	24.22 3.84

Table 2 (

 2 A. thaliana, C. elegans, D. melanogaster, S. cerevisiae), Table3(several chromosomes fromWe still observe some large values for L max but L ave and L perc are smaller, reflecting a different distribution of the LCP values. In this case, we have a very limited number of large repeats: as an example for H. sapiens chr. 1, L perc that corresponds to the percentile LCP value is equal to 171 while L max is equal to 67,631. All values between 171 and 67,631 are representing only 1% of all the LCP values. The remainder mostly consists in small repeats such as Short Tandem Repeats, micro and minisatellites or other well known biological structures. Once again, the large R FW values for several chromosomes confirm that a Markovian model of order one is not suitable for every DNA sequences. However, for some other chromosomes, R FW values are close to 1, indicating that despite the low Markovian model, the theoretical value can be a pretty good approximation for less repetitive texts.

H. sapiens) and Table

4

(several chromosomes from M. musculus), generally contain a smaller percentage of repeats than prokaryote genomes but they are representative among the different species. Hence, on average, LCP values are much smaller than for the five most repeated prokaryote genome sequences.

Table 6

 6 Values for various natural language texts. etexts is a concatenation of several texts from Gutenberg project. The other texts correspond to the Wikipedia corpus in the given language.

	Text	Length	L max L ave L perc	R n R FW	I r
	etexts	91,070,340	1,727	12	32 1.32 • 10 -7	1.7 0.82
	Afrikaans	68,989,658 34,205	66	383 9.57 • 10 -7	10.1 3.00
	Af. (cleaned)	68,743,934	6,831	31	295	4.5 • 10 -7	4.7 1.98
	Basque 139,868,091 13,879	45	406 3.22 • 10 -7	6.3 2.16
	Bosnian 122,215,463 18,314	40	318 3.27 • 10 -7	5.9 2.28
	Estonian 213,877,916 26,119	39	284 1.82 • 10 -7	5.4 2.16
	Latin	85,735,379	8,099	38	346 4.43 • 10 -7	5.6 2.15
	Occitan	70,250,160 11,003	64	914 9.11 • 10 -7	9.4 2.64

http://ncbi.nlm.nih.gov/pmc/articles/PMC1769404/bin/1471-2105-7-541-S1.pdf

Downloaded from http://www.gutenberg.org

Downloaded on March 30 th 2009 at http://download.wikimedia.org

The Wikipedia languages have been selected because of the length of the corpora which can fit within our computer main memory, each of them representing more than 70,000,000 bytes.

http://web.unipmn.it/∼manzini/lightweight/

We furthermore studied the distribution of L ave when the size of the text grows and we have seen that it follows a logarithmic function not only for random texts but also for the texts we mentioned above.

We conclude that the number of elements to be reordered when updating a suffix array-related structure is small and generally logarithmic even on DNA sequences or natural language texts. It is therefore reasonable to update the index rather than building a new one when insertion, deletion or substitution operations occur with the frequencies we measured for both natural language texts and genomic sequences (eg. see Fig. 6 in [21]).


Yasmina Abdeddaïm
email: y.abdeddaim@esiee.fr

Damien Masson
email: d.masson@esiee.fr

Scheduling Self-Suspending Periodic Real-Time Tasks Using Model Checking

In this paper, we address the problem of scheduling periodic, possibly self-suspending, real-time tasks. We provide schedulability tests for PFP and EDF and a feasibility test using model checking. This is done both with and without the restriction of work-conserving schedules. We also provide a method to test the sustainability w.r.t the execution and suspension durations of the schedulability and feasibility properties within a restricted interval. Finally we show how to generate an on-line scheduler for such systems when they are feasible.

I. INTRODUCTION

A real-time task can suspend itself when it has to communicate, to synchronize or to perform external input/output operations. Classical models neglect self-suspensions, considering them as a part of the task computation time [START_REF] Liu | Real-Time Systems[END_REF]. Models that explicitly consider suspension durations exist but their analysis is proved difficult: in [START_REF] Ridouard | Negative results for scheduling independent hard real-time tasks with self-suspensions[END_REF], the authors present three negative results on systems composed by hard real-time selfsuspending periodic tasks scheduled on-line. We propose in this paper the use of model checking on timed automata to address these three negative results: 1) the scheduling problem for self-suspending tasks is NP-hard in the strong sense, 2) classical algorithms do not maximize tasks completed by their deadlines and 3) scheduling anomalies can occur at run-time. Result 1) means that there cannot exist a non-clairvoyant online algorithm that takes its decisions in a polynomial time and always successfully schedules a feasible self-suspending task set. We so propose to use model checking to generate off-line a feasible scheduler for each specific instances of the problem, i.e. for each task sets. Result 2) implies that traditional on-line schedulers are not optimal, whereas our approach is. Result 3) points out that changing the properties of a feasible task set in a positive way (e.g. reducing an execution or a suspension duration, extending a period) can affect its feasibility. We prove that our approach is sustainable w.r.t execution and suspension durations. Related work: In [START_REF] Lakshmanan | Scheduling self-suspending real-time tasks with rate-monotonic priorities[END_REF], the authors prove that the critical scheduling instant characterization is easier in the context of sporadic real-time tasks. They provide, for systems scheduled under a rate-monotonic priority assignment rule, a pseudopolynomial response-time test. The rest of the literature on self-suspending tasks focus on the multiprocessor context [START_REF] Liu | Task scheduling with self-suspensions in soft real-time multiprocessor systems[END_REF], [START_REF]Improving the schedulability of sporadic self-suspending soft real-time multiprocessor task systems[END_REF]. Our approach addresses the problem for periodic tasks, and is not restricted to RM. The timed automata approach has been already used to solve job shop scheduling problems [START_REF] Abdeddaïm | On optimal scheduling under uncertainty[END_REF], [START_REF] Fehnker | Scheduling a steel plant with timed automata[END_REF]. In [START_REF] Fersman | Task automata: Schedulability, decidability and undecidability[END_REF], the authors present a model based on timed automata to solve real-time scheduling problems, but this model cannot be applied to self-suspending tasks because of the way the preemptions are modeled. Moreover, it intrinsically considers work-conserving schedules and cannot handle uncertain times. Contributions: In this paper, we propose a timed-automatabased model for hard real-time periodic self-suspending tasks. We show how to use model checking to obtain both a necessary and sufficient feasibility test, and a schedulability test for classical scheduling policy (RM, DM, EDF). When these algorithms fail to schedule a feasible system, we show how to generate an appropriate scheduler. We also consider the cases where both the execution time and the suspension time of each task are uncertain. The generated schedulers then have an important property: the feasibility of a task set is sustainable w.r.t the execution and suspension durations.

Sect. II presents the task model, Sect. III introduces the selfsuspending task automaton, Sect. IV exposes how to check the feasibility and the schedulability with PFP and EDF, and how to generate a sustainable scheduler. Sect. V presents experiments and finally we conclude.

II. SELF-SUSPENDING TASK MODEL

Task Model and assumptions: We consider the problem of scheduling a system Σ = {τ 1 , ..., τ n } of n independent possibly self-suspending periodic tasks synchronously activated on one processor. A task which does not suspend itself is characterized by the tuple τ i = (C i , T i , D i ) where C i is its execution time, T i its period and D i its relative deadline with D i ≤ T i . A task which suspends itself is characterized by the tuple τ i = (P i , T i , D i ) where P i is its execution pattern. An execution pattern is a tuple

P i = (C 1 i , E 1 i , C 2 i , E 2 i , ..., C m i )
where each C j i is an execution time and each E j i a suspension time. When these parameters can take values within an interval, we refer to the tasks as uncertain tasks and the execution pattern becomes

P i = ([C 1 i l , C 1 i u ], [E 1 i l , E 1 i u ], [C 2 i l , C 2 i u ], ..., [C m i l , C m i u ]
). To simplify the notations, when there is no ambiguity, the task id is not mentioned and the tuples becomes τ = (C, T, D) etc. This model and notations are inspired by existing literature on self suspending tasks [START_REF] Liu | Task scheduling with self-suspensions in soft real-time multiprocessor systems[END_REF], [START_REF] Liu | Real-Time Systems[END_REF], [START_REF] Lakshmanan | Scheduling self-suspending real-time tasks with rate-monotonic priorities[END_REF], [START_REF] Ridouard | Negative results for scheduling independent hard real-time tasks with self-suspensions[END_REF]. Sustainability: The schedulability of a task set with a given algorithm is said sustainable w.r.t. a parameter when a schedulable task set remains schedulable when this parameter is changed in a positive way. The sustainability is an important property since it permits to study the worst case scenario. In our task model, execution and suspension durations can be bounded within intervals. In the remaining of the paper, we

i = m cτ = C -1 ∧ xτ = 1 xτ ≤ 1 cτ = C -1 ∧ xτ = 1 ∧ i < m xτ = E e p xτ ≤ E f xτ := 0 s sp id C := C 0 ∧ i := 0 E := E i ∧ xτ := 0 i := i + 1 ∧ C := C i ∧ cτ := 0 xτ := 0 ∧ cτ := 0 ∧ z := 1 cτ < C ∧ xτ = 0 ∧ z = 0 xτ := 0 ∧ z := 1 cτ := cτ + 1 ∧ xτ := 0 ∧ z := 0 cτ < C -1 ∧ xτ = 1 Figure 1.
Self-Suspending Task. UPPAAL model is available in [START_REF] Abdeddaïm | Uppaal and uppaal-tiga implementations[END_REF] say that the schedulability is sustainable when the task set is feasible with all the possible values in the intervals. By extension we say that an algorithm is sustainable when the schedulability of a task set with this algorithm is sustainable.

III. THE SELF-SUSPENDING TIMED AUTOMATON MODEL

A timed automaton [START_REF] Alur | Automata for modeling real-time systems[END_REF], is an automaton augmented with a set of real variables X , called clocks. The firing of a transition can be controled by a clock constraint called guard and a location can be constrained by a staying condition called invariant. A clock valuation is a function that associates to a clock x its value v(x) and a configuration is a pair (q, v) where q is a state of the automaton and v a vector of clock valuations. A run is a sequence of timed and discrete transitions, timed transitions representing the elapse of time in a state, and discrete ones the transitions between states. Synchronous communication between timed automata can be done using input actions (a?) and output actions (a!). Self-Suspending Timed Automaton: Let τ = (P, T, D) be a self-suspending task. We associate to τ a timed automaton A τ (Fig. 1) with one clock x τ and of states Q = {s, e, f, p, sp, id}. State s is the waiting state where the task is active, e the execution state, p the state where the task is preempted, f is the one where the task has finished its execution, sp the suspension state and id is an idle state used to model a possible idle step. To compute the execution time of a task, we use the clock x τ and a variable c τ as follows: the automaton can stay in the execution state exactly one time unit and the variable c τ keeps track of how many time units have been performed. This is represented by an invariant x τ ≤ 1 on state e and a loop transition that increments c τ . The guard c τ < C -1 restricts the number of loop transitions. The preemption of the task is modeled using transitions from e to p and from p to e. Note that the modeling of preemption at any time is not possible using timed automata because clock variables cannot be stopped. Thus, it is supposed in this model that a task can be preempted only at integer times. The transition from e to id is enabled when the total time spent in the active state is equal to C and the number of performed step is less than m. To model the suspension of a task we use the state sp. The automaton can stay in the state sp exactly E time units. This is modeled using an invariant x τ ≤ E on the state and a guard x τ = E. The task finishes when all the steps have been computed. This is modeled using the guard i = m from state e to state f . To make this task periodic of period T , we use a second timed automaton A T with one state and a loop transition enabled every T time units. This transition is labeled with an output action T ! and synchronizes with the transition from waiting state (f ) to active state (s) of the automaton A τ . Finally, we introduce a third automaton A D , which models the deadline of the task using an output action D! fired every D time units. We call the tuple (A τ , A T , A D ) a A self-suspending automata model.

IV. SCHEDULABILITY USING MODEL CHECKING

Model checking is an automatic verification technique used to prove formally whether a model satisfies a property or not. In this section, we present how we use CTL [START_REF] Kozen | Logics of Programs[END_REF] model checking to test the feasibility of a task set, its schedulability with PFP and EDF, and the sustainability of a schedule.

A. Feasibility and Schedulability

Let Σ = {τ 1 . . . τ n } be a finite set of self-suspending tasks. We associate to every task τ i a self-suspending automata model

(A i τ , A i T , A i D )
. We use a global variable proc to ensure, using a guard, that only one task is executing at once. We introduce a new state ST OP i for each automaton A i τ , which is reached if a task misses its deadline. For every non final state of A i τ , a transition labeled by an input action D i ? leads to state ST OP i . These transitions synchronize with the deadline automaton. Thus, if an instance of a task does not terminate before its deadline, the automaton goes to the state ST OP i .

Proposition 1 (feasibility): Let Σ = {τ 1 . . . τ n } be a set of self-suspending tasks. Σ is feasible iff CTL Formula 1 is satisfied.

φ Sched : EG¬(

i∈[1,n] ST OP i ) 1 (1)
Sketch of Proof: Proposition 1 states that the selfsuspending problem is feasible iff there exists a feasible scheduling run i.e a run that never reached a STOP state. Suppose that the scheduling problem Σ is feasible and Formula 1 is not satisfied. If the problem is feasible, then there exists a schedule where all the instances of all tasks never miss their deadline. This schedule corresponds in the self-suspending automaton to a feasible scheduling run. This contradicts the hypothesis that Proposition 1 is not satisfied. Suppose now that Formula 1 is satisfied and the scheduling problem is not feasible. If the formula is not satisfied, then all scheduling runs are not feasible i.e all the scheduling run lead to a STOP state. This contradicts the fact that the problem is feasible, the contradiction comes from the fact that the automaton capture all possible behaviors of task instances.

An on-line scheduling algorithm can be obtained using a feasible run satisfying Formula 1 by simply reading sequentially the configurations of the feasible run until reaching the one where all active tasks have terminated their execution without missing their deadline.

Proposition 2 (PFP Schedulability): Let Σ be a set of selfsuspending tasks sorted according to a priority function, with τ 1 ≤ . . . ≤ τ n . Σ is feasible according to a fixed priority scheduler iff CTL Formula 2 is satisfied.

Cu := C 0 u ∧ Cl := C 0 l ∧ i := 0 cτ ≥ Cl -1 ∧ xτ = 1 ∧ i < m cτ ≥ Cl -1 ∧ xτ = 1 cτ < Cu ∧ xτ = 0 ∧ z = 0 xτ := 0 ∧ cτ := 0 ∧ z := 1 xτ ≤ 1 xτ := 0 ∧ z := 1 xτ ≥ El Eu := E i u ∧ El := E i l ∧ xτ := 0 xτ ≤ Eu i := i + 1 ∧ Cu := C i u ∧ Cl := C i l ∧ cτ := 0 cτ := cτ + 1 ∧ xτ := 0 ∧ z := 0 cτ < Cu -1 ∧ xτ = 1
Figure 2. Uncertain Self-Suspending Task. UPPAAL-TIGA model is available at [START_REF] Abdeddaïm | Uppaal and uppaal-tiga implementations[END_REF]. Uncontrollable transitions are represented using dashed lines.

φ F P :EG¬( i∈[1,n-1] (s i j∈[i,n] e j ) ∨ i∈[1,n-1] (p i j∈[i,n] e j )∨ ( i∈[1,n] ST OP i )) (2)
Formula 2 states that there exists a feasible run where, in all the configurations, a task cannot be in its execution state if a less priority task is active i.e. in state s or p. The idea of the proof is similar to the one of Proposition 1.

Proposition 3 (EDF Schedulability): Let Σ be a set of selfsuspending tasks. Σ is schedulable according to EDF iff CTL Formula 3 is satisfied.

φ edf :EG¬( i∈[1,n] j =i∈[1,n] (s i ∧ e j ∧ P ij ) i∈[1,n] j =i∈[1,n] (p i ∧ e j ∧ P ij ) ∨ ( i∈[1,n] ST OP i )) (3)
P ij is a state of an observer automaton reachable when x Di -x Dj > D i -D j with x Di and x Dj the clocks of the deadline automata A i D and A j D . Formula 3 states that there exists a run where, in all the configurations, a task cannot be in its execution state if a task with a closer deadline is active. The idea of the proof is similar to the one of Proposition 1.

B. Sustainable Scheduler

In a timed game automaton (TGA) [START_REF] Maler | On the synthesis of discrete controllers for timed systems[END_REF], the set of transitions is split into controllable (∆ c ) and uncontrollable (∆ u ) ones. Solving a timed game consists in finding a strategy f s.t. a TGA supervised by f always satisfies a given formula. Fig. 2 represents our model adapted to uncertain self-suspending tasks. The start and preemption transitions are controlled by the scheduler, while the transitions from state e to f , from e to id and from sp to p are controlled by the environment. The guard c τ ≥ C l -1 on transitions e to id and e to f models that the task can terminate after C l time units and the guard c τ < C u -1 models that the task terminates before C u time units. The invariant x τ ≤ C u and the guard x τ ≥ E l from sp to p models a suspension of a duration within [E l , E u ].

Definition 1 (Scheduling strategy): Let A be a set of n uncertain self-suspending task automata. A scheduling strategy is a function f from the set of configurations of A to the set ∆ 1 c ∪ . . . ∆ n c ∪ {λ}. s.t. if f ((q, v)) = e ∈ ∆ c then execute the controllable transition e and if f ((q, v)) = λ then wait in the configuration (q, v).

Proposition 4 (Sustainability): Let Σ = {τ 1 . . . τ n } be a finite set of uncertain self-suspending tasks. A sustainable scheduling algorithm exists for Σ iff there exists a scheduling strategy f s.t the self-suspending timed game automata model of Σ supervised by f satisfies the safety CTL Formula 4.

φ sust : AG¬( i∈[1,n] ST OP i ) (4)
Sketch of Proof: Formula 4 states that there exists a strategy function f s.t. for every configuration of the task set and every possible execution duration or suspension duration, there is a way to avoid the STOP states. Let Σ be a set of uncertain self-suspending tasks. Suppose that there exists a sustainable feasible scheduler for Σ. Then, the strategy f can be simply computed by defining a strategy that mimic the decisions of the sustainable scheduler. Consider now that there exists a feasible scheduling strategy for the timed game satisfying Formula 4, then there exists a way for all the tasks to not miss their deadline whatever are the execution and the suspension durations. Thus this strategy can be used as a sustainable scheduling algorithm.

Algorithm 1 Scheduling Strategy Algorithm 1: (q, v) ← (q 0 , v 0 ), t ← 0 2: while q = q 0 or t = 0 do 3:

while f ((q, v)) = λ or no task finished or no end of suspension do

4:

Wait: increase v

5:

end while 6:

if f ((q, v)) = tr ∈ ∆ j c then 7:
(q k , v k ) is the successor of (q, v) while taking the transition tr

8:

if ∃q j k = q j and q j k = e j then 9:

execute task τ j at time t ← t k 10:

end if

11:

if ∃q j k = q j and q j k = p j then 12:

preempt task τ j at time t ← t k

13:

end if

14:

if ∃q j k = q j and q j k = sp j then 15:

suspend task τ j at time t ← t k

16:

end if

17:

end if

18:

if a task τ j has finished then 19:

(q k , v k ) is the configuration (q, v) where q j k ← f j , v k (x j τ ) ← 0 20: end if

21:

if a task τ j has terminate the suspension then

22:

(q k , v k ) is the configuration (q, v) where q j k ← p j , v k (x j τ ) ← 0

23:

end if

24:

(q, v) ← (q k , v k )

25: end while

To provide a sustainable scheduling algorithm we need to construct a feasible strategy if one exists. Such a strategy is finite, because of the decidability of the timed game problem [START_REF] Henzinger | Discrete-time control for rectangular hybrid automata[END_REF]. Then an on-line scheduler is an algorithm that executes the pre-computed strategy. This is formalized by Algorithm 1 2 . According to the actual configuration, the scheduling algorithm can decide 1) to stay in this configuration, i.e to continue the execution of a task or let the processor idle (lines 3-5) ; or 2) to execute, preempt or suspend a task (lines 6-17). Finally when an execution or a suspension terminates, the algorithm computes the new configuration (lines 18-25).

Note that if the presented methods generate possibly idle schedules, we can compute work-conserving schedules by specifying this property in CTL Formulas 1, 2, 3 and 4.

Figure 3 .

 3 Figure 3. Feasible schedule exists but neither or EDF is able to find it

The operator A means for all paths, E there exists a path and G globally in the future[START_REF] Kozen | Logics of Programs[END_REF].

i is the state of the automaton j in the configuration (q i , v i ) ; q 0 is the initial configuration ; t is an additional global clock which is never reset ; t i the valuation of the clock t in the configuration (q i , v i ).

00 00 11 11 00 00 11 11 00 00 11 11 00 00 11 11 00 00 11 11 00 00 11 11 00 00 11 11 00 00 11 11 00 00 11 11 00 00 11 11 00 00 11 11 00 00 11 11 00 11 00 11 0 0 1 1 0 0 1 1 0 0 1 1 0 0 1 1 0 0 1 1 0 0 1 1 0 0 1 1 0 0 1 1 0 0 1 1 0 0 1 1 0 0 1 1 0 0 0 1 1 1 0 0 0 1 1 1 0 0 0 1 1 1 0 0 0 1 1 1 0 0 0 1 1 1 0 0 0 1 1 1 0 1 0 1 0 1 0 1 00 11 00 11 00 11 00 11 00 11

V. EXPERIMENTS

We used UPPAAL [START_REF] Larsen | Uppaal in a nutshell[END_REF] and UPPAAL-TIGA [START_REF] Behrmann | Uppaal-tiga: Time for playing games![END_REF] to implement our model [START_REF] Abdeddaïm | Uppaal and uppaal-tiga implementations[END_REF]. We tested it on two examples. The first is composed by two regular self-suspending tasks. The second is composed by three uncertain self-suspending tasks. Fig. 3 and Fig. 4 present the obtained results. White squares represent suspension durations and hatched ones execution durations. Example 1 (regular self-suspending tasks) In this experiment, we have modeled the system Σ = {τ 1 , τ 2 } with τ 1 = ((1, 4, 1), 7, 7) and τ 2 = ((1, 3, 1), [START_REF] Abdeddaïm | On optimal scheduling under uncertainty[END_REF][START_REF] Abdeddaïm | On optimal scheduling under uncertainty[END_REF]. We have first used Formula 2 with RM priority assignments. The property is not verified, this result permits us to conclude that the task set is not schedulable according to RM. The same result is obtained with the inversed priority assignment. Sub-Fig. 3(a) and 3(b) validate these results: we see that τ 2 effectively misses a deadline at time 6 with inverse RM, and that τ 1 misses a deadline at time 7 with RM. We have then used Formula 3 to test the feasibility with EDF. The property is not verified, this can be confirmed by Sub-Fig. 3(c) that shows that τ 2 misses a deadline at time 42 under EDF. Finally we have used Formula 1 to test the unconstrained feasibility. The property is verified, thus a feasible schedule exists for this task problem. Using the produced feasible scheduling run, we are effectively able to produce the schedule presented by Sub-Fig. 3(c) (TAAS stands for Timed-Automata-Assisted Scheduler). Example 2 (uncertain self-suspending tasks) In this experiment, we have modeled the system Σ = {τ 1 , τ 2 , τ 3 } with τ 1 = ((2, 2, 4), 10, 10), τ 2 = ((2, 8, 2), 20, 20) and τ 3 = ((2), [START_REF] Abdeddaïm | Uppaal and uppaal-tiga implementations[END_REF][START_REF] Abdeddaïm | Uppaal and uppaal-tiga implementations[END_REF], where τ 1 has the highest priority and τ 3 the lowest. We first have verified Formula 2: the property is verified, the system is then feasible with a fixed priority scheduler. We then have modeled the system [START_REF] Liu | Real-Time Systems[END_REF][START_REF] Ridouard | Negative results for scheduling independent hard real-time tasks with self-suspensions[END_REF], 4), 10, 10). We have verified Formula 4 restricted to fixed priority schedulers on our model, the property is not verified. We conclude that feasibility with a fixed priority scheduler is not sustainable for this system. Subfigure 4(a) presents the schedule obtained with Σ. Sub-figure 4(b) presents the schedule with Σ * where the third instance of τ 1 executes with the pattern P 1 = (

). It results in a deadline missed for τ 3 at time 49. However, we have tested Formula 4. The outcome is positive in both cases: valid schedules restricted and non restricted to work-conserving ones. The feasibility of the system is then sustainable (within the intervals

) in the general case and with a work-conserving scheduler. Indeed, there exists a simple way to enforce the sustainability: forcing the system to insert idle times when a task completes earlier than it was supposed to. Fig. 4(c) shows the resulting schedule of this strategy. Fig. 4(d) presents a work-conserving feasible schedule which can be obtained using the strategy generated by UPPAAL-TIGA.

VI. CONCLUSION

In this paper, we present how to use model checking to solve a difficult scheduling problem: the scheduling of periodic self-suspending tasks. We provide a feasibility test and schedulability tests with PFP and EDF. We also provide a method to test the sustainability of schedules w.r.t the execution and suspension durations. This is done both with the restriction of work-conserving schedules and in the general case. Finally our approach permits to generate a scheduler for such systems.

Work in progress:

We first have to implement the scheduler generation and to formalize the memory complexity of generated on-line schedulers. Then we have to improve the way our automata handle preemptions. Finally we have to extend our model to consider multiprocessor platforms and tasks sporadic activation and compare the results with the ones presented in [START_REF] Lakshmanan | Scheduling self-suspending real-time tasks with rate-monotonic priorities[END_REF].