Mathieu Giraud
email: mathieu.giraud@lifl.fr

Asymptotic behavior of the numbers of runs and microruns ✩

Keywords: combinatorics, maximal repetitions, runs, asymptotic behavior, maximum number of runs

. However, very few properties are known for the ρ(n)/n function. We show here by a simple argument that lim n →∞ ρ(n)/n exists and that this limit is never reached. We further study the asymptotic behavior of ρ p (n), the maximal number of runs with period at most p. Finally, we provide the first exact limits for some microruns. For example, we have lim n →∞ ρ 14 (n)/n = 15/17.

Introduction

The study of repetitions is an important field of research, both for word combinatorics theory and for practice, with applications in domains like computational biology or cryptanalysis. The notion of run (also called maximal repetition or m-repetition [START_REF] Kolpakov | Maximal repetitions in words or how to find all squares in linear time[END_REF]) allows a compact representation of the set of all tandem periodicities, even fractional, in a string. The proper counting of those runs is important for all algorithms dealing with repetitions.

Since the work of Kolpakov and Kucherov in [START_REF] Kolpakov | Maximal repetitions in words or how to find all squares in linear time[END_REF][START_REF] Kolpakov | On maximal repetitions in words[END_REF], it is known that ρ(n), the maximum number of runs in a string, is linear in the length n of the string. They gave the first algorithm computing all runs in a linear time, but without an actual constant. 1: Values of ρ(n) for binary strings, from [START_REF] Kolpakov | Maximal repetitions in words or how to find all squares in linear time[END_REF].

Upper bounds have been recently provided by Rytter (5n) [START_REF] Rytter | The number of runs in a string : improved analysis of the linear upper bound[END_REF] and Puglisi, Simpson, and Smyth (3.48n) [START_REF] Puglisi | How many runs can a string contain?[END_REF]. A 1.6n bound was obtained by Crochemore and Ilie [START_REF] Crochemore | Maximal repetitions in strings[END_REF]. They count separately the microruns, that is the runs with short periods, and the runs with larger ones. They show that the number of microruns with period at most 9 satisfies ρ 9 (n) ≤ n. For larger runs, they prove that

ρ ≥p (n) ≤ 2 p ∞ i=0 2 3 i n = 6 p • n.
Crochemore, Ilie, and Tinta extended those results with massive computations, bringing down the upper bound from 1.6n to 1.048n [START_REF] Crochemore | Towards a solution to the runs conjecture[END_REF].

A lower bound of αn, with α = 3/(1 + √ 5) = 0.927..., has been given by [START_REF] Frantisek Franek | The maximum number of runs in a string[END_REF] then [START_REF] Franek | An asymptotic lower bound for the maximal-number-of-runs function[END_REF]. In [START_REF] Frantisek Franek | The maximum number of runs in a string[END_REF], Franek, Simpson, and Smyth propose a sequence of strings (x n ) with increasing lengths such that lim n →∞ r(x n )/|x n | = α, where r(x) is the number of runs in the string x. In [START_REF] Franek | An asymptotic lower bound for the maximal-number-of-runs function[END_REF], Franek and Yang show that α is an asymptotic lower bound by showing that there exists a whole family of asymptotic lower bounds arbitrarily close to α. Recently, Matsubara et al. provided an 174719/184973 = 0.9445.. lower bound by repeating a large run-rich string [START_REF] Matsubara | New lower bounds for the maximum number of runs in a string[END_REF].

In fact, very few properties are known for the ρ(n)/n function [START_REF] Franek | An asymptotic lower bound for the maximal-number-of-runs function[END_REF][START_REF] Smyth | The maximum number of runs in a string[END_REF]. In this paper, after giving some definitions (Section 2), we show by a simple rewriting argument that ℓ = lim n →∞ ρ(n)/n exists and that this limit is never reached (Section 3.1), proving that

ρ(n) n ≤ ℓ - 1 4n .
Section 3.2 proves the convergence of ρ(n)/n even in the case of a fixed alphabet, for example for binary strings. In Section 4.1, we further study the asymptotic behavior of ρ p (n), the number of runs with short periods, showing that ℓ p = lim n →∞ ρ p (n)/n exists and that, for some constant z p ,

ℓ p - z p n ≤ ρ p (n) n ≤ ℓ p ≤ ℓ.
Moreover, we provide a simple way to exactly count some microruns (Section 4.2). We give in Table 3 the first exact limits ℓ p for microruns on binary strings with p ≤ 14. Section 5 gives some concluding remarks.

Definitions

Let x = x 1 x 2 . . . x n be a string over an alphabet. Let p ≥ 1 be an integer. We say that x has a period p if for any i with 1 ≤ i ≤ n -p, x i+p = x i . We denote by x[i . . j] the substring x i x i+1 . . . x j . A run is an interval [i . . j]:

• such that x[i . . j] has period p ≤ (j -i + 1)/2,
• that is maximal: if they exist, neither

x i-1 = x i-1+p , nor x j+1 = x j+1-p ,
• and such that x[i . . i + p -1] is primitive: it is not an integer power of another string.

We define by r p (x) the number of runs of period ≤ p in x, called microruns in [START_REF] Crochemore | Maximal repetitions in strings[END_REF], and by r(x) = r ⌊|x|/2⌋ (x) the total number of runs in x. For example, the four runs of x = atattatt are [4 . . 5] (tt), [7 . . 8] (tt), [1 . . 4] (atat) and [2 . . 8] (tattatt), and thus r 1 (x) = 2, r 2 (x) = 3, and r 3 (x) = r(x) = 4. Given an integer n ≥ 2, we now consider all strings of length n. We define as

ρ p (n) = max{r p (x) | |x| = n}
the maximum number of runs of period ≤ p in a string of length n. Then we define as

ρ(n) = max{r(x) | |x| = n} = ρ ⌊n/2⌋ (n)
the maximum number of runs in a string of length n. Kolpakov and Kucherov gave in [START_REF] Kolpakov | On maximal repetitions in words[END_REF] some values for ρ(n) (Table 1). Table 4, at the end of this paper, shows some values for ρ p (n). Note that r(x) = ρ(|x|) does not imply that r p (x) = ρ p (|x|) for all p: for example, r(aatat) = 2 = ρ(5) but r 1 (aatat) = 1 < ρ 1 (5) = 2. Finally, we can define values r ≥p (x) and ρ ≥p (n) for macroruns, that is runs with a period at least p. Again, r(x) = ρ(|x|) does not imply that r ≥p (x) = ρ ≥p (|x|). For example, r ≥2 (aatt) = 0 < ρ ≥2 (4) = 1 = r ≥2 (atat).

On the number of runs

Rewritings and asymptotic behavior of the number of runs

Franek et al. [START_REF] Frantisek Franek | The maximum number of runs in a string[END_REF][START_REF] Franek | An asymptotic lower bound for the maximal-number-of-runs function[END_REF] list some known properties for ρ(n):

• For any n, ρ(n + 2) ≥ ρ(n) + 1 • For any n, ρ(n + 1) ≤ ρ(n) + ⌊n/2⌋
• For some n, ρ(n + 1) = ρ(n)

• For some n, ρ(n + 1) = ρ(n) + 2

We add the following two simple properties.

Proposition 1. The function ρ is superadditive: for any m and n, we have

ρ(m + n) ≥ ρ(m) + ρ(n).
Proof. Take two strings x and y of respective lengths m and n such that r(x) = ρ(m) and r(y) = ρ(n). Let ȳ be a rewriting of y with characters not present in x. (See below for a discussion on the size of the alphabet.) Then xȳ is a string of length m + n containing exactly the runs of x and the rewritten runs of y. Thus ρ(m

+ n) ≥ r(xȳ) = r(x) + r(y) = ρ(m) + ρ(n).
For any t ≥ 1, we have in particular ρ(tn) ≥ tρ(n).

Proposition 2. For any n, ρ(4n) ≥ 4ρ(n) + 1

Proof. Take a string x of length n with r(x) = ρ(n). Let x be a rewriting of x with characters not present in x. Then r(xxxx) ≥ 4r(x) + 1.

We can now state our main result:

Theorem 1. ρ(n)/n converges to its upper limit ℓ. Moreover, the limit is never reached, as for any n we have

ρ(n) n ≤ ℓ - 1 4n .
Proof. Let ℓ be the upper limit of ρ(n)/n. This limit is finite because of [START_REF] Kolpakov | On maximal repetitions in words[END_REF]. Given ε, there is a n 0 such that ρ(n 0 )/n 0 ≥ ℓ -ε/2. For any n ≥ n 0 , let be t = ⌊n/n 0 ⌋. Then we have ρ(n)/n ≥ ρ(tn 0 )/n ≥ tρ(n 0 )/n by Proposition 1, thus ρ(n)/n ≥ t/(t + 1) • ρ(n 0 )/n 0 . Let be t 0 such that t 0 /(t 0 + 1)

• ρ(n 0 )/n 0 ≥ ρ(n 0 )/n 0 -ε/2. Then, for any n ≥ t 0 n 0 , we have ρ(n)/n ≥ ℓ -ε, thus ℓ = lim n →∞ ρ(n)/n. Finally, Proposition 2 gives ℓ ≥ ρ(4n)/4n ≥ ρ(n)/n + 1 4n .
The proof of convergence of f (n)/n when f is superadditive is known as Fekete's Lemma [START_REF] Fekete | Über die Verteilung der Wurzeln bei gewissen algebraischen Gleichungen mit ganzzahligen Koeffizienten[END_REF][START_REF] Van Lint | A course in combinatorics[END_REF]. This convergence result was an open question of [START_REF] Franek | An asymptotic lower bound for the maximal-number-of-runs function[END_REF]. In fact, the motivation of [START_REF] Franek | An asymptotic lower bound for the maximal-number-of-runs function[END_REF] was the remark that "the sequence |x i | (of [START_REF] Frantisek Franek | The maximum number of runs in a string[END_REF]) is only "probing" the domain of the function ρ(n) and r(x i ) is "pushing" the value of ρ(n) above αn in these "probing" points". Then Franek and Yang [START_REF] Franek | An asymptotic lower bound for the maximal-number-of-runs function[END_REF] prove that every α -ε is an actual asymptotic lower bound by building specific sequences. With Propositions 1 and 2 and Theorem 1, it is now sufficient to study bounds on any (ρ(n i )/n i ) sequence (for a growing sequence (n i )) to give bounds on ρ(n)/n.

Note that this convergence does not imply monotonicity. In fact, if ℓ < 1, then ρ(n)/n is asymptotically non monotonic, as there will be in this case an infinity of n's such that ρ(n + 1) = ρ(n). Note also that, although Proposition 1 and 2 require to double the alphabet size, the alphabet remains finite: the proof of Theorem 1 only requires to double once this alphabet size. Moreover, it is possible to prove Proposition 1 without rewriting in a larger alphabet, thus proving the convergence of ρ(n)/n when considering only binary strings. This second proof, more elaborated, is given in the next section.

The bound ℓ -1 4n can be improved. For example, with a rewriting similar to the one used in Proposition 2, it can be shown that ρ(2n 2 ) ≥ (2n + 1)ρ(n), giving by successive iterations ρ(n)/n ≤ ℓ -1 2n . This has not been reported here to keep the proof simple.

Concerning microruns with period at most p, Proposition 1 still holds:

Proposition 3. For any p, m, and n, we have

ρ p (m + n) ≥ ρ p (m) + ρ p (n).
Thus for any p, ρ p (n)/n converges to its upper limit ℓ p .

The proof is the same as above. On the contrary, Proposition 2 may be not true for microruns. For example, for any n, ρ 1 (n) = ⌊n/2⌋, and thus for any even n, we have

ρ 1 (n)/n = ℓ 1 = 1/2.
Finally, Theorem 1 is fully valid for macroruns. Moreover, taking the following inequality to the limit, we have ℓ ≥p ≥ ℓ/p. Proposition 4. For any p and n, ρ ≥p (pn) ≥ ρ(n).

Proof. Take a string x = x 1 x 2 . . . x n of length n with r(x) = ρ(n). Let x be the string x p 1 x p 2 . . . x p n of length pn. Then r ≥p (x) ≥ r(x).

A proof of Proposition 1 for fixed alphabets

Here we prove Proposition 1 without rewriting in a larger alphabet, thus proving the convergence of ρ(n)/n when considering only binary strings. This proof is borrowed and simplified from one part of a proof of Franek et al. (Theorem 2 of [START_REF] Frantisek Franek | The maximum number of runs in a string[END_REF]). A key observation is that some runs of x and y are merged in xy only when a string z 2 is both a suffix of x and a prefix of y (case a 2 on Figure 1). We first have this property : Proposition 5. Let Σ be an alphabet with |Σ| ≥ 2, and let x and y be strings on Σ such that |y| ≥ 1. Then there exists strings x ′ and y ′ on Σ such that

|x ′ | + |y ′ | = |x| + |y|, |y ′ | < |y| and r(x ′ ) + r(y ′ ) ≥ r(x) + r(y).
Proof. Let w be the longest string, eventually empty, such that w is a suffix of x and a prefix of y. Thus x = uw and y = wv for some strings u and v. Let x ′ = uwv and y ′ = w. Clearly |x ′ | + |y ′ | = |x| + |y| and |y ′ | ≤ |y|. Without loss of generality, we assume that y is not a suffix of x. (If it is not the case, we rewrite y into ȳ using an isomorphism of Σ onto itself.) Thus |y ′ | < |y|. Now we consider the runs of period p that were counted in r(x) + r(y). The runs with 2p characters ("a square") completely included in w were counted once in r(x) and once in r(y). Such runs are counted again once in r(x ′ ) and once in r(y ′ ). By definition of w, all the others runs counted in r(x) and r(y) are counted exactly once in r(x ′ ), without being merged.

To prove Proposition 1, we take two strings x 0 and y 0 of respective lengths m and n such that r(x 0 ) = ρ(m) and r(y 0 ) = ρ(n). Applying recursively Proposition 5 gives a finite sequence of pairs of strings (x 0 , y 0 ), (x 1 , y 1 ), . . . (x t , y t ) with r(x i ) + r(y i ) ≥ r(x i-1 ) + r(y i-1 ) and |y 0 | > |y 1 | > . . . > |y t | = 0 for some t.

a 1 ) a 2 ) b) q q b) q q 2p -1 2p -1 x y
Figure 1: a 1 ) Run with at least two periods included in x. a 2 ) Run with at least two periods included in x, and at least two periods included in y. b) "New runs" between x and y. To bound the new runs with period q ≤ p, it is sufficient to consider strings of length 4p -2.

Note that 4p characters would be required to exactly count the new runs.

Finally |x t | = |x 0 | + |y 0 | = m + n, and thus ρ(m + n) ≥ r(x t ) ≥ r(x 0 ) + r(y 0 ) = ρ(m) + ρ(n), proving Proposition 1.
Note that the proof of Franek et al. in [START_REF] Frantisek Franek | The maximum number of runs in a string[END_REF] was in a different context, and that no result leading to our Proposition 1 was stated as such in their paper.

On the number of microruns

In the following sections, p is fixed and we study the asymptotic behavior of the number of microruns ρ p (n)/n beyond the result of Proposition 3. In section 4.1, the idea is to bound the new runs created by the concatenation of two strings. In section 4.2, the idea is to count exactly the new runs created by the concatenation of a string and a character. Both sections provide new bounds or exact limits on the number of some microruns.

New runs obtained by string concatenation

Let x and y be two strings, and s be a run of xy with period q ≤ p. Then s is exactly in one of the following two cases (Figure 1):

• a) s has at least two periods included in x, or at least two periods included in y;

• b) s has strictly less than two periods included in x, and strictly less than two periods included in y.

We call the runs in the case b) the new runs between x and y, and we denote by NR p (x, y) the number of such runs. Then r p (xy) ≤ r p (x) + r p (y) + NR p (x, y), the inequality coming from the fact that a run from x can be merged with a run from y (case a 2 on Figure 1). We can bound the number of new runs, and thus have an upper bound on r p (xy): Proposition 6. Let z p = max{NR p (x, y) | |x| = |y| = 2p -1} the maximum number of new runs with period q ≤ p between strings of length 2p -1. Then, for every strings x and y of any length, we have NR p (x, y) ≤ z p . Proof. Any new run with period q ≤ p has at most 2q -1 ≤ 2p -1 characters in x, and in y (Figure 1).

p (m + n). Then ρ p (m + n) = r p (xy) ≤ r p (x) + r p (y) + NR p (x, y) ≤ ρ p (m) + ρ p (n) + z p .
Table 2 provides some values of z p for binary strings. An immediate bound on z p is z p ≤ z p-1 + 2. Knowing bounds on z p helps to characterize the asymptotic behavior of the number of microruns: Theorem 2. For any n, we have ℓ p ≤ ρ p (n)/n + z p /n, and thus

ℓ p - z p n ≤ ρ p (n) n ≤ ℓ p ≤ ℓ.
Proof. By Proposition 7, for any t ≥ 1, we have ρ p (tn) ≤ tρ p (n) + (t -1)z p . Thus ρ p (tn)/tn ≤ ρ p (n)/n + t-1 t z p /n. Taking this inequality to the limit, as t goes to infinity, gives the result.

Thus we know that the convergence of ρ p (n)/n to ℓ p is faster than z p /n. Note that we do not have a similar result for ρ(n), as we do not have a convenient way to bound ρ(m + n) like in Proposition 6.

As a side result of Theorem 2, we have new bounds of the number of some microruns. For example, brute-force computations give for binary strings z 9 = 7 and ρ 9 (34) = 26, thus ℓ 9 ≤ 33/34 = 0.970. For binary strings, this result is better than Lemma 2 of [START_REF] Crochemore | Maximal repetitions in strings[END_REF] which proved the n bound by the count of amortizing positions for centers of runs. The next section further improves this bound and finds exact values for some ℓ p 's.

The exact number of microruns

In this section, we propose to count exactly the number of microruns, by considering the concatenation of a string and a single character. Let x be a string, α ∈ Σ a character, and s be a run of xα with period q ≤ p. Then s is exactly in one of the following two cases (Figure 2) :

a 1 ) q q a 1 ) p p b) p p 2p x α
Figure 2: a 1 ) Runs of x, extended in α. b) "New runs" between x and α. Unlike in Fig. 1, there is no case a 2 ) where some runs are merged. To count the new runs with period q ≤ p, it is sufficient to consider the suffix of x of length 2p.

• a) s has at least two periods included in x;

• b) s has strictly less than two periods included in x.

As in the previous section, we call the runs in the case b) the new runs between x and α, and we denote by NR p (x, α) the number of such runs. As there is here no merging of runs, we have r p (xα) = r p (x) + NR p (x, α). In fact, the last 2p characters of x are sufficient to know NR p (x, α):

Proposition 8. Let x be a string with |x| = n ≥ 2p. Then we have NR p (x, α) = NR p (x[n -2p + 1 . . . n], α).
Proof. Any new run with period q ≤ p has at most 2q -1 ≤ 2p -1 characters in x. Any run of x with period q ≤ p extending in α has at a least 2 periods in the last 2p characters of x (Figure 2). Knowing the last 2p characters of x is thus sufficient to tell apart the two kinds of runs.

If v is a string of length 2p and n ≥ 2p, we define f n p (v) = max |x|+|v|=n r p (xv) as the maximum number of runs of all the strings ending with the suffix v. The function f n+1 p can be entirely determined from the functions f n p and NR p :

Proposition 9. If |w| = 2p -1, α ∈ Σ, and n ≥ 2p, then f n+1 p (wα) = max β∈Σ f n p (βw) + NR p (βw, α) .
Proof. To compute f n+1 p (wα) = max |x|+|wα|=n+1 r p (xwα), we suppose, without loss of generality, that the string x is of length at least one and we write x = yβ, where y is a string and β ∈ Σ. In particular, for every t ≥ 0, we have f n0+tk p = f n0 p + ts, that is ρ p (n 0 + tk) = ρ p (n 0 )+ts and finally lim t →∞ ρ p (n 0 +tk)/(n 0 +tk) = s/k. As ρ p (n)/n converges (Proposition 3), the Theorem is proved.

By computing f n p functions for successive n and by checking the additive periodicity condition of Theorem 3, one can have exact values of ℓ p for small p's. Table 3 lists results for p ≤ 14 on binary strings. Note that the periodicity on ρ p (n) can appear before the periodicity on f n p . For example, as soon as n ≥ 35, ρ 9 (n + 13) = ρ 9 (n) + 11, but the periodicity on f n 9 only starts at n 0 = 51.

Using the result of Crochemore and Ilie's Proposition 1 [START_REF] Crochemore | Maximal repetitions in strings[END_REF] for large runs, we get an upper bound on ρ(n)/n. For binary strings, the exact value ℓ 14 = 15/17 gives:

ℓ ≤ ℓ 14 + ℓ ≥15 ≤ 15 17 + 6 15 = 1.282...
Thus the number of runs in a binary string of length n is not more than 1.29n. This result was better than the 1.6n bound published in [START_REF] Crochemore | Maximal repetitions in strings[END_REF], but the better bound of 1.048n has now been published [START_REF] Crochemore | Towards a solution to the runs conjecture[END_REF]. Nevertheless, the values we give in the Table 3 are the first known exacts limits for such microruns.

Perspectives

The results on the asymptotic behavior of the functions ρ and ρ p of Theorems 1 and 2 simplify the research on lower and upper bounds. Moreover, the application of the Theorem 3 provides the first exact limits for the number of some microruns. We hope that these results will bring a better understanding of the number of runs and be a step towards proving the conjecture of [START_REF] Kolpakov | Maximal repetitions in words or how to find all squares in linear time[END_REF] (ℓ ≤ 1).

As Theorem 2 and 3 provide upper bounds or exact limits for some microruns, they can be used to bound the total number of runs. In both cases, this would require large evaluations of z p or f n p (w) values that could be improved by a more precise analysis, for example by taking inspiration from the methods of Crochemore and Ilie. Moreover, a better analysis could improve their 6/p bound on the number of macroruns ℓ ≥p . For the lower bound, it remains to be shown if one can find strings with more runs than those of [START_REF] Frantisek Franek | The maximum number of runs in a string[END_REF][START_REF] Franek | An asymptotic lower bound for the maximal-number-of-runs function[END_REF]. Although Theorem 1 also provides a way to have a lower bound on ρ(n)/n, all the computations we ran gave not better bounds than the 0.9445... bound of [START_REF] Matsubara | New lower bounds for the maximum number of runs in a string[END_REF].

Now an important question is if the actual value of ℓ can be found with such a separation between microruns and macroruns. The inequality ℓ ≤ ℓ p + ℓ ≥p+1 may be strict for some p. If this inequality is strict for several p's, the conjecture may be impossible to prove by this way if one choose a bad splitting period p.

Another open question is if one of the constants ℓ p = lim n →∞ ρ p (n)/n is equal to ℓ, or if, more probably, the limit ℓ is obtained by considering asymptotically runs with any period. Finally, it remains to be proven if strings on binary alphabets can always achieve the highest number of runs.

Proposition 7 .

 7 For any m and n, ρ p (m + n) ≤ ρ p (m) + ρ p (n) + z p . Proof. Let x and y be two strings such that |x| = m, |y| = n, and r p (xy) = ρ

Theorem 3 .

 3 = max |yβ|+|wα|=n+1 r p (yβwα) = max |yβ|+|wα|=n+1 (r p (yβw) + NR p (yβw, α)) = max |yβ|+|wα|=n+1 (r p (yβw) + NR p (βw, α)) (Proposition 8) = max β∈Σ max |y|+|βw|=n r p (yβw) + NR p (βw, α) = max β∈Σ f n p (βw) + NR p (βw, α) Once all the |Σ| 2p+1 values of NR p are computed, the above equation can be used to recursively determine any f n p function in O(n • |Σ| 2p+1 ) time. Then any number of microruns ρ p (n) = max |v|=2p f n p (v) follows. When n grows, an additive periodic behavior can emerge: If for some n 0 , k and s with 2p ≤ n 0 < n 0 + k, we have f n0+k p = f n0 p + s then, for any n ≥ n 0 , we have f n+k p = f n p + s, and ℓ p = s/k. Proof. If the property f n+k p = f n p + s is true for some n, then f n+1+k p (wα) = max β∈Σ f n+k p (βw) + NR p (βw, α) = max β∈Σ (f n p (βw) + s) + NR p (βw, α) = f n+1 p (wα) + s and the property is true for n + 1. By induction, it is true for every n ≥ n 0 .

Table

	n	5 6 7 8 9 10 11 12 13 14 15 16 17 18
	ρ(n) 2 3 4 5 5 6	7	8	8 10 10 11 12 13
	n	19 20 21 22 23 24 25 26 27 28 29 30 31
	ρ(n) 14 15 15 16 17 18 19 20 21 22 23 24 25
	Preprint submitted to Information and Computation	February 16, 2009

Table 2 :

 2 Values for zp for binary strings with worst-case examples of length ≤ 4p -2.

Table 4 :

 4 Values of ρp(n) for binary strings. For each n, the value in bold shows the smallest period p such that ρp(n) = ρ(n).

	n	ρ(n)	1, 2	3	4	5	6	7	8	9	10	11	12	13	14
	5	2	2
	6	3	3	3
	7	4	3	4
	8	5	4	4	5
	9	5	4	5	5
	10	6	5	6	6	6
	11	7	5	6	7	7
	12	8	6	7	8	8	8
	13	8	6	8	8	8	8
	14	10	7	9	9	9	9	10
	15	10	7	9	9	10	10	10
	16	11	8	10	10	11	11	11	11
	17	12	8	11	11	11	11	12	12
	18	13	9	12	12	12	12	13	13	13
	19	14	9	12	13	13	13	14	14	14
	20	15	10	13	13	14	14	14	15	15	15
	21	15	10	14	14	14	14	15	15	15	15
	22	16	11	15	15	15	15	15	16	16	16	16
	23	17	11	15	16	16	16	16	17	17	17	17
	24	18	12	16	16	17	17	18	18	18	18	18	18
	25	19	12	17	17	18	18	18	18	18	18	19	19
	26	20	13	18	18	18	18	19	19	19	19	19	19	20
	27	21	13	18	18	19	19	20	20	20	20	21	21	21
	28	22	14	19	19	20	20	21	21	21	21	21	21	22	22
	29	23	14	20	20	21	21	21	22	22	22	22	22	23	23
	30	24	15	21	21	21	21	22	23	23	23	23	23	24	24
	31	25	15	21	21	22	22	23	24	24	24	24	24	25	25
	32	26	16	22	22	23	23	24	25	25	25	25	25	26	26
	33	27	16	23	23	24	24	25	26	26	26	26	26	27	27
	34	27	17	24	24	25	25	26	26	26	26	26	26	27	27
	35	28	17	24	24	25	25	26	27	27	27	27	27	28	28

= f n p + s for some s and k > 0. The value ℓ 13 required three hours of computation on a standard 2 GHz workstation. This time is almost entirely spent in the initial computation of the 2 2p+1 values of the function NRp, obtained by aggregate calls to mreps [START_REF] Kolpakov | mreps: Efficient and flexible detection of tandem repeats in dna[END_REF]. The successive computations of f n p are done in a few seconds.


Mathieu Giraud
email: mathieu.giraud@lifl.fr

Jean-Stéphane Varré
email: jean-stephane.varre@lifl.fr

Parallel Position Weight Matrices Algorithms

Position Weight Matrices (PWMs) are broadly used in computational biology. The basic problem, SCAN, aims to find the occurrences of a given PWM in large sequences. Some other PWM tasks share a common NP-hard subproblem, SCOREDISTRIBUTION. The existing algorithms rely on the enumeration on a large set of scores or words, and they are mostly not suitable for parallelization. We propose a new algorithm, BUCKETSCOREDISTRIBUTION, that is both very efficient and suitable for parallelization. We bound the error induced by this algorithm. We realized a GPU prototype for SCAN and BUCKETSCOREDISTRIBUTION with the CUDA libraries, and report for the different problems speedups of 21× and 77× on a Nvidia GTX 280.

Introduction

Position Weight Matrices (PWMs). Position Weight Matrices (PWMs) are broadly used in computation biology to model conserved sequence patterns. The most common application of PWMs is about gene regulation: the transcription of a gene is controlled by regulatory proteins that bind to transcription factor binding sites (TFBSs) on DNA. The TF-BSs are located mostly in non-coding regions preceding the genes. Discovering these motifs is rather difficult because of their very low information content. The reference databases JASPAR [START_REF] Sandelin | JASPAR: an open-access database for eukaryotic transcription factor binding profiles[END_REF] and TRANSFAC [START_REF] Wingender | TRANSFAC: an integrated system for gene expression regulation[END_REF] respectively contain 123 and 856 matrices of TFBSs. New sequencing technologies enable large-scale mapping of DNA-protein interactions : those "ChIP-Seq" methodologies produce new collections of sequences and matrices [START_REF] Shendure | Next-generation DNA sequencing[END_REF], [START_REF] Robertson | Genome-wide profiles of STAT1 DNA association using chromatin immunoprecipitation and massively parallel sequencing[END_REF].

Given a finite alphabet Σ and a positive integer m, a PWM M is a matrix with |Σ| rows and m columns (Figure 1). The coefficient M (p, x) gives the score at position p for the letter x in Σ. The PWM defines a function from Σ m to R, that associates a score to each word u = u 1 u 2 . . . u p of Σ m :

Score M (u) = m p=1 M (p, u p ),
Let α be a score threshold. We say that M has an occurrence in a text T at position k if Score M (T k . . . T k+m-1 ) ≥ α.

The most recurrent task is to predict binding sites in a large DNA sequence, that is to look for occurrences of a PWM given a text and a score threshold. This basic task is often involved into a more general analysis pipeline. After the occurrences have been found, a filter is applied in order to refine the results. Tools which compute occurrences often output the statistical significance, called the P-value, of each occurrence. The P-value Pv M (s) is the probability that the background model achieves a score at least equal to s: it is the proportion of words u (randomly chosen according to the observed letter frequencies) whose score is greater than s. If the background model is chosen with identically and independently distributed character symbols, the P-value is simply

Pv M (s) = | {u ∈ Σ m | Score M (u) ≥ s} | |Σ| m
but other background models can be used by weighting the words u with their relative probability in the background model.

To decide if a PWM occurs at a position in a text, the score threshold has to be chosen. Generally, a P-value p is chosen independently from the matrix and the background model: this P-value reflects the expected number of occurrences that will be found in the text [START_REF] Claverie | The statistical significance of nucleotide position-weight matrix matches[END_REF]. Then the score threshold is computed for a given matrix M : the goal is to find the score α such that Pv M (α) = p. The two problems of computing the P-value for a score and of computing the score from a given P-value are NP-hard [START_REF] Zhang | Computing exact p-values for DNA motifs[END_REF], [START_REF] Touzet | Efficient and accurate p-value computation for position weight matrices[END_REF].

Another task that attracted interest in the past few years is the design of a method to compare matrices. Column-tocolumn comparisons with correlation coefficients have been proposed by several authors [START_REF] Schones | Similarity of position frequency matrices for transcription factor binding sites[END_REF]- [START_REF] Gupta | Quantifying similarity between motifs[END_REF]. A better method is to consider that two matrices are similar if their occurrences are almost the same [START_REF] Pape | Natural similarity measures between position frequency matrices with an application to clustering[END_REF]. Being able to compare matrices serve several goals. Firstly, it can be used to detect if two matrices are similar. This is useful to remove redundancy from databases or to test if a new matrix has a similar one in a database [START_REF] Kielbasa | Measuring similarities between transcription factor binding sites[END_REF]- [START_REF] Pape | Natural similarity measures between position frequency matrices with an application to clustering[END_REF]. Secondly, it can be used to improve the SCAN problem (see page 3). Computation bottleneck in PWM algorithms. The exponential nature of some PWM problems is a limiting factor for using matrices of medium or large length. Score threshold computations for matrices whose length is greater than 15 usually require several seconds, and hours or days for matrices of length greater than 20. The JASPAR and TRANSFAC databases already contain almost 200 matrices of length 15 to 30, that is 20% of their total number of matrices. Moreover, the new techniques using data from next-generation sequencers should produce longer matrices.

In this paper, we give some parallel solutions to the common PWM problems. We propose a prototype implementation on graphic processing units (GPUs) to test the ability to compute with longer matrices. We think that GPUs are a first step toward new massively many-cores architectures: we will discuss how our algorithms can be applied to other parallelization techniques. General-purpose computation on GPU. Everyone can have some teraflops of cheap computing power with the recent Graphics Processing Units (GPUs). GPUs were used in bioinformatics since 2005 for phylogenetic studies [START_REF] Charalambous | Initial experiences porting a bioinformatics application to a graphics processor[END_REF], then for multiple sequence alignment based on an optimized Smith-Waterman implementation [START_REF] Liu | GPU-ClustalW: using graphics hardware to accelerate multiple sequence alignment[END_REF]. The CUDA libraries, first released in 2007 [START_REF]Nvidia CUDA programming guide 2[END_REF], have deeply simplified the development on GPUs. Recent papers provide speedups on bioinformatics applications involving suffix trees [START_REF] Schatz | High-throughput sequence alignment using graphics processing units[END_REF] or again Smith-Waterman comparisons [START_REF] Manavski | CUDA compatible GPU cards as efficient hardware accelerators for Smith-Waterman sequence alignment[END_REF].

The current Nvidia architectures [START_REF]Nvidia CUDA programming guide 2[END_REF] offer two level of parallelism. For the coarse-grained level, several multiprocessors execute blocks of independent computations. Each multiprocessor is then a kind of large SIMD device, able to process several different fine-grained threads at a given time. All those threads are executing exactly the same instructions: if a divergence inside a conditional expression occurs, the two branches are serialized. A 16 KB shared memory is available for the threads in a same block. This local memory is very fast and should be used to maximize the efficiency. Contents. We investigate the four following problems in terms of parallel programming, especially on GPU architectures.

• SCAN. Given a matrix M with a score threshold α and a text T , find all the occurrences of M in T . • SCORETOPVALUE. Given a matrix M and a score s, compute the P-value Pv M (s). • PVALUETOSCORE. Given a matrix M and a P-value p, compute the score s such that Pv M (s) = p.

• COMPARE. Given two matrices M and M of same length with respective score thresholds α and α , compute the number of words with a score greater than α for M and than α for M . Section 2 addresses the SCAN problem, for which a simple parallelization is very efficient. Section 3 is related to the three problems SCORETOPVALUE, PVALUETOSCORE and COMPARE. Those three problems share a common NPhard sub-problem, SCOREDISTRIBUTION, whose current solutions are not suitable for parallelization. We propose a new algorithm, BUCKETSCOREDISTRIBUTION, that is more efficient than the existing algorithms (both in speed and in precision) and suitable for parallelization. Section 4 reports our prototype implementation of the three first algorithms with the CUDA libraries [START_REF]Nvidia CUDA programming guide 2[END_REF], and discuss our solution compared to other parallelizing techniques.

Looking for occurrences of a matrix

Finding the occurrences of a matrix M of length m in a text T of length n may be done by a naive algorithm in O(mn) time: for each position k of the text, Score M (T k . . . T k+m-1 ) is computed. Optimizations. In 2000, [START_REF] Wu | Fast probabilistic analysis of sequence function using scoring matrices[END_REF] proposed an improvement: for a given word u ∈ Σ m , the computation of Score M (u) can be stopped as soon as

Score M [1..j] (u 1 . . . u j ) < α -max M [j + 1..m] with max M [j + 1..m] = m k=j+1 max M [k], where max M [k]
is the maximal score of the k-th column of the matrix and α the score threshold. Indeed, when the above condition is met, the score of u cannot be greater than α. This property does not change the O(mn) worst-time complexity, but gives an average O(m n) time complexity, where m is the average stop position. We reference this algorithm as the Lookahead Strategy Algorithm (LSA for short).

In 2006, [START_REF] Liefooghe | Large scale matching for position weight matrices[END_REF] proposed to precompute an index for the SCAN of one or several matrices. The main idea was to split the matrix into sub-matrices called slices of length (typically 7 or 8, depending on the available memory). For each slice, the |Σ| scores for each word are computed and stored into a table. The time complexity remains O(m n), but with only O(m n/ ) memory accesses. When several matrices are scanned, the table is organized such that scores for the set of matrices are in the same memory location, thus avoiding memory latencies. On a scan involving a large set of matrices, a practical 8× speedup is obtained compared to the LSA. For a single matrix, no significant speedup is obtained. We call this idea the Slices Strategy.

Lastly, similarities between matrices can lead to another improvement when searching for occurrences: one can avoid to look for occurrences of each matrix but only for occurences of one representative matrix [START_REF] Kielbasa | Measuring similarities between transcription factor binding sites[END_REF], [START_REF] Liefooghe | Large scale matching for position weight matrices[END_REF]. Occurrences of the other matrices are then computed only for positions where the representative matrix occurred.

Preprocessing and indexation. Ideas can also be borrowed from the algorithms searching a pattern in a text, as with the classical Aho-Corasick [START_REF] Aho | Efficient string matching: An aid to bibliographic search[END_REF], Knuth-Morris-Pratt (KMP) [START_REF] Knuth | Fast pattern matching in strings[END_REF], or Boyer-Moore [START_REF] Boyer | A fast string searching algorithm[END_REF] algorithms and their variants. This is more difficult than in the usual pattern matching case, as the combination of the scores does not always allow large shifts in the matrices. With KMP, [START_REF] Liefooghe | Self-overlapping occurrences and Knuth-Morris-Pratt algorithm for weighted matching[END_REF] obtains a 2× -3× improvement on the LSA method. With Aho-Corasick, [START_REF] Pizzi | Fast search algorithms for position specific scoring matrices[END_REF] obtains a 2× -5× speedup. Nevertheless, the size of the automaton is such that it does not fit in memory for matrices longer than 15.

Instead of indexing the matrices, another way to speedup the SCAN problem is to preprocess the text: [START_REF] Dorohonceanu | Accelerating Protein Classification Using Suffix Trees[END_REF] uses suffix trees (speedup 2× -5×), [START_REF] Beckstette | Fast index based algorithms and software for matching position specific scoring matrices[END_REF] uses suffix arrays and [START_REF] Freschi | Using sequence compression to speedup probabilistic profile matching[END_REF] compress the text (speedup 2× -5×). Some parallel implementations of suffix trees have been reported [START_REF] Schatz | High-throughput sequence alignment using graphics processing units[END_REF], but those parallelizations are especially difficult due to the nonlocality of memory accesses and the structure of the tree.

In fact, those evolved data structures do not fit well with pattern matching with errors. The PWM SCAN problem is far more difficult, as it can be seen as a generalized pattern matching with a complex error function. Other solutions could use seed-based indexing, in particular spaced seeds that handle better errors [START_REF] Brown | Bioinformatics Algorithms: Techniques and Applications[END_REF].

A simple parallel implementation. The parallelization of the SCAN problem can be done easily by splitting the different positions of the text across several threads (Figure 2). At a given time, each thread computes the score of one word.

The Slices Strategy could bring a very small improvement. On the contrary, the LSA as well as KMP-like improvements do not apply on SIMD-like architectures such as GPU. Their efficiency is indeed based on a variable number of iterations in the most inner loops. As even close words do not always lead to the same number of iterations, different threads with different words will diverge most of the time, thus providing a very bad parallel performance.

Computing score threshold, P-value and comparing matrices

Methods for score threshold and P-value computation. The computation of the P-value (SCORETOPVALUE) can be done using probability generating functions or dynamic programming [START_REF] Claverie | The statistical significance of nucleotide position-weight matrix matches[END_REF]- [START_REF] Touzet | Efficient and accurate p-value computation for position weight matrices[END_REF], [START_REF] Staden | Methods for calculating the probabilities of finding patterns in sequences[END_REF], [START_REF] Rahmann | Dynamic programming algorithms for two statistical problems in computational biology[END_REF]. In both cases, the time complexity is O(S), where S is the number of possible different scores. Let Q M (s) be the probability to achieve exactly the score s. For a given score s , the P-value is obtained with the equation:

Pv M (s ) = s≥s Q M (s)
Computing SCORETOPVALUE (that is the score associated to a given P-value) can be done by adding the values of Q M (s ) from the maximal score until the desired P-value is obtained.

Let M [1.

.i], 0 ≤ i ≤ n, denote the matrix consisting of the columns 1 to i of M . The matrix M [1..0] is the empty matrix. The Q M [1..i] score distribution can be expressed from the Q M [1..i-1] score distribution by the following dynamic programming algorithm, where p(x) is the probability of the letter x in the background model [START_REF] Claverie | The statistical significance of nucleotide position-weight matrix matches[END_REF]:

Q M [1..0] (s) = 1 if s = 0 0 otherwise Q M [1..i] (s) = x∈Σ Q M [1..i-1] (s -M (i, x)) × p(x)
If the matrix has non-negative integer coefficient values, then S, the number of possible different scores, is bounded by

m i=1 max M [i].
It follows that known algorithms are pseudo-polynomial. However, PWMs are built from logratios and do not fulfill this constraint: S can be as large as |Σ| m , and thus the worst-case time complexity is O(|Σ| m ). Some methods [START_REF] Touzet | Efficient and accurate p-value computation for position weight matrices[END_REF], [START_REF] Pape | Natural similarity measures between position frequency matrices with an application to clustering[END_REF], [START_REF] Beckstette | Fast index based algorithms and software for matching position specific scoring matrices[END_REF] round the coefficients of the matrix to maintain S low. Even if those methods claim to compute exact P-values, this rounding induces an error on the P-value and on the score threshold [START_REF] Touzet | Efficient and accurate p-value computation for position weight matrices[END_REF]. As an example, the implementation of MOSTA [START_REF] Pape | Natural similarity measures between position frequency matrices with an application to clustering[END_REF] rounds each column of the matrix to the nearest multiple of ε = 0.05, and thus the total score has an error of at most mε. A new algorithm to compute the score distribution. We propose to compute the score distribution by splitting the matrix M in N slices M 1 , M 2 , . . . M N , and by combining the score distributions of the slices (Figure 3):

Q M (s) = s1+...+s N =s (Q M1 (s 1 ) × . . . × Q M N (s N ))
In the following algorithm, we use tables with B elements called buckets. For any slice M i , the scores are in the range [min M i , max M i ], and we store this distribution in a table with B buckets, thus rounding down the scores to the nearest multiple of ∆ Mi , where ∆ Mi = (max M i -min M i )/B. We denote by s the discretized score of s, S M the set of all possible scores the matrix M can achieve and S M the set of all possible discretized scores the matrix M can achieve.

step k + 1 step k N threads N 1 2 T Score(M, u 2 ..u m+1 ) . . . . . . Score(M, u 1 ..u m ) 1 N Score(M, u N ..u N +m-1 ) . . . . . . Score(M, u N +2 ..u N +m+1 ) Score(M, u N +1 ..u N +m ) Score(M, u 2N ..u 2N +m-1 ) 1 N + 1 2 N N 2

Algorithm BUCKETSCOREDISTRIBUTION

• For each slice M i , compute a score distribution Q Mi by enumeration of 4 m/N words, rounding the scores of each word. The result is stored in a table Q Mi with B entries. For a score s in S Mi ,

Q Mi [s] = s ∈ S M i s = s Q Mi (s ) (1)
Note that we round the score of each word, and not of each column, thus keeping the error low. • A score distribution for M is then computed by recursively merging the score distributions of the N slices (Figure 3). For a score s in S M1⊕2 ,

Q M1⊕2 [s] = s1 ∈ S M 1 s2 ∈ S M 2 s1 + s2 = s Q M1 [s 1 ] × Q M2 [s 2 ]
All the tables have B buckets. One merge can be done in O(B 2 ) time, for a total time of at most O(N B 2 ).

The merge can also be computed in O(N B log B) time using rapid convolution algorithms [START_REF] Knuth | The Art of Computer Programming[END_REF], but here this step is not a limiting factor.

Assuming that m is large enough, the time complexity of BUCKETSCOREDISTRIBUTION is O(N 4 m/N ), enabling the study of matrices N times larger than the previous methods.

Once Q M 1⊕2⊕•••⊕N is known, the final step of computation, in O(B) time, depends on the problem. For SCORE-TOPVALUE, we obtain an approximate P-value,

Pv M [s], maximal 3∆ M 2∆ M ∆ M M 2 M = M 1⊕2⊕3⊕4 M 1⊕2 M 3⊕4 M 3 M 1 M 4 error
Pv M [s] = s ∈ S M s ≥ s Q M [s ]
For PVALUETOSCORE, the score threshold s in S M is the biggest one such that

p ≤ s ∈ S M s ≥ s Q M [s ]
Precision evaluation. We now bound the error induced by the score discretization and the error induced by the convolution to show that they are similar to the errors of other algorithms. At the first step, the maximum error when discretizing the scores for a slice M i is ∆ Mi . When combining two slices M i and M j into M i⊕j , the equation ( 1) is not more valid: now the maximum total error is

∆ Mi + ∆ Mj + ∆ Mi⊕j
where ∆ Mi⊕j is the maximum error when discretizing the result. As the combined scores are in the range

[min M i + min M j , max M i + max M j ], we have ∆ Mi⊕j = ∆ Mi + ∆ Mj .
The maximum total error is thus 2∆ Mi⊕j . With N slices, the maximal total error on M is

1 + log N ∆ Mi = 1 + log N ∆ M with ∆ M = (max M -min M )/B. This O((log N )/B)
maximal error is on the scores. The actual error on the number of words depends on s and on the score distribution of M : it is at most the number of words in the 1 + log N buckets (gray area on the Figure 4). For the PVALUETOSCORE and the SCORETOPVALUE problems, the error on the P-value is thus

Pv M [ s ] ≤ Pv M (s) ≤ Pv M [ s + 1 + log N ∆ M ] words 3∆ M
error on score error on This new algorithm is perfectly suited for parallelization, as the word enumerations can be split across different independent computations. In our GPU prototype implementation, the enumeration of 4 m/N words is split on 4 β blocks with 4 τ threads by block, leaving 4 µ words to enumerate within each thread, with µ = m/N -(β + τ ) (Figure 5, on the left). At a given time, all the threads of a same block are enumerating the same µ leftmost characters.

From one word to another, the score is evaluated only on the modified positions, thus bringing no divergence between the threads of a same block. Then each thread increments one of the B buckets (Figure 5, on the right). As several threads can increment the same bucket at a given time, atomic instructions or similar mechanisms must be used there. The final merging operations, in O(N B 2 ) time, can be performed on the host. Similarity between two matrices. The computation of similarity between two matrices is very similar. Given two matrices 1 M and M of length m with their respective score thresholds α and β, the goal is to measure TP M M , the number of true positive words u ∈ Σ m such that Score M (u) ≥ α 1. If the two matrices have different lengths m and m , one has to choose an alignment, i.e. a shift between the two matrices. Then it is sufficient to pad the two matrices with zeros to obtain two matrices of same length. When one wants to compute such a similarity score, one has to compute TP for all possible m + m shifts. and Score M (u) ≥ β. This number can be computed with the following equation:

TP M M (α, β) = s≥α,s ≥β Q M M (s, s )
where Q M M (s, s ) is now the probability to achieve exactly the score s with M and the score s with M . False positives, true negatives and false negatives are computed in a similar way. The score distribution is now expressed by the following recurrence [START_REF] Liefooghe | Large scale matching for position weight matrices[END_REF]:

Q M [1..0] M [1..0] (s, s ) = 1 if s = 0 or s = 0 0 otherwise Q M [1..i] M [1..i] (s, s ) = x∈Σ Q M [1..i-1] M [1..i-1] (s -M (i, x), s -M (i, x)) × p(x)
The same BUCKETSCOREDISTRIBUTION algorithm applies, but now the B buckets induce a maximal O((log N )/ √ B) error on words.

Results and discussion

Testing environment. We benchmarked the three SCAN, SCORETOPVALUE and PVALUETOSCORE algorithms with the CUDA 2.0 libraries from Nvidia [START_REF]Nvidia CUDA programming guide 2[END_REF]. The sources or our implementation are available at http://bioinfo.lifl.fr/cudapwm. Two GPU were tested: the Nvidia GeForce 8800 (16 × 8 cores, 1.3 GHz, 768 MB RAM), and the Nvidia GTX 280 (30 × 8 cores, 1.3 GHz, 1 GB RAM).

The host system and the CPU benchmarks were produced on an Intel Core 2 Duo 6600 (2.40 GHz) with 3 GB RAM and 4 MB cache (the CPU has two cores, but only one core was used in the benchmarks). The compiler was Nvidia nvcc used with the -O3 option.

Benchmarks were done on real data (Table 1 and Figure 8) and on random data (Figures 6 and7). We found no significant difference in terms of speedups between those datasets.

Algorithm SCAN

For the parallel SCAN, the target sequence was always in the shared memory, enabling different threads in the same block to work on the same data (see Figure 2). Figure 6 and Table 1 detail the results. Using the GPU implies an overhead of 0.27s on a 225 Mbp chromosome (Table 1). This overhead is mostly due to the transfer of the sequence data to the GPU (even large collections of matrices are small compared to megabytes of sequence data). The overhead is negligible starting with matrices of length 30, or as soon as a collection of matrices is scanned, as there is in this case only one transfer of the sequence.

The Slices Strategy does not bring here a significant improvement. As expected, the strategies with some divergence in the threads (LSA and KMP) do not provide any speedup (results not shown): the naive algorithm is here the more efficient algorithm to parallelize. bases, release hg18). Times are in seconds.

Figures 6 details the speedup of the GPU implementations compared to a 1-thread CPU implementation. The best speedups, 16× on GeForce 8800 and 21× on GTX 280, are obtained starting from 30 • 10 9 positions computed. Those speedups should be compared to the 2× to 8× speedups of the methods cited on page 3, and to a maximum 8× speedup using 128-bit SIMD instructions using a 16-bit precision.

Algorithms based on BUCKETSCOREDISTRI-BUTION

Implementation and speedups. We required that the bucket table fit in the shared memory, leading to B = 3600 32bit buckets for one matrix. Figure 7 details the results for PVALUETOSCORE with N = 4. The computation times are the same for the SCORETOPVALUE problem. Due to some overheads in initializations and in memory transfers between the host and the GPU (results not shown, similar to Table 1), a significant speedup starts only from a slice size of = 48, that is from the enumeration of 4 12 words. Because of its increased number of cores, the GTX 280 gives a 80% increased speedup on GeForce 8800. For = 65, the CPU takes 1426 seconds: the speedup is 3.6× for the GeForce 8800 and 9.8× for the GTX 280. The main limitation in those speedups is the bucket incrementation: most of the time, several threads are incrementing the same bucket. On the GeForce 8800, this incrementation is serialized between the threads of a same block at the end of each score computation. An improved algorithm could have here better results. On the GTX 280, the atomic instructions in shared memory give an additional 8× speedup compared to the serialized incrementation. For = 65, the total speedup is thus 77×. Precision. Figure 8 details the errors on P-value when computing PVALUETOSCORE for all matrices of JASPAR with length ≥ 10. We accept a relative error of 0.1, suitable for the usage of PWMs in computational biology (for a Pvalue of 10 -5 , the error will be at most 10 -6 ). Here 83 out of 101 matrices (of which all matrices of length ≥ 13) fulfill this precision for a P-value of 10 -5 . The relative error decreases when the P-value is higher and when the matrix is longer. In both cases, this is because the position of the score s is more "on the left" relative to the shape of the score distribution. Other parallelization techniques. On any multi-cores architecture, the BUCKETSCOREDISTRIBUTION 4 m/N enumerations can be split among several cores. Even on general purpose CPUs, SIMD techniques can benefit from BUCK-ETSCOREDISTRIBUTION, as only one common memory access is needed between all the 4 τ "threads" that process the same leftmost characters at a given iteration (see page 5). Nevertheless, at 16-bit precision, the maximum theoretical speedup using SSE 128-bit SIMD extensions is 8×, far beyond the 77× speedup achieved on GTX 280.

Conclusion and perspectives

We proposed a parallel SCAN of a Position Weight Matrix (PWM) in a text, and a new algorithm, BUCKETSCOREDIS-TRIBUTION, that computes the score distribution of a PWM and that resolves the SCORETOPVALUE, PVALUETOSCORE and COMPARE problems. The error on P-value induced by the BUCKETSCOREDISTRIBUTION algorithm for the SCORETOPVALUE and PVALUETOSCORE problems is not more than O((log N )/B), where B is the number of buckets used to store the distributions. The BUCKETSCOREDIS-TRIBUTION can be adapted to any parallel architecture, as it involves large blocks of independent computations. On a GPU, threads simultaneously enumerate neighbor words without divergence. The only bottleneck was in the buckets incrementation. The best speedup is here obtained through the use of atomic instructions of the GTX 280. Further work could be done on benchmarking the COM-PARE problem on real data to see if the O((log N )/ √ B) error is suitable for clustering applications. Other perspectives include testing and improving the BUCKETSCOREDIS-TRIBUTION algorithm on other many-cores architectures, in particular through the new OpenCL standard [START_REF]The Khronos Group, OpenCL 1.0 specification[END_REF].

Figure 1 .

 1 Figure 1. A Position Weight Matrix (PWM) modeling the CREB1 transcription factor binding site (from the JASPAR database). The coefficients are log-odds ratios of letter frequencies: they indicate affinities between letters and positions at the binding site.

Figure 2 .

 2 Figure 2. Parallel GPU SCAN. N threads compute the scores of words from position 1 to N , then all threads are shifted by N characters onto the sequence.

Figure 3 .

 3 Figure 3. Algorithm BUCKETSCOREDISTRIBUTION

Figure 4 .

 4 Figure 4. Error on P-value when computing BUCK-ETSCOREDISTRIBUTION Parallelization.This new algorithm is perfectly suited for parallelization, as the word enumerations can be split across different independent computations. In our GPU prototype implementation, the enumeration of 4 m/N words is split on 4 β blocks with 4 τ threads by block, leaving 4 µ words to enumerate within each thread, with µ = m/N -(β + τ ) (Figure5, on the left). At a given time, all the threads of a same block are enumerating the same µ leftmost characters. From one word to another, the score is evaluated only on the modified positions, thus bringing no divergence between the threads of a same block. Then each thread increments one of the B buckets (Figure5, on the right). As several threads can increment the same bucket at a given time, atomic instructions or similar mechanisms must be used there. The final merging operations, in O(N B 2 ) time, can be performed on the host. Similarity between two matrices. The computation of similarity between two matrices is very similar. Given two matrices 1 M and M of length m with their respective score thresholds α and β, the goal is to measure TP M M , the number of true positive words u ∈ Σ m such that Score M (u) ≥ α

Figure 5 .

 5 Figure 5. Parallel GPU BUCKETSCOREDISTRIBUTION. Each thread enumerates 4 µ words.

Figure 6 .

 6 Figure 6. Results for SCAN (random matrices, random sequence of length 500 • 10 6 ).

Figure 7 .

 7 Figure 7. Results for PVALUETOSCORE, with N = 4 slices on random matrices (other N give similar speedups).

Figure 8 .

 8 Figure 8. Relative error on P-value when computing PVALUETOSCORE with N ∈ [2, 4] for all 101 matrices of JASPAR of length ≥ 10.

Table 1 .

 1 Results for parallel SCAN on GeForce 8800, on the chromosome 1 of the human genome (225 • 10 6

	matrix	length init + I/O kernel	total
	JASPAR 0075	5	0.27	0.27	0.54
	JASPAR 0023	10	0.27	0.34	0.61
	JASPAR 0106	20	0.27	0.48	0.75
	random	40	0.27	2.04	2.31
	random	80	0.27	2.88	3.15
	JASPAR collection	1304	0.27	36.60 36.87

Acknowledgements

This research was carried out with a grant "Action Incitative LIFL" and through the "NVIDIA Professor Parternship" program. We thank Jérémie Allard (INRIA Lille) for giving us access to his GTX 280 card, and anonymous reviewers for comments on a previous version on the manuscript.


Peter Steffen
email: psteffen@techfak.uni-bielefeld.de

Robert Giegerich
email: robert@techfak.uni-bielefeld.de

Mathieu Giraud
email: mathieu.giraud@lifl.fr

GPU Parallelization of Algebraic Dynamic Programming

Algebraic Dynamic Programming (ADP) is a framework to encode a broad range of optimization problems, including common bioinformatics problems like RNA folding or pairwise sequence alignment. The ADP compiler translates such ADP programs into C. As all the ADP problems have similar data dependencies in the dynamic programming tables, a generic parallelization is possible. We updated the compiler to include a parallel backend, launching a large number of independent threads. Depending on the application, we report speedups ranging from 6.1× to 25.8× on a Nvidia GTX 280 through the CUDA libraries.

Introduction

Dynamic programming in bioinformatics. In biological sequence analysis, there arise numerous combinatorial optimization problems that are solved by dynamic programming. Pattern matching in DNA or protein sequences, comparison for local or global similarity, and structure prediction from RNA sequences are frequent tasks, as well as the modeling of families of proteins and RNA structures with the widely used Hidden Markov Models (HMMs) and stochastic context free grammars (SCFG), respectively [START_REF] Durbin | Biological Sequence Analysis[END_REF]. The scoring schemes associated with these optimization problems can be quite sophisticated. The thermodynamic model for RNA structure prediction, for example, has more than thousand parameters. This requires elaborate case analysis. Objective functions often ask for more than a single answer, such as the best non-overlapping pattern hits to a genome above a certain score threshold. Finally, biological sequences tend to be long (from 77 characters for a tRNA, 10000 for a gene, 3 * 10 6 for a bacterial genome, to the 3 * 10 9 nucleotides of a mammalian genome such as human or mouse). The time and space requirements for a dynamic programming algorithm are often limiting factors for the problems the biologists need to solve. The development of reliable and efficient dynamic programming algorithms in bioinformatics is a recurring challenge, in sharp contrast to the simplicity suggested by the textbook examples of dynamic programming which we use to teach computer science students.

Algebraic dynamic programming. In all these optimization problems, the logical problem decomposition follows the decomposition of the input sequence into subwords. It has been observed early that the resulting dynamic programming recurrences strongly resemble those of a Cocke-Younger-Kasami [START_REF] Aho | The Theory of Parsing, Translation and Compiling[END_REF] type parsing algorithm [START_REF] Searls | Linguistic approaches to biological sequences[END_REF]. Pursuing this analogy, we have developed an algebraic style of dynamic programming (ADP) over sequential data. The search space of the optimization problem at hand is described by a yield grammar, which is a regular tree grammar generating a tree language, and implicitly a context-free language as the set of leaf sequences of these trees. Scoring and optimization are described by an evaluation algebra, which interprets the tree operators as functions that compute local score contributions, and hence solve larger problems when given optimal solutions of smaller ones, consistent with the general paradigm of dynamic programming. This leads to a complete specification of dynamic programming algorithms on a rather high level of abstraction. General-purpose computation on GPU. For a few years, issues with heat dissipation have prevented the processors from having higher frequencies. One of the answers to maintain the Moore's Law is the use of parallel processing with massively manycore architectures. Graphic processing units (GPUs) are a first step towards those architectures, and recent trends blur the line between such GPUs and multi-core processors.

GPUs were used in bioinformatics since 2005 for phylogenetic studies [START_REF] Charalambous | Initial experiences porting a bioinformatics application to a graphics processor[END_REF], then for multiple sequence alignment based on an optimized Smith-Waterman implementation [START_REF] Liu | GPU-ClustalW: using graphics hardware to accelerate multiple sequence alignment[END_REF]. The CUDA libraries, first released in 2007 [START_REF][END_REF], have deeply simplified the development on GPUs. Recent papers provide speedups on applications involving suffix trees [START_REF] Michael C Schatz | Highthroughput sequence alignment using graphics processing units[END_REF] or again Smith-Waterman comparisons [START_REF] Ligowski | An efficient implementation Smith-Waterman algorithm on GPU using CUDA, for massively parallel scanning of sequence databases[END_REF][START_REF] Liu | CUDASW++: optimizing Smith-Waterman sequence database searches for CUDA-enabled graphics processing units[END_REF][START_REF] Svetlin | CUDA compatible GPU cards as efficient hardware accelerators for Smith-Waterman sequence alignment[END_REF], error correction in DNA short-reads sequencing [START_REF] Shi | Accelerating error correction in high-throughput short-read DNA sequencing data with CUDA[END_REF], computation with position weight matrices [START_REF] Mathieu | Parallel position weight matrices algorithms[END_REF], RNA folding [START_REF] Rizk | GPU accelerated RNA folding algorithm[END_REF], and neighbor-joining trees for multiple sequence alignments [START_REF] Liu | Parallel reconstruction of neighbor-joining trees for large multiple sequence alignments using CUDA[END_REF].

The current Nvidia architectures [START_REF][END_REF] offer two levels of parallelism. For the coarse-grained level, several multiprocessors execute blocks of independent computations. Each multiprocessor is then a kind of large SIMD device, able to process several different fine-grained threads at a given time. All those threads are executing exactly the same instructions: if a divergence in a condition occurs, the branches of the condition are serialized.

Contents. In this paper we describe an approach to extend the ADP compiler such that it generates parallel programs (Figure 1). This approach has been implemented in the ADP compiler with the CUDA libraries [START_REF][END_REF]. Our new contribution is thus a generic method to create parallel CUDA programs for bioinformatics applications -classical and yet-to-be written ones. The next section presents background on Algebraic Dynamic Programming. Section 3 presents the GPU parallelization of ADP. Section 4 gives some results and discussion: depending on the application, we get speedups ranging from 6.1× to 25.8× on a Nvidia GTX 280.

Algebraic Dynamic Programming

We briefly introduce the Algebraic dynamic programming (ADP) methodology on a simple Nussinov type RNA secondary structure prediction problem, the maximization of the number of base pairs [START_REF] Nussinov | Algorithms for loop matchings[END_REF]. Section 4 reports results for several other applications. See [START_REF] Giegerich | A discipline of dynamic programming over sequence data[END_REF] for a complete presentation of the ADP method. gucaugcaguguca (...)((...)).. When designing a dynamic programming algorithm in algebraic style, we need to specify four constituents: the input alphabet, the search space, the scoring of the candidates, and the objective function.

Alphabet. The input RNA sequence is a string over the alphabet A = {a, c, g, u}.

Search space. Given an input sequence w ∈ A * , the search space is the set of all possible secondary structures the sequence w can form. In the ADP terminology, the elements of the search space for a given input sequence are called candidates. Figure 2 gives example of candidates for RNA folding. This tree representation of candidates exists for any application of dynamic programming [START_REF] Steffen | Versatile and declarative dynamic programming using pair algebras[END_REF]. To describe the candidates, the ADP methodology uses the notion of tree grammar. Figure 3 shows the grammar nussinov78, origin of our two example trees. For each sequence w ∈ A * , the grammar defines a search space P G (w) that is the set of all parses of the sequence w for G. Scoring and objective. Given an element of the search space as a tree t ∈ P G (w), we need to score this element. In our example we are only interested in counting base pairs, so scoring is very simple: the score of a tree is the number of pairnodes in t. For the two candidates of Figure 2 we obtain scores of 3 (t 1 ) and 4 (t 2 ). Moreover, we need to choose one or several solutions from the pool of candidates. For this purpose we add an objective function h which chooses one or more elements from a list of candidate scores.

= s + 1 split(s,s') = s + s' h([]) = [] h([s1, . . . , sr]) = [max 1≤i≤r si]
Scoring schemes with objective functions are called evaluation algebras in ADP. The above example is the evaluation algebra bpmax for maximizing the number of base pairs. The flexibility of the algebraic approach lies in the fact that we don't have to stop with definition of one algebra: simply define another algebra and get other results for the same search space. We use the notation E(t) to indicate the value obtained from t under evaluation with algebra E. All that is left to do is to evaluate the candidates in a given algebra, and make our choice via the objective function h. For example, candidates t 1 and t 2 of Figure 2 are evaluated by algebra bpmax, with h(bpmax(t 1 ), bpmax(t 2 )) = [max(3, 4)] = [START_REF] Charalambous | Initial experiences porting a bioinformatics application to a graphics processor[END_REF].

This example was fairly simple: complete RNA folding algorithms are typically based on energy minimization, and include energies of stacking regions (or helices), bulge loops, internal loops, hairpin loops and multiple loops. Figure 4 shows an excerpt of the real RNAfold.lhs grammar that includes the full Turner model [START_REF] Mathews | Expanded sequence dependence of thermodynamic parameters improves prediction of RNA secondary structure[END_REF]. The grammar can be read as a standard context-free grammar. The operator ~~~connects succeeding symbols and the operator ||| divides alternative productions for a non-terminal. The symbol <<< denotes application of an algebra function and ... denotes application of the evaluation function h. Finally, the operator with denotes the use of a filter function, that means that (base ~~~closed ~~~base) 'with' basepairing is only successful if the two bases can form a base pair.

Automatic Parallelization of ADP

Principle. A compiler that translates ADP programs into C was previously developed [START_REF] Giegerich | Challenges in the compilation of a domain specific language for dynamic programming[END_REF]. This task includes some advanced optimization techniques, see [START_REF] Steffen | Compiling a Domain Specific Language for Dynamic Programming[END_REF] for a detailed overview. With the option -cuda, the compiler is now switched into the CUDA code generation mode. The compiler uses the same backend both for CPU and GPU code generation and only differs in the following parts:

1. The dynamic programming tables need to be stored both on the host and on the global memory of the GPU. The compiler generates code to synchronize these tables.

2. For each dynamic programming table, the compiler generates a calculation function. This is the same function both for CPU and GPU mode, so the only change is that it is declared to be executed as a GPU kernel. Figure 5 shows the CUDA code for the dynamic programming main loop in the RNA secondary structure prediction program. The kernel function, calc all, contains the calls for the calculation of the six dynamic programming tables.

3. In CPU mode, all table elements are calculated sequentially with increasing subword length. This order of computation has to be changed to enable parallelization. For the RNA secondary structure prediction program, the calculation of a table element (i, j) depends only on the elements in the triangle in the lower left (see Figure 6, on the left). So all elements in one diagonal can be calculated in parallel. The whole dynamic programming table is then calculated in a loop over all diagonals (see Figure 5). This approach can be generalized to all dynamic programming algorithms over sequence data: in all ADP grammars, all results are combined from results of shorter subsequences. Therefore, the calculation of a table element (i, j) depends only on results that lie between the indices (i, j).

All these changes are done automatically by the compiler and do not require any changes to the ADP grammar. The number of blocks and threads used for the calculation can be configured as a parameter at runtime.

__global__ static void calc_all(int diag, int n) { int i = blockIdx.x*blockDim.x+threadIdx.x; int j = i + diag; if ((i <= n) && (j <= n)) {
calc_closed(i, j); calc_initstem(i, j); calc_struct(i, j); calc_block(i, j); calc_comps(i, j); } } static void mainloop() { for (int diag=0; diag<=n; diag++) { (...) calc_all <<< grid, threads >>> (diag, n); } } Fig. 5. Kernel and main CUDA code for the dynamic programming main loop corresponding to the grammar shown on Figure 4. Note that each kernel thread also computes inner loops for the folding calculations. Actual codes are available on the ADP website (http://bibiserv.techfak.uni-bielefeld.de/adp/cuda.html).

Window mode. It does not make any sense to fold a complete genome as a single RNA molecule. This remark is the same for other applications: for example, a thermodynamic matcher [START_REF] Reeder | Locomotif: From graphical motif description to RNA motif search[END_REF] looks for some small sub-sequences (50 to some Fig. 6. Left: data dependencies for RNA secondary structure prediction. The computation of the table element (i, j) needs the O((j -i) 2 ) elements in the underlying triangle. Right: window mode. With a large genome (of size n), we just need to fold sequences on small windows (of size w).

hundreds bases) matching a given structural pattern in a large sequence. In those applications, we need to compute only some diagonals above the main diagonal (Figure 6, on the right). The option -cudaw sets the ADP compiler in window mode. Whereas the CPU version sequentially computes all windows, the GPU version loads a large sequence into the GPU and launches a large number of independent threads, thus increasing the parallelism. This is also done automatically by the compiler and does not require any changes to the source program.

Results

Table 1 shows the results on three different applications with RNA sequences: RNAfold (see previous section), pknotsRG (detection of pseudo-knots [START_REF] Reeder | Design, implementation and evaluation of a practical pseudoknot folding algorithm based on thermodynamics[END_REF]), and a tRNA thermodynamic matcher. The program pknotsRG predicts RNA secondary structures including a restricted class of pseudoknots. The thermodynamic matcher was created by the graphical tool Locomotif [START_REF] Reeder | Locomotif: From graphical motif description to RNA motif search[END_REF]. Practical 9.9×, 14.5× and 6.1× speedups are obtained on those real applications with a GTX 280. In these speedups, the main bottlenecks are in memory transfers, as only the global memory of the GPU is used.

In the original RNAfold.lhs grammar, a part of production calc closed is in fact computed for only 6 out of the 16 possible basepairs (filter 'with' stackpairing on Figure 4). This brings a large divergence between the threads and breaks the GPU SIMD model. To confirm this fact, we tested a special version of this grammar, RNAfold-bp.lhs, that computes for every basepair the full recurrence equations (penalizing non-pairs): the speedup with the GTX 280 is almost doubled. This indicates that a similar speedup would result for the calculation of stochastic grammars, since here arbitrary base pairs are considered. It should be noted that our speedups are lower than the best possible ones. For example, Rizk and Lavenier [START_REF] Rizk | GPU accelerated RNA folding algorithm[END_REF] developed an optimized GPU RNAfold implementation: in particular, they pack together the 6/16 computations corresponding to the production calc closed. They obtain a 17× speedup on a GTX 280 against one core of a 2.66 GHz Xeon (applied on a whole sequence, without window mode), whereas our best speedup without window mode is only 2.8× (results not shown). However, as our approach is generic, it can be applied on several algorithms with few efforts to the user.

On pknotsRG, runs with w = 160 get a little smaller speedup than with w = 80. As w is fixed, this does not limit the scalability of our approach: the input data size, n, can always grow with the same speedup.

Current limitations. Whereas grammars involving several sequences can be encoded in the ADP formalism, the ADP compiler now only works for one input sequence. Removing this limitation would allow to study other dynamic programming problems, as for example Smith-Waterman sequence alignment or RNA co-folding [START_REF] Giegerich | A discipline of dynamic programming over sequence data[END_REF]. Finally, for some grammars (including the tRNA matcher), the ADP automatic table design generates some recursive functions, and those functions cannot be compiled with the CUDA libraries (there is no stack on the current cards). This automatic table design is removed through omitting the -cto option, but, in this case, the grammar should specify precisely which symbols of the grammar are to be "tabulated".

We implemented a parallel GPU CUDA backend for the ADP compiler, which works out-of-the-box for several grammars dealing with RNA sequences. The new ADP compiler and some example codes are available on the ADP website (http://bibiserv.techfak.uni-bielefeld.de/adp/cuda.html). Our approach is generic and requires few efforts to the user, even if the speedups are not the best ones that could be obtained by manually optimized implementations. We plan to remove the limits explained above. Other perspectives include the following points.

Shared memory. The ADP compiler could be improved to better use the memory hierarchy of the card. In the Nvidia architecture, a 16 KB shared memory is available for the threads in the same block. This local memory is very fast and should be used to maximize the efficiency, for example in storing portions of some dynamic programming tables. This memory is not used in our current implementation. Of course, the best usage of the shared memory depends on the application: for now, we did not find a generic way to determine from the grammar this best usage. Some hints given in the grammar file could indicate to the compiler which dynamic programming tables should be handled in this way.

Static evaluation of grammars. We plan to test other grammars, in bioinformatics as in other domains. Which grammars are efficient to parallelize, and why?

Other targets. We plan to test the ADP methodology on other manycore architectures, in particular through the new OpenCL standard [START_REF] Giraud | OpenCL 1.0 specification[END_REF]. Again, the fact that the ADP methodology is generic allows to write portable solutions.

Fig. 1 .

 1 Fig. 1. ADP workflow. The goal of this study is to conceive and implement an automatic GPU parallelization (bottom).

Fig. 2 .

 2 Fig.2. Two candidates in the search space for the best secondary structure for the sequence gucaugcaguguca.

Fig. 3 .

 3 Fig. 3. Tree grammar nussinov78 consists of one production with four alternatives. Symbol Z denotes the axiom of the grammar.

  bpmax = (nil, right, pair, split, h) where nil(s) = 0 right(s,b) = s pair(a,s,b)

Table 1 .

 1 Time (real times, in seconds) for executing different ADP grammars. CPU versions are compiled with adpc, and executed on a 2.4 GHz Core2 processor (PC1, 1 core used). and on a 3.0 GHz Xeon X5472 processor (PC2, 1 core used). GPU CUDA versions are compiled adpc -cudaw, and executed on a GeForce 8800 (PC1) and on a GTX 280 (PC2). Because of its increased number of cores and of its better handling of uncoalesced memory loads, the GTX 280 gives better speedups than the GeForce 8880. Note that the performance of the CPU does not impact the times reported for the GPU versions. For example, for RNAfold, the 19.22s for the PC1 GPU include only 0.20s of non-kernel computations, mainly for traceback in the DP matrix. Tests on RNAfold and tRNA-matcher were done on the 160 kbp genome of Candidatus Carsonella ruddii (Genbank reference NC 008512). Tests on pknotsRG were done on the first 20 kbp of the same genome.

	Grammar,	PC1		PC2
	window size,	Core2 + GeForce 8800	Xeon + GTX 280
	time complexity	CPU GPU speedup CPU GPU speedup
	RNAfold-bp.lhs	-w 80 O(w 2 n) 176.09 19.22	9.1× 133.77 5.18	25.8×
	RNAfold.lhs	-w 80 O(w 2 n) 43.43 8.08	5.4×	35.57 3.59	9.9×
	tRNA-matcher.lhs -w 100 O(w 2 n) 52.46 6.76	7.8×	43.60 3.01	14.5×
	pknotRG.lhs	-w 80 O(w 3 n) 26.82 10.64	2.5×	23.54 3.25	7.2×
	pknotRG.lhs	-w 160 O(w 3 n) 188.68 87.65	2.2× 166.27 27.22	6.1×

Acknowledgements

Part of this research was done during P. Steffen's stay in Université Lille 1. This research was carried through the "NVIDIA Professor Partnership" program.


The problem of diagnosis with respect to physical parameters for changes in structures

George V. Moustakides


N Keeley

R S Mackintosh

C Beck

Breakup Coupling Effects on Near-Barrier 6 Li, 7 Be and 8 B + 58 Ni Elastic Scattering Compared

New data for near-barrier 6 Li, 7 Be and 8 B + 58 Ni elastic scattering enable a comparison of breakup coupling effects for these loosely-bound projectiles. Coupled Discretised Continuum Channels (CDCC) calculations suggest that the large total reaction cross sections for 8 B + 58 Ni are dominated by breakup at near-barrier energies, unlike 6 Li and 7 Be where breakup makes a small contribution. In spite of this, the CDCC calculations show a small coupling influence due to breakup for 8 B, in contrast to the situation for 6 Li and 7 Be. An examination of the S matrices gives a clue to this counter-intuitive behaviour.

INTRODUCTION

Recent data [1] for near-barrier 6 Li, 7 Be and 8 B + 58 Ni elastic scattering allow some interesting comparisons for these weakly-bound nuclei. Optical model fits find much larger total reaction cross sections (σ R ) for 8 B than for 6 Li or 7 Be, even when "reduced" [2]; while the reduced σ R for other weakly-bound projectiles lie on a universal curve, those for 8 B and 6 He are significantly larger [1]. The low 8 B → 7 Be + p breakup threshold (0.1375 MeV) suggests a dominant contribution to the direct part of σ R . This is not automatic: for 6 He with an α + 2n breakup threshold of 0.973 MeV, 1n-and 2n-stripping are the main contributors to σ R at near-barrier energies. However, the weakly-bound proton in 8 B experiences Coulomb barrier and charge polarisation effects tending to suppress transfer.

CDCC calculations [3] find that breakup does dominate the direct component of σ R for 8 B: as the cross sections are large -of the order of 100 mb or more -one might expect an equally important coupling effect on the elastic scattering angular distribution. However, this is not the case [3]. We thus have an apparent paradox: 6 Li, with a relatively small breakup cross section, exhibits an important breakup coupling effect on the elastic scattering (see e.g. [4]) whereas 8 B, with a large breakup cross section, shows only a modest coupling effect. A comparison of S matrices obtained from CDCC calculations for 6 Li, 7 Be and 8 B + 58 Ni provides a clue to this behaviour. Preliminary dynamic polarisation potentials (DPPs) are also presented.

CALCULATIONS

Calculations were performed with the code FRESCO [5]: only a brief outline is given here. The 6 Li, 7 Be and 8 B nuclei were modelled as α + d, α + 3 He and 7 Be + p clusters, respectively. The 7 Be core was treated as inert but its non-zero spin was retained. Interaction potentials were obtained by Watanabe-type folding of global optical potentials, with a 6 Li potential as surrogate for 7 Be, the well-depths being adjusted to give the best fit to the data. The 6 Li and 7 Be calculations were similar to those in [4] and [6], but with finer continuum binning for 7 Be. The 8 B calculations included couplings to the L = 0, 1, 2 and 3 continuum and the 0.774 MeV 1 + and 2.32 MeV 3 + resonances. Good fits to all the data were obtained. Due to lack of space we show only results for the same values of E c.m. -V B for each system, where V B is the nominal Coulomb barrier, taking as our "benchmark" the 8 B data at E lab = 29.26 MeV. This procedure yields values of E lab = 19.04 and 24.12 MeV for 6 Li and 7 Be + 58 Ni, respectively. In this way effects due to differences in projectile charge should be minimised.

Results are presented in Fig. 1: the coupling effect is much stronger for 6 Li and 7 Be, with 6 Li → α + d and 7 Be → α + 3 He breakup thresholds of 1.47 and 1.59 MeV, respectively, an order of magnitude larger than the 8 B → 7 Be + p threshold. The 6 Li → α + d process has the additional peculiarity that it cannot proceed via dipole breakup. If we include population of the bound 1/2 -state in 7 Be (considering breakup as an inelastic excitation) the total breakup cross sections for both 6 Li and 7 Be are about a factor of three smaller than for 8 B. To obtain a clue to this apparent paradox, we show in Fig. 2 the modulus and argument of the J-weighted S matrices [7] obtained from full and no-coupling calculations. The coupling effect on |S| is almost negligible for 8 B and largest for 6 Li, but qualitatively similar for all three nuclei: a decrease of |S| at small L and an increase at large L. By contrast, for arg(S) the coupling effect is greatest for 8 B, smallest for 7 Be and intermediate for 6 Li.

DISCUSSION

For protons and other light particles, changes in |S| correspond to changes in the imaginary part of the potential, while changes in arg(S) correspond to changes in the real part. While this simple picture is not so clear-cut in the presence of strong absorption (as here) it provides a useful guide. Thus, the coupling effect on |S| suggests reduced absorption at small L, switching to increased absorption at large L. The effect on arg(S) suggests repulsion at small L and attraction at large L. These effects are qualitatively similar for all three nuclei. The fact that the coupling effect on both the elastic scattering and |S| is so small for 8 B suggests that, paradoxical as it may seem for a coupling producing such a large cross section, its effective imaginary potential is small.

DPPs may be obtained by inversion of the S matrix, see e.g. [8]. In Fig. 3 we show the results of such a procedure for 8 B and 7 Be. Those for 7 Be are preliminary; we expect the final DPPs to be somewhat smoother. While the DPPs are qualitatively similar, short-range repulsion and long-range attraction combined with surface absorption (this behaviour seems to be universal, see e.g. [9]), the details are very different. The small imaginary DPP for 8 B is particularly striking, confirming the conclusions inferred from the S matrices. The surface repulsion for 8 B is also much smaller than for 7 Be, although for radii larger than about 9 fm it is significantly larger than for 7 Be, having a longer, more repulsive tail. Our results show that a large cross section is no guarantee of a large coupling effect. The S matrices and DPPs shed some light on this, but it remains to be explained at a more fundamental level.

Figure 1 .

 1 Figure 1. CDCC calculations for 8 B, 7 Be and 6 Li + 58 Ni at E lab = 29.26, 24.12 and 19.04 MeV. Solid and dashed curves denote full and no-coupling results, respectively.

Figure 2 .

 2 Figure 2. |S| and arg(S) from CDCC calculations for 8 B, 7 Be and 6 Li + 58 Ni at E lab = 29.26, 24.12 and 19.04 MeV. Solid and dashed curves denote full and no-coupling results, respectively.

Figure 3 .

 3 Figure 3. DPPs from CDCC calculations for 8 B (solid curves) and 7 Be (dotted curves).


Hugo Gimbert

Randomized Strategies are Useless in Markov Decision Processes

We show that in a Markov decision process with arbitrary payoff mapping, restricting the set of behavioral strategies from randomized to deterministic does not influence the value of the game nor the existence of almost-surely or positively winning strategies. As a corollary, we get similar results for Markov decision processes with partial observation.

Definitions

We use the following notations throughout the paper. Let S be a countable set. The set of finite (resp. infinite) sequences on S is denoted S * (resp. S ω ). and S ω denotes the set of infinite sequences u ∈ S N . A probability distribution on S is a function δ : S → R such that ∀s ∈ S, 0 ≤ δ(s) ≤ 1 and s∈S δ(s) = 1. The set of probability distributions on S is denoted D(S). Definition 1 (Markov Decision Processes). A Markov decision process M = (S, A, (A(s)) s∈S , p) is composed of a countable set of states S, a countable set of actions A, for each state s ∈ S, a set A(s) ⊆ A of actions available in s, and transition probabilities p : S × A → D(S) .

In the sequel, we only consider Markov decision processes with finitely many states and actions.

An infinite history in M is an infinite sequence in (SA) ω . A finite history in M is a finite sequence in S(AS) * . The first state of an history is called its source, the last state of a finite history is called its target. A strategy in A is a function σ : S(AS) * → D(A) such that for any finite history s 0 a 1 • • • s n , and every action a ∈ A, (σ(s 0 a 1

• • • s n )(a) > 0) =⇒ (a ∈ A(s n )).
We are especially interested in strategies of the following kind.

Definition 2 (Deterministic strategies). A strategy σ is deterministic if for every finite history h and action a, (σ(h)(a) > 0) ⇐⇒ (σ(h)(a) = 1).

Given a strategy σ and an initial state s ∈ S, the set of infinite histories with source s is naturally equipped with a σ-field and a probability measure denoted P σ s . Given a finite history h and an action a, the set of infinite histories in h(AS) ω and ha(SA) ω are cylinders that we abusively denote h and ha. The σ-field is the one generated by cylinders and P σ s is the unique probability measure on the set of infinite histories with source s such that for every finite history h with target t, for every action a ∈ A and for every state r,

P σ s (ha | h) = σ(h)(a) , (1)
P σ s (har | ha) = p(r|t, a) . (2)
For n ∈ N, we denote S n and A n the random variables S n (s

0 a 1 s 1 • • • ) = s n and A n (s 0 a 1 s 1 • • • ) = a n .
Some strategies are better than other ones, this is measured by mean of a payoff function. Every Markov decision process comes with a bounded and measurable function f : (SA) ω → R, called the payoff function, which associates with each infinite history h a payoff f (h).

Definition 3 (Values and guaranteed values).

Let M be a Markov decision process with a bounded measurable payoff function f : (SA) ω → R. The expected payoff associated with an initial state s and a strategy σ is the expected value of

f under P σ s , denoted E σ s [f ].
2 Randomized strategies are useless

Randomizing his own behaviour is useless when there is no adversary to fool. This is the intuitive interpretation of the following theorem:

Theorem 4. Let M be a Markov decision process with a bounded measurable payoff function f : (SA) ω → R, x ∈ R and s a state of M. Suppose that for every deterministic strategy σ, E σ s [f ] ≤ x. Then the same holds for every randomized strategy σ.

Proof. For simplifying the notations, suppose that for every state s there are only two available actions 0, 1 and for every action a ∈ {0, 1} there are only two successor states L(s, a) and R(s, a) distinct and chosen with equal probability 1 2 . Let σ be a strategy and s an initial state. We define a mapping

f s,σ : {L, R} ω × [0, 1] ω → (SA) ω
that will be used for proving that P σ s is a product measure. With every infinite word u ∈ {L, R} ω and every sequence of real numbers x = (x n ) n∈N ∈ [0, 1] ω between 0 and 1 we associate the unique infinite play

f s,σ (u, x) ∈ (SA) ω = s 0 a 1 s 1 • • • such that s 0 = s, for every n ∈ N if u n = L then s n+1 = L(s n , a n+1 ) otherwise s n+1 = R(s n , a n+1 ) and for every n ∈ N, if σ(s 0 a 1 • • • s n )(0) ≥ x n then a n+1 = 0 otherwise a n+1 = 1.
We equip {L, R} ω with the σ-field generated by cylinders and the natural head/tail probability measure denoted µ 1 . We equip [0, 1] ω with the σ-field generated by cylinders I 0 ×I 1 • • •×I n ×[0, 1] ω where I 1 , I 2 , . . . , I n are intervals of [0, 1], and the associated product of Lebesgue measures denoted µ 2 .

Then P σ s is the image by f s,σ of the product of measures µ 1 and µ 2 , i.e. for every measurable set of infinite plays A,

P σ s (A) = (µ 1 × µ 2 )(f -1 s,σ (A)) . (3)
This holds for cylinders hence for every measurable A. Now:

E σ s [f ] = p∈(SA) ω f (p)dP σ s = (u,x)∈{L,R} ω ×[0,1] ω f (f σ,s (u, x))d(µ 1 × µ 2 ) = x∈[0,1] ω u∈{L,R} ω f (f σ,s (u, x))dµ 1 dµ 2
where the first equality is by definition of E σ s [f ], the second equality is a basic property of image measures and the third equality is Fubbini's theorem, that we can apply since f is bounded and the measures are probability measures.

Once x is fixed, the behaviour of strategy σ is deteministic. Formally, for every x ∈ [0, 1] let σ x be the deterministic strategy defined by σ

x (s 0 a 1 • • • s n ) = 0 if and only if σ(s 0 a 1 • • • s n )(0) ≥ x n . Then for every y ∈]0, 1[ ω and u ∈ {L, R} ω , f σx,s (u, y) = f σ,s (u, x) hence: E σx s [f ] = u∈{L,R} ω f (f σ,s (u, x))dµ 1 ,
and finally:

E σ s [f ] = x∈[0,1] ω E σx s [f ] dµ 2 ,
hence the theorem, since for every x, strategy σ x is deterministic.

Applications

We provide an extension of Theorem 4 to Markov decision processes with partial observation.

A Markov decision process with partial observation is similar to a Markov decision process except every state s is labelled with a color col(s) and strategies should depend only on the sequence of colors. Formally, a strategy is said to be observational if for every finite plays s 0

• • • s n and t 0 • • • t n , if col(s 0 • • • s n ) = col(t 0 • • • t n ) then σ(s 0 • • • s n ) = σ(t 0 • • • t n ).
Corollary 5. Let M be a Markov decision process with a bounded measurable payoff function f : (SA) ω → R, x ∈ R and s a state of M. Suppose that for every deterministic observational strategy σ, E σ s [f ] ≤ x. Then the same holds for every randomized observational strategy σ.

Proof. Fix an initial state s. Consider the Markov decision process whose state space is the set of finite sequences a 0 c 0 a 1 • • • a n c n ∈ (AC) * of colors interleaved with actions. The initial state is the empty sequence. From state a 0 c 0 a 1 c 1 • • • a n c n , playing action a leads to state a 0 c 0 a 1 c 1 • • • a n c n ac with probability:

P σ s (A n+1 = a, col(S n+1 ) = c | A 0 C 0 A 1 • • • A n C n = a 0 c 0 a 1 • • • a n c n ) ,
and the payoff associated with an infinite play is defined by:

g(a 0 c 0 a 1 c 1 • • • ) = E σ s [f | A 0 C 0 A 1 C 1 • • • = a 0 c 0 a 1 c 1 • • • ] ,

where in both definitions σ is any deterministic strategy such that for every

The state space of this new Markov decision process is countable therefore we can apply Theorem 4 to it, which immediately gives us the result.


Boyd Michailovsky

On the etymology of Nepali jhyāl 1

Nepali jaskel/jaskelo 'side door, window' is proposed as the immediate source of jhyāl 'window', whose etymology has long been disputed. The origin of jaskelo remains unknown.

Introduction

Nepali jhyāl 'window' appears to be a Nepalese lexical innovation, surrounded as it is, in neighboring Indo-Aryan languages (and even in the Nepali of Dorjiling) by khirkī /khiḍkī . Turner proposed no etymology in his great Nepali dictionary (1965[1931]), and the word does not appear in his Indo-Aryan comparative dictionary (1966,1985).

Other lexicographers have tried to fill this gap. Bal Candra Sharma, in the first Royal Nepal Academy dictionary (1962/1963), declared that the word was of Newar origin (cf. modern Newar jhyāː ), a hypothesis which gave rise to lively controversy. In a characteristically polemical article, "jhyāl kasko [whose is window?] ", Mahananda Sapkota (V.S. 2027 = 1970 C.E.) refuted the idea of Newar origin and identified jhyāl as a descendant of Sanskrit jāla 'net'. This etymology appeared in the second Royal Academy Dictionary (Pokharel 1983/4).

On the semantic side, Sapkota's argument for jāla as the source for 'window' depends on the topological resemblance between lattice-windows, or even woven bamboo shutters, and nets. There is support for this connection in Sanskrit itself; Monier-Williams (1970:419-420) cites 'lattice-window' as a secondary meaning of jāla. Turner, in his comparative dictionary (1966, no. 5213), cites one or two instances where modern descendants of Sanskrit jāla have the secondary meaning 'window'. But he was not moved to include Nepali jhyāl among these.

One point raised by Sapkota against Newar origin was quickly refuted by the Newar scholar Indra Mali. Newar origin does not imply the unmotivated addition of final l by Nepali speakers. Modern Newar jhyāː is a regular reflex of older Newari jhyāl, which itself continues an still older form jhyār [START_REF] Mali | jhyāl' śabdako ek adhyayan[END_REF]). 2 My purpose is not to review this old controversy in detail, but to introduce a new piece of evidence. In the end, the best argument against an etymology is a more convincing one.

The controversy over jhyāl was fresh in my mind when, in 1972, in a Hayu village in Ramechhap, among houses which had no windows, I learned the Nepali word jaskelo, used for the small side door on the downhill side of the house. Turner (1965:211) has the following: jaskel, in jaskel dailo a small gate at the side of the house, wicket-gate. jaskelo, s. A small window. [cf. jaskel.] (I have since heard jaskelo used for small windows in Achham.) This, I thought, was the origin of jhyāl .

The word jaskel is not attested in other Indo-Aryan languages or in Newar.

Semantically, the relationship between 'door' and 'window' is a classic one. As Buck (1949:469) remarks concerning Indo-2 I do not know the precise phonetic significance of the final -r in this written form.

European generally, "Words for 'window' are connected with those for 'door', 'light', 'wind', and 'eye'." Nepali khirkī and widespread related words in India (Turner 1966 no. 3770) are examples.

Phonetically, I do not know of any exact parallel to the change from jaskelo to jhyāl, but it appears less difficult than the complexification of jāl. 3 In particular, the e of jaskel is a plausible source for the yā of jhyāl.

Paper presented at the

th Annual Conference of Linguistic Society ofNepal, Kathmandu, November, 2005.

Accepting jaskel as the source of jhyāl leaves a number of questions unanswered, but it does appear to rule out any role for jāla 'net'. The origin of jaskel itself, and the date and the authors of its transformation to jhyāl are questions that remain for further research.