O Bénichou

V Calvez

N Meunier

R Voituriez

Front acceleration by dynamic selection in Fisher population waves

Keywords:

We introduce a minimal model of population range expansion in which the phenotypes of individuals present no selective advantage and differ only in their diffusion rate. We show that such neutral phenotypic variability can yield alone phenotype segregation at the front edge even in absence of genetic noise, and significantly impact the dynamical properties of the expansion wave. We present an exact asymptotic traveling wave solution and show analytically that phenotype segregation accelerates the front propagation. The results are compatible with field observations such as invasions of cane toads in Australia or bush crickets in Britain.

I. INTRODUCTION

The combination of random mutations and natural selection plays a key role in evolution. In the case of large well-mixed populations, it is found typically that mutations that eventually become dominant must present a selective advantage (i-e a larger effective growth rate) over the wild type [START_REF] Crow | An Introduction to Population Genetics Theory[END_REF]. More recently, it appeared that the case of expanding populations, which are a common phenomenon in biology, can lead to strikingly different behaviors. It was observed that in population range expansions even neutral mutations, i-e mutations that present no selective advantage, can prevail at least transiently due to large fluctuations at the edge of the population traveling front, and lead to spectacular gene segregation phenomena [START_REF] Edmonds | [END_REF][3][4]. Such population waves are exemplified by the hypothesized migration of humans from Africa [5], or invasions of species such as cane toads in Australia [6] or bush crickets in Britain [7], and can now be observed at a smaller scale on the example of migrations of microorganisms in a Petri dish [3]. Following these observations, numerical models have been developed in the ecology community to analyze quantitatively range expansions [8][9][10][11][12][13]; these studies give further support to the mechanism of phenotype selection, which was dubbed spatial sorting in [14].

The standard mathematical models for population waves are generally based on the Fisher-Kolmogorov (FKPP) equation [15,16], which is the simplest non linear equation that combines diffusion and growth. Beyond population genetics, this equation has been widely used for example in ecology, epidemiology [START_REF] Murray | Mathematical Biology[END_REF] or chemistry [START_REF] Van Saarloos | [END_REF]. Fisher waves are still actively investigated from the theoretical point of view and despite this effort, exact results, which mostly concern the propagation speed [19,20], remain elusive especially in the presence of noise or in space dimensions larger than one.

In this letter, we introduce a minimal model of popu-lation dynamics based on a modified FKPP equation in which mutations are neutral and affect only the mobility properties of individuals. Our approach is inspired by the observations of range expansions of cane toads in Australia [6] or bush crickets in Britain [7], as well as several numerical studies [8][9][10][11][12][13][14] which highlighted the following features: (i) The leading part of the expanding front is mostly populated by the fastest phenotypes, that is longer-legged individuals in the case of cane toads [6] and longer-winged individuals in the case of bush crickets. (ii)

The proportion of fast phenotypes is the highest in locations of new arrivals and then declines over time. (iii) The speed of the population wave increases with time. Our analysis is based on a modified FKPP equation with an additional mutation term that allows for variations of the diffusion rate of individuals, for which we present an exact asymptotic traveling wave solution and calculate analytically the front speed. We show that such phenotypic variability, when taken into account, can yield alone the segregation of the fastest phenotypes at the front edge in absence of genetic noise, and significantly accelerate the propagation speed. Importantly, the results reproduce qualitatively the observations (i)-(iii) above made on the basis of field experiments and numerical simulations. The novelty of our approach lies on the modified FKPP equation that we propose in Eq.( 1), which enables an analytical derivation of the main properties of phenotypic segregation in range expansions. We analyze in this paper this equation in its minimal form, keeping in mind that it could be modified to account more faithfully for real systems. Given the wide spectrum of situations described by the FKPP equation, potential applications for example in the context of chemical kinetics are also expected. The leading edge of the front is mostly populated by fastest individuals (with diffusion rate D * > D0) that are selected dynamically, while the population behind the front is characterized by the steady state value D0.

II. MODEL

We consider the following 1d model (see Fig. 1). We assume that each individual of a growing population is characterized by a diffusion coefficient D, which is subject to mutations, and we denote by n(D, x, t) the density of such individuals of phenotype D at position x and time t. We will also make use of the marginal density N (x, t) = ∞ 0 n(D, x, t)dD and the probability distribution of the trait D at (x, t) defined by P x (D, t) ≡ n(D, x, t)/N (x, t). We assume that mutations, which only affect D, are neutral and do not affect the population growth rate s. Mutations are described phenomenologically by a density current j[n(D, x, t)] in the phenotype space (an explicit example is given below). The dynamics is then given by the following modified FKPP equation :

∂ t n = D∂ 2 x n + sn(1 -N ) -∂ D j. (1)
The diffusion term on the right-hand side accounts for the mobility of the individuals, which is effectively modelled by a diffusion process of coefficient D. The logistic term sn(1 -N ) accounts for the birth and death of individuals, where the growth rate s is independent of the phenotype D. Note that here we omit number fluctuations of the birth and death process. Last, the divergence of the density current -∂ D j accounts for the dynamics of mutations.

To get insight into the effect of the mutation term, we first consider the case of a well-mixed population in a bounded range (for example x ∈ [0, 1]). The stationary state (see also [21] for the analysis of a related problem) is then given by N (x, t) = 1 for all x, with P x (D, t) ≡ P 0 (D) independent of x and defined by j = 0. The function P 0 (D) therefore characterizes the steadystate distribution of the trait D in a well-mixed population. We will consider the representative case of a distribution P 0 that is centered around a mean value D 0 with a typical width ∆D. Without loss of generality, a simple choice for the density current is then

j = -α∂ D n(x, D, t) (2)
where α is the mutation rate, together with reflecting boundary conditions (j = 0) at D = D -≡ D 0 -∆D/2 and D = D + ≡ D 0 + ∆D/2. With this choice of j, the distribution P 0 is then uniform over the range [D -, D + ] and the different phenotypes have no selective advantage. Another example of interest is given by

j = -α [β(D -D 0 )n(x, D, t) + ∂ D n(x, D, t)] , (3)
where α is the mutation rate and the extra term β accounts here for a selective advantage of phenotype D 0 . This choice, together with the zero flux condition at D = D -and D = D + , yields P 0 (D) ∝ exp -β(D -D 0 ) 2 /2 : at stationary state in a bounded range, the trait D 0 is now favored, as is the case in many real situations. Alternative choices are possible to account more faithfully for field experiments, and would lead qualitatively to similar results. We will focus below on the minimal choice of Eq.( 2), which captures the main features of the model and keep calculations analytically tractable. We now turn to the case of a range expansion and assume that at t = 0 the population is seeded with the uniform distribution P 0 in a bounded interval of the x axis. We show that despite the fact that mutations are neutral, a traveling wave solution characterized by a non homogeneous distribution of D, which we calculate in the regime u ≡ x -v * t 1, emerges. We determine the velocity v * and show that it is larger than the velocity v 0 = 2 √ sD 0 expected from the classical FKPP model with an homogeneous population of diffusion rate D 0 .

Let us introduce u ≡ x-vt where v is to be determined, and look for a traveling front so that N (u → -∞) = 1 and N (u → +∞) = 0. Notations are kept unchanged in the co-moving frame of velocity v and all functions of x and t are assumed to depend on u only. Eq. ( 1) then rewrites

0 = v∂ u n + D∂ 2 u n + sn(1 -N ) + α∂ 2 D n. (4)
The key point of the following analysis is that one can show self-consistently that in the limit u → ∞, one has P x (D, t) ≡ P u (D) → P * (D), where P * (D) is independent of u. The dependences on u and D in n(D, u) are therefore factorized at the leading edge of the front and one can write in this limit: n(D, u) P * (D)N (u). This factorized form of the solution can be shown to hold in a mathematically rigorous framework which will be presented elsewhere. We give here the main ingredients of derivation. In the regime u → ∞, n is small and Eq.(4) can be linearized, which suggests the following ansatz n(D, u) = P * (D) exp(-γu), (5) with γ > 0. Eq.( 4) then rewrites where it should be noted that P * (D) depends on γ. The speed v of a front of prescribed profile characterized by γ is then obtained by integration over D and follows:

0 = (-γv + Dγ 2 + s)P * (D) + α∂ 2 D P * (D), (6)
v = D * γ + s/γ, (7)
where D * , which depends on γ, is determined below. Eq.( 6), completed by ∂ D P * | D=D-= ∂ D P * | D=D+ = 0 for the choice (2) of the density current j, explicitly defines the function P * and shows that it is independent of u in the limit u → ∞, therefore proving that the factorized form of Eq. ( 5) holds true. We stress that this argument holds for any density current j, which would only modify the ordinary differential equation ( 6) defining P * . In the explicit example of Eq. ( 2), the existence of a non zero solution for Eq. ( 6) implies that :

Ai[1, µ(D * -D -)]Bi[1, µ(D * -D + )] (8) -Bi[1, µ(D * -D -)]Ai[1, µ(D * -D + )] = 0
where µ = (γ 2 /α) 1/3 and Ai and Bi denote Airy functions. This last equation, which can not be solved analytically to the best of our knowledge, implicitly defines D * as a function of γ, and completes the definition of v in Eq.( 7). The distribution P * is then given by

ZP * = Ai[µ(D * -D)]- Ai[1, µ(D * -D -)] Bi[1, µ(D * -D -)] Bi[µ(D * -D)], (9)
where Z is a normalization constant that ensures that D+ D-P * (D)dD = 1. This result, together with the decoupling form of Eq.( 5) explicitly characterizes the asymptotics of a front of prescribed profile characterized by γ. In practice, it can be showed following classical arguments [START_REF] Van Saarloos | [END_REF] that in the case of a localized initial condition in space the slowest solution is selected. The critical velocity v * (and the corresponding γ * ) is therefore obtained as the minimum of the function v(γ) defined by Eqs.( 7) and ( 8 While in practice v * can be determined by a numerical analysis of Eqs.( 7) and ( 8), a useful and very accurate approximate of v * , which is in fact an exact upper bound, denoted v u * , can be obtained by assuming D * constant in Eq. (7). Under this hypothesis the minimization of v is realized for γ = γ u * ≡ s/D * and yields v u * = 2 √ sD * , where D * is determined by Eq. ( 8) taken at γ = γ u * . These asymptotic results have been checked numerically by solving Eq. ( 1) using a classical Euler numerical scheme, the diffusion part being considered implicitly and the reaction part explicitly, and show that the upper bound v u * provides a very accurate estimate of the exact value v * . Figure 2 shows an example of the obtained population wave which reaches a stationary state of constant velocity after a transient acceleration regime (Fig. 2a). Figure 2b indicates that the population is the most heterogeneous at the edge of the front where the segregation of the fastest phenotypes is observed, while variations of the trait D decrease behind the front. More quantitatively, we find that after the transient accelera-tion regime the wave advances at the predicted velocity v * (see Fig. 3a). Numerical results clearly validate the asymptotic result of Eq. (10) as seen in Fig. 3b,c, which show that both the dependences on u and D of the density n(D, u) are captured by our analysis. Remarkably, the approximate value v u * is numerically extremely close to the exact result v * for all parameters analyzed. Last, Fig. 3d confirms that the segregation of phenotypes with larger diffusion rate D is more important at the leading edge of the front, where it is characterized by the predicted distribution P * (D), while the distribution P u (D) reaches its steady state value P 0 (D) (which is uniform with the choice of Eq.( 2)) at positions x far behind the front.

III. DISCUSSION

Several comments are in order. First, the asymptotic solution of Eq. ( 10) and the determination of the front velocity v * are exact. Similar results, and in particular the exact asymptotic decoupling of the u and D dependences can be obtained along the same lines for other choices of the density current j defined in Eq. ( 2). In particular, the important case where a trait D 0 presents a selective advantage, as modeled by the current defined in Eq. ( 3), can be analyzed in the same way and leads to qualitatively similar results. Second, one finds that for s/α 1, D * → D 0 and P u (D) P 0 (D) = 1/∆D for all u. As expected, when mutations are extremely fast the population is locally well-mixed everywhere with diffusion rate D 0 and the dynamics follows the usual FKPP equation with a front speed v 0 = 2 √ sD 0 . In the more realistic case of slow mutations, s/α is large and the dynamics selects faster individuals at the lead-ing edge so that P * (D) favors larger values of D. In the extreme case of very slow mutations (s/α → ∞) one has D * → D + and P * (D) δ(D -D + ). The front can then be significantly faster than the expected velocity v 0 for a well-mixed population of diffusion rate D 0 since v * → 2 s(D 0 + ∆D/2) > v 0 . Third, it should be emphasized that this segregation of fastest phenotypes, which leads to an acceleration of the wave propagation, is purely dynamical. It is obtained in the low noise limit, and therefore strikingly differs from the mechanism of segregation of neutral mutations studied in [START_REF] Edmonds | [END_REF][3][4]. Last, and importantly, our results qualitatively reproduce the observations (i)-(iii) stated in introduction. Indeed, we found that: (i) fastest phenotypes are segregated at the leading edge of the front, as shown in Figs. 2b,3c. (ii) The proportion of fast phenotypes is higher at the edge and decreases behind the front as shown in Figs. 2b,3d. (iii) The speed of advance of the population wave increases with time, as seen in Fig. 2a, and is larger than the expected propagation speed for a well-mixed population.

To conclude, we have proposed a minimal model of population range expansion in which mutations are neutral and affect only the mobility of individuals. We have shown that such neutral phenotypic variability can yield alone phenotype segregation at the front edge in absence of genetic noise, and significantly impact the dynamical properties of the expansion wave. An exact traveling wave solution can be obtained asymptotically, and shows that phenotype segregation accelerates the front propagation. The results are compatible with available data of invasions of cane toads in Australia or bush crickets in Britain. Applications beyond the field of population dynamics could be expected.

FIG. 1 :

 1 FIG.1: Fisher population waves for (a) a homogeneous population with diffusion rate D0 and for (b) a population with mutations that affect the diffusion rate so that D ∈ [D0 -∆D/2, D0 + ∆D/2]. It is found in case b that a traveling wave solution exists and propagates faster than the classical velocity of a FKPP wave v0 = 2 √ sD0 (case a). The leading edge of the front is mostly populated by fastest individuals (with diffusion rate D * > D0) that are selected dynamically, while the population behind the front is characterized by the steady state value D0.

FIG. 2 :

 2 FIG. 2: Traveling wave solution of Eq.(1) obtained numerically for D0 = 1.5, ∆D = 1, s = 1 and α = 0.01. a Kymograph showing the propagation of the marginal density N (x, t). b Distribution of phenotypes in the traveling wave at a given time t in the stationary state (the front is approximately at x = 90). A larger heterogeneity corresponding to phenotype segregation is found at the leading edge of the front.

FIG. 3 :

 3 FIG. 3: Numerical simulations of Eq.(1) for D0 = 1.5, ∆D = 1 and s = 1. a. Front velocity as a function of time for different values of the mutation rate α. The front is accelerated until it reaches the exact predicted velocity in the stationary state (plain line). b. Shape of the front (log scale) at the leading edge in the stationary state. The theoretical prediction (plain lines) is compared to numerical simulations (symbols and dashed lines) for different values of the mutation rate α. c. Probability distribution of the trait D at the leading edge in the stationary state, where δD ≡ D -D-. The theoretical prediction of Eq. (10) (plain lines) is compared to numerical simulations (symbols and dashed lines). d. Probability distribution of the trait D in the stationary state at different positions x, where δD ≡ D -D-. Here α = 0.01, x f denotes the position of the front and x-= x f -50 and x+ = x f + 10. The theoretical prediction P * (plain line) holds at the leading edge x+ while the distribution is homogeneous at x-, which is the steady state value P0 for the choice of Eq.(2).

Acknowledgments N.M. wishes to thank C. Mouhot for having brought his attention on this subject and for helpful discussions.


Emeric Bouin

Vincent Calvez

Nicolas Meunier

Sepideh Mirrahimi

Benoît Perthame

Gaël Raoul

Raphaël Voituriez

Invasion fronts with variable motility: phenotype selection, spatial sorting and wave acceleration

Invasion fronts in ecology are well studied but very few mathematical results concern the case with variable motility (possibly due to mutations). Based on an apparently simple reaction-diffusion equation, we explain the observed phenomena of front acceleration (when the motility is unbounded) as well as other quantitative results, such as the selection of the most motile individuals (when the motility is bounded). The key argument for the construction and analysis of traveling fronts is the derivation of the dispersion relation linking the speed of the wave and the spatial decay. When the motility is unbounded we show that the position of the front scales as t 3/2 . When the mutation rate is low we show that the canonical equation for the dynamics of the fittest trait should be stated as a PDE in our context. It turns out to be a type of Burgers equation with source term.

Résumé

Fronts d'invasion avec motilité variable : répartition des phénotypes et accélération de l'onde. Les fronts d'invasion en écologie ont été largement étudiés. Cependant peu de résultats mathématiques existent pour le cas d'un coefficient de motilité variable (à cause des mutations). A partir d'un modèle minimal de réactiondiffusion, nous expliquons le phénomène observé d'accélération du front (lorsque la motilité n'est pas bornée), et nous démontrons l'existence d'ondes progressives ainsi que la sélection des individus les plus motiles (lorsque la motilité est bornée). Le point clé pour la construction des fronts est la relation de dispersion qui relie la vitesse de l'onde avec la décroissance en espace. Lorsque la motilité n'est pas bornée nous montrons que la position du front suit une loi d'échelle en t 3/2 . Lorsque le taux de mutation est faible, nous montrons que, dans notre contexte, l'équation canonique pour la dynamique du meilleur trait est une EDP. C'est une équation de type Burgers avec terme source.

Version française abrégée

Dans cette note nous étudions un modèle simple décrivant des fronts invasifs en écologie, pour lesquels la mobilité des individus est sujette à variations. Le modèle, issu de [START_REF] Bénichou | Front acceleration by dynamic selection in Fisher population waves[END_REF], est le suivant,    ∂ t n(t, x, θ) = θ∂ 2 xx n(t, x, θ) + rn(t, x, θ) (1 -ρ(t, x)) + α∂ 2 θθ n(t, x, θ) , x ∈ (-∞, +∞),

θ > 0 , ρ(t, x) = n(t, x, θ) dθ , n(t, -∞, θ) = N (θ) , n(t, +∞, θ) = 0 . (1)
Les conditions aux limites sont complétées ci-dessous. Nous étudions la dynamique de cette équation sous différents régimes. Dans un premier temps nous étudions le problème de propagation de front à vitesse constante. Cela circonscrit l'espace des traits θ ∈ (0, Θ), qui doit être borné Θ < +∞. Nous démontrons que la situation est similaire au cas de l'équation de Fisher-KPP. Il existe une vitesse minimale c * de propagation des ondes de réaction-diffusion (Théorème 1). La relation de dispersion, analogue à celle de Fisher-KPP (un trinôme du second degré en l'occurrence) est donnée via la résolution d'un problème spectral. On lit sur la distribution des phénotypes (le vecteur propre) que les fortes motilités sont favorisées, conclusion opposée au cas des domaines bornés en espace [START_REF] Dockery | The evolution of slow dispersal rates: a reaction diffusion model[END_REF]. Ce même problème spectral intervient lorsqu'on étudie la propagation du front en régime asymptotique hyperbolique (t, x) → (t/ε, x/ε) dans la limite WKB de l'optique géométrique, en suivant l'ansatz

n ε (t, x, θ) = exp(u ε (t, x)/ε)N ε (t, x, θ).
Dans un second temps, nous considérons un espace des traits non borné, Θ = +∞. Dans ce cas le front accélère sans cesse et nous montrons heuristiquement que la loi de propagation du front est naturellement x ∼ α 1/4 r 3/4 t 3/2 . Nous exhibons une solution particulière qui confirme cette heuristique (Proposition 2).

Dans un troisième temps, nous étudions le régime de mutations rares, et nous écrivons une équation canonique pour l'évolution du trait sélectionné localement à l'avant du front en régime asymptotique. La dérivation formelle de cette équation conduit à une équation de transport de type Burgers, avec terme source (Proposition 3). La partie transport est dûe à la progression du front qui "transporte" les individus et donc le trait sélectionné localement. Le terme source est dû à la pression de sélection qui tend à faire augmenter la valeur du trait sélectionné localement.

Introduction

Recently, several works have addressed the issue of front invasion in ecology, where the motility of individuals is subject to variability [START_REF] Champagnat | Invasion and adaptive evolution for individual-based spatially structured populations[END_REF][START_REF] Arnold | Existence of nontrivial steady states for populations structured with respect to space and a continuous trait[END_REF]. It has been postulated that selection of more motile individuals can occur, even if they have no advantage regarding their reproductive rate (spatial sorting) [START_REF] Simmons | Changes in dispersal during species' range expansions[END_REF][START_REF] Kokko | From individual dispersal to species ranges: perspectives for a changing world[END_REF][START_REF] Ronce | How does it feel to be like a rolling stone? Ten questions about dispersal evolution[END_REF][START_REF] Shine | An evolutionary process that assembles phenotypes through space rather than through time[END_REF]. This phenomenon has been described in the invasion of cane toads in northern Australia [START_REF] Phillips | Invasion and the evolution of speed in toads[END_REF]. It has been shown that the speed of the front increases, coincidentally with significant changes in toads morphology. Up to now, only numerical simulations have been proposed to address this issue. Here we analyse a simple model of Bénichou et al [START_REF] Bénichou | Front acceleration by dynamic selection in Fisher population waves[END_REF] which contains the basic features of this process: spatial mobility, logistic reproduction, and variable motility. It is given by equation ( 1) where the last term on the right hand side accounts for modifications of the dispersal rate θ of individuals due to mutations. We consider that mutations are random and they act as a diffusion process in the phenotype space. When needed, we impose Neumann boundary conditions in the variable θ and far-field conditions in the variable x.

Phenotype selection and spatial sorting in the traveling wave

We first consider bounded dispersal rates, say θ ∈ (0, Θ). Following [START_REF] Bénichou | Front acceleration by dynamic selection in Fisher population waves[END_REF], we seek a traveling wave solution to the equation ( 1) connecting 0 to the uniform stationary state N (θ) ≡ Θ -1 . For x large, we make the following ansatz n(t, x, θ) = exp(λ(x -ct))Q(θ), where c > 0 is the speed of the wave, λ < 0 is the spatial decay and Q(θ) denotes the phenotypic distribution of the individuals at the edge of the front. The dispersion relation is equivalent to the following spectral problem: Given a spatial decay rate λ < 0, find c(λ) and an eigenvector Q(θ, λ) such that

     λc(λ) + θλ 2 + r Q(θ, λ) + α∂ 2 θθ Q(θ, λ) = 0 , ∂ θ Q(0, λ) = ∂ θ Q(Θ, λ) = 0 , ∀θ Q(θ, λ) ≥ 0 , Q(θ, λ) dθ = 1 . (2)
The invasion front speed c(λ) is such that 0 is the principal eigenvalue of this spectral problem. Like the Fisher-KPP equation, there is a minimal speed c * > 0 associated with a critical spatial decay λ * < 0.

Theorem 1 (Front propagation and spatial sorting) For all c ≥ c * , there exists formally a traveling front solution n(t, x, θ)

= N (x -ct, θ). It satisfies N (z, θ) ∼ exp(λz)Q(θ, λ) as z → +∞.
The phenotypic distribution at the edge of the front Q(θ, λ) is unbalanced towards more motile individuals.

More precisely, we have

θ edge (λ) := θQ(θ, λ) dθ > Θ 2 , θ edge (λ) -→ α→0 Θ . ( 3
)
This problem is closely related to front propagation in kinetic equations [START_REF] Bouin | A kinetic eikonal equation[END_REF]. There is a natural extension of this result for other mutation operators as integral operators α G(θ, θ )n(t, x, θ ) dθ . In this case, the solution of the spectral problem is deduced from the Krein-Rutman theorem.

Interestingly, we can measure the asymmetry of the phenotypic distribution Q(θ, λ). The relevant quantity here is the mean diffusion coefficient at the edge of the front θ edge (λ). In order to show (3), we integrate the spectral problem (2) over (0, Θ). We get λc + λ 2 θ edge + r = 0. Dividing by Q(θ, λ), and integrating again the spectral problem we get after integration by parts,

(λc + r)Θ + Θ 2 2 λ 2 + α ∂ θ Q(θ, λ) Q(θ, λ) 2 dθ = 0 .
Hence, the mean diffusion coefficient is given by

θ edge (λ) = Θ 2 + α Θλ 2 ∂ θ Q(θ, λ) Q(θ, λ) 2 dθ > Θ 2 . (4)
Individuals at the edge of the front are more motile than at the back of the front. There, the population is homogeneous since N (-∞, θ) ≡ Θ -1 and the averaged diffusion coefficient is θ back = Θ/2. The estimate (4) measures how far the phenotypic distribution differs from the uniform distribution. Finally we find easily that as α vanishes, the distribution Q(θ, λ) concentrates as a Dirac mass on θ = Θ hence the second statement in (3): the most motile individuals are selected. The conclusion is exactly the opposite on bounded spatial domains [START_REF] Dockery | The evolution of slow dispersal rates: a reaction diffusion model[END_REF]. We can derive additional information for the minimal speed c * . Differentiating (2) we obtain

(λc (λ) + c(λ) + 2θλ) Q(θ, λ) + λc(λ) + θλ 2 + r ∂ ∂λ Q(θ, λ) + α∂ 2 θθ ∂ ∂λ Q(θ, λ) = 0 .
Definition of c * ensures that c (λ * ) = 0. We use the notation f := f (θ)Q * (θ) dθ. Since the operator in ( 2) is self-adjoint, we obtain after multiplication by Q(θ, λ) and integration, the relation

c * Q * + 2λ * θQ * = 0 , c * = -2λ * θQ * Q * , (5)
Recalling that λ * c * + (λ * ) 2 θ + r = 0, we can eliminate λ * and we get the following expression for c * ,

(c * ) 2 = 4r θ 1 -1 - θ Q * θQ * 2 -1 > 4r θ .
In other words, the usual formula for the KPP wave speed underestimates the actual minimal speed.

Spatial sorting and the invasion front

Next, we focus on the invasion front. It is natural to perform the hyperbolic rescaling (t, x) → (t/ε, x/ε) in order to catch the motion of the front [START_REF] Freidlin | Geometric optics approach to reaction-diffusion equations[END_REF][START_REF] Evans | A PDE approach to geometric optics for certain semilinear parabolic equations[END_REF]. The new equation writes after rescaling

ε∂ t n ε (t, x, θ) = ε 2 θ∂ 2 xx n ε (t, x, θ) + rn ε (t, x, θ) (1 -ρ ε (t, x)) + α∂ 2 θθ n ε (t, x, θ
) . We perform the partial WKB ansatz in the x variable: n ε (t, x, θ) = exp(u ε (t, x)/ε)N ε (t, x, θ), with the renormalization N ε (t, x, θ) dθ = 1. As ε → 0, the first order expansion leads to solve

∂ t u 0 (t, x)N 0 = θ|∂ x u 0 (t, x)| 2 N 0 (t, x, θ) + r 1 -ρ 0 (t, x) N 0 (t, x, θ) + α∂ 2
θθ N 0 (t, x, θ) . The edge of the front is delimited by the area {u 0 (t, x) < 0}. On this set we have ρ 0 (t, x) = 0 by construction. Therefore we shall solve again the spectral problem (2) for N 0 ≥ 0. Consequently the motion of the front is driven by the eikonal equation built on the effective speed c(λ),

max u 0 , ∂ t u 0 + ∂ x u 0 • c(∂ x u 0 ) = 0 .
The rigorous derivation of this Hamilton-Jacobi equation requires some work. We need basically refined a priori estimates on (u ε ) ε . We formally show the main argument leading to establish the viscosity limit u 0 of u ε in the set {u 0 (t, x) < 0}. We leave the complete proof for future perspectives. Let v 0 be a C 2 test function such that u 0 -v 0 has a strict maximum at (t 0 , x 0 ). The function u ε -v 0 has a maximum at (t ε , x ε ), with (t ε , x ε ) close to (t 0 , x 0 ). Plugging v 0 into the equation satisfied by (u ε , N ε ), namely

∂ t u ε (t, x) -θ|∂ x u ε (t, x)| 2 -εθ∂ 2 xx u ε (t, x) -r + O(ε) N ε (t, x, θ) = α∂ 2 θθ N ε (t, x, θ)
, we obtain at (t ε , x ε ):

∂ t v 0 (t ε , x ε ) -θ|∂ x v 0 (t ε , x ε )| 2 -r + O(ε) N ε (t ε , x ε , θ) ≤ α∂ 2 θθ N ε (t ε , x ε , θ) .
Therefore, N ε is a non-negative, non-trivial subsolution of the spectral problem [START_REF] Dockery | The evolution of slow dispersal rates: a reaction diffusion model[END_REF]. From the characterization of the principal eigenvalue, we have

∂ t v 0 (t ε , x ε ) + ∂ x v 0 (t ε , x ε ) • c(∂ x v 0 (t ε , x ε )) + O(ε) ≤ 0 .
Passing to the limit ε → 0, we obtain that v 0 satisfies

∂ t v 0 + ∂ x v 0 • c(∂ x v 0 ) ≤ 0 at (t 0 , x 0 ). Therefore u 0 is a viscosity sub-solution of the eikonal equation ∂ t u 0 + ∂ x u 0 • c(∂ x u 0 ) = 0 in
the interior of the set {u 0 (t, x) < 0}. The same argument shows that it is also a supersolution and thus a viscosity solution.

Front acceleration

In the case where the set of dispersal rates is unbounded, say θ ∈ (0, +∞), then we cannot solve the spectral problem [START_REF] Dockery | The evolution of slow dispersal rates: a reaction diffusion model[END_REF]. There is no intrinsic speed of propagation, and the front is accelerating as time goes on. Heuristically, we expect the averaged diffusion coefficient to grow linearly with time θ ∼ ( √ αr) t (as for the Fisher-KPP equation set in the phenotype space). Hence the invasion front should scale as x ∼ α 1/4 r 3/4 t 3/2 since the speed is given by c ∼ θ r. Therefore we perform the asymptotic scaling (t, x, θ) → t/ε, x/ε 3/2 , θ/ε , in order to catch the motion of the front asymptotically. The equation writes after rescaling

ε∂ t n ε = ε 2 θ∂ 2 xx n ε + rn ε (1 -ρ ε (t, x)) + ε 2 α∂ 2 θθ n ε . (6)
We perform the WKB ansatz in both variables (x, θ): n ε (t, x, θ) = exp (u ε (t, x, θ)/ε). We derive formally the Hamilton-Jacobi equation for u 0 in the limit ε → 0:

∂ t u 0 -θ|∂ x u 0 | 2 -α|∂ θ u 0 | 2 = r 1 -ρ 0 (t, x) , (7)
where ρ 0 (t, x) is the Lagrange multiplier associated with the constraint max θ u 0 (t, x, θ) ≤ 0. In particular, when the constraint is inactive, i.e. max θ u 0 (t, x, θ) < 0, we have ρ 0 (t, x) = 0. The behavior of the phase function u 0 (t, x, θ) is as follows: the nullset {u 0 (t, x, θ) = 0} propagates in space and phenotype under the action of space mobility, and mutations, combined with growth of individuals. We are able to compute explicitely the evolution of this set for a particular initial data, using Lagrangian Calculus.

Proposition 2 (Front acceleration) For the particular initial data, u 0 (0, x, θ) = 0 if (x, θ) = (0, 0) -∞ if (x, θ) = (0, 0) the location of the nullset {u 0 (t, x, θ) = 0} is given by the following implicit formula (see Fig. 1)

x 2 = 1 9α 2 √ rα t -θ 2θ + 2 √ rα t 2 .
Sketch of proof. The Hamiltonian is given by H((x, θ), (p x , p θ )) = θ|p x | 2 +α|p θ | 2 +r, and the corresponding Lagrangian writes L((x, θ),

(v x , v θ )) = v 2 x /(4θ) + v 2 θ /(4α) -r.
The system of characteristics is given by Ẋ(t) = (2θ(t)p x (t), 2αp θ (t)), and Ṗ (t) = -(0, |p x (t)| 2 ). Using Lagrangian formulation, we deduce after some calculation that, with Z the solution to equation Z 3 + (12θ/α)Z + 24x/α = 0, u 0 (t, x, θ) is given by

u 0 (t, x, θ) = - 1 4αt θ + α 4 Z 2 2 .
This enables to compute the nullset {u 0 (t, x, θ) = 0}. 2

The far-right point of the curve is attained for θ edge = ( √ rα) t. This determines the location of the front. Hence the position of the front in space is exactly, X edge (t) = 4 3 α 1/4 r 3/4 t 3/2 .

Adaptive dynamics at the edge of the front

Interestingly enough, the same equation ( 7) can be derived in the context of adaptive dynamics, a theory which studies mutation-selection processes. It is generally assumed that the mutation process is so slow that mutants can replace the resident species before new mutants arise, if they are better adapted to their environment. This yields a canonical equation which gives the dynamical evolution of the selected trait in the population [START_REF] Dieckmann | The dynamical theory of coevolution: a derivation from stochastic ecological processes[END_REF][START_REF] Champagnat | From individual stochastic processes to macroscopic models in adaptive evolution[END_REF]. Recently, PDE-based methods have been successfully used to derive such a canonical equation for continuous mutation-selection processes [START_REF] Diekmann | The dynamics of adaptation: an illuminating example and a hamilton-jacobi approach[END_REF][START_REF] Barles | Concentrations and constrained Hamilton-Jacobi equations arising in adaptive dynamics[END_REF][START_REF] Lorz | Dirac mass dynamics in multidimensional nonlocal parabolic equations[END_REF]. Here we extend this theory to the case of front invasion coupled with a basic mutation process.

The only difference with ( 1) is that mutations are assumed to be rare (α → ε 2 α):

∂ t n(t, x, θ) = θ∂ 2 xx n(t, x, θ) + rn(t, x, θ) (1 -ρ(t, x)) + ε 2 α∂ 2 θθ n(t, x, θ) .
It is natural to perform a long time rescaling t → t/ε at the scale of evolutionary changes. Then it is useful to rescale space accordingly x → x/ε in order to catch the motion of the front (otherwise it would travel at order O(1/ε)). With these changes of scales, we end up again with equation ( 6), resp. (7) in the WKB limit. We restrict to the edge of the front, namely sup θ u 0 (t, x, θ) < 0, ρ 0 (t, x) = 0. We seek a canonical equation for the locally selected trait θ(t, x) such that u 0 (t, x, θ(t, x)) = sup θ u 0 (t, x, θ).

Proposition 3 (Formal derivation of the canonical equation) The locally selected trait θ(t, x) formally satisfies a Burgers type equation with a source term,

∂ t θ(t, x) -2 θ(t, x)∂ x u 0 ∂ x θ(t, x) = |∂ x u 0 | 2 -∂ 2 θθ u 0 . (8)
The speed of the transport equation is -2θ(t, x)∂ x u 0 . It coincides with the local minimal speed of the traveling front, see e.g. [START_REF] Simmons | Changes in dispersal during species' range expansions[END_REF]. The positive source term accounts for the evolutionary drift which pushes the population towards higher motility (numerical simulations not shown). This equation may yield shock wave singularities, as for the classical Burgers equation, because more motile populations, when located behind less motile populations, will invade them.

Proof. We start from the first order condition ∂ θ u 0 (t, x, θ(t, x)) = 0. We differentiate this relation with respect to t and x, respectively, ∂ 2 tθ u 0 (t, x, θ(t, x)) + ∂ 2 θθ u 0 ∂ t θ(t, x) = 0 , ∂ 2 xθ u 0 (t, x, θ(t, x)) + ∂ 2 θθ u 0 ∂ x θ(t, x) = 0 . On the other hand, we differentiate equation [START_REF] Ronce | How does it feel to be like a rolling stone? Ten questions about dispersal evolution[END_REF] with respect to θ, ∂ 2 θt u 0 -|∂ x u 0 | 2 -2θ∂ x u 0 ∂ 2 θx u 0 -2α∂ θ u 0 ∂ 2 θθ u 0 = 0 . Evaluating the latter at θ = θ(t, x) yields

∂ 2 θt u 0 -2θ(t, x)∂ x u 0 ∂ 2 θx u 0 = |∂ x u 0 | 2 .
Combining these calculations, we conclude that θ satisfies equation [START_REF] Shine | An evolutionary process that assembles phenotypes through space rather than through time[END_REF]. 2

Figure 1 .

 1 Figure 1. (a) Phenotype selection. At the edge of the front the phenotypic distribution is not uniform: more motile individuals are selected. We solved numerically the relation dispersion for (r, α, Θ) = (1, 1, 1). (b) Front acceleration. The set {u 0 = 0} is plotted in the phase space (x, θ), for successive times. The far-right point X edge (t) determines the location of the front.

Email addresses: emeric.bouin@ens-lyon.fr (Emeric Bouin), vincent.calvez@ens-lyon.fr (Vincent Calvez), nicolas.meunier@parisdescartes.fr (Nicolas Meunier), mirrahimi@cmap.polytechnique.fr (Sepideh Mirrahimi), benoit.perthame@upmc.fr (Benoît Perthame), g.raoul@damtp.cam.ac.uk (Gaël Raoul), voiturie@lptmc.jussieu.fr (Raphaël Voituriez).

Acknowledgements. S. M. benefits from a 2 year "Fondation Mathématique Jacques Hadamard" (FMJH) postdoc scholarship. She would like to thank Ecole Polytechnique for its hospitality.


Johan Oudinet
email: oudinet@kit.edu

Alain Denise
email: denise@lri.fr

Marie-Claude Gaudel

A new dichotomic algorithm for the uniform random generation of words in regular languages

Keywords: random generation, regular languages, automata

We present a new algorithm for generating uniformly at random words of any regular language L. When using floating point arithmetics, its bit-complexity is O(q log 2 n) in space and O(qn log 2 n) in time, where n stands for the length of the word, and q stands for the number of states of a finite deterministic automaton of L. We implemented the algorithm and compared its behavior to the state-of-the-art algorithms, on a set of large automata from the VLTS benchmark suite. Both theoretical and experimental results show that our algorithm offers an excellent compromise in terms of space and time requirements, compared to the known best alternatives. In particular, it is the only method that can generate long paths in large automata.

Introduction

The problem of randomly and uniformly generating words from a regular language was first addressed by [START_REF] Hickey | Uniform random generation of strings in a contextfree language[END_REF], as a particular case of context-free languages. Using the so-called recursive method [START_REF] Wilf | A unified setting for sequencing, ranking, and selection algorithms for combinatorial objects[END_REF][START_REF] Flajolet | A calculus for the random generation of labelled combinatorial structures[END_REF], they gave an algorithm in O(qn) space and time for the preprocessing stage and O(n) for the generation, where n denotes the length of the word to be generated, and q denotes the number of states of a deterministic finite automaton of L. Later, [START_REF] Goldwurm | Random generation of words in an algebraic language in linear binary space[END_REF] showed that the memory space can be reduced to O(q), by using a parsimonious approach to keep in memory only a few coefficients. These results are in terms of arithmetic complexity, where any number is supposed to take O(1) space and each basic arithmetic operation takes O(1) time. As for the bit complexity, the above formulas must be multiplied by O(n) due to the exponential growing of the involved coefficients according to n.

If floating point arithmetic is used [START_REF] Denise | Uniform random generation of decomposable structures using floating-point arithmetic[END_REF] for the [START_REF] Hickey | Uniform random generation of strings in a contextfree language[END_REF] algorithm, almost uniform generation can be performed in O(qn log n) bit complexity for the preprocessing stage (in time and memory), and O(n log n) for the generation. Meanwhile, the parsimonious version of Goldwurm cannot be subject to floating point arithmetic, because of the numerical instability of the involved operations [START_REF] Oudinet | Random exploration of models[END_REF].

Another technique, the so-called Boltzmann generation method [START_REF] Duchon | Boltzmann samplers for the random generation of combinatorial structures[END_REF], is well fitted for approximate size generation: it makes possible to generate words of size between (1 -ε)n and (1 + ε)n, for a fixed value ε, in average linear time according to n [START_REF] Flajolet | Boltzmann sampling of unlabelled structures[END_REF]. As for exact size generation, the average complexity of generation is in O(n 2 ), although it can be lowered to O(n) if the automaton of L is strongly connected. Bolztmann generation needs a preprocessing stage whose complexity is in O(q k log k n), for some constants k ≥ 1 and k ≥ 1 whose precise values are not given in [START_REF] Duchon | Boltzmann samplers for the random generation of combinatorial structures[END_REF] and subsequent papers, to our knowledge.

Recently, [START_REF] Bernardi | A linear algorithm for the random sampling from regular languages[END_REF] developed a new divide and conquer approach for generating words of regular languages, based on the recursive method. If using floating point arithmetic, their algorithm runs in O(qn log(qn)) in the worst case, with a preprocessing in O(q 3 log n log 2 (qn)) time and O(q 2 log n log(qn)) space. Moreover the average complexity of the generation stage can be lowered to O(qn) if using a bit-by-bit random number generator.

Here we present a new algorithm named dichopile, also based on a divideand-conquer approach, although drastically different from the above one. When using floating point arithmetics, its bit-complexity is O(q log 2 n) in space and O(qn log 2 n) in time.

The paper is organized as follows. In Section 2 we present the dichopile algorithm and we compute its complexity in space and time. Section 3 is devoted to experiments: we compare the running time and space requirements of C++ programs that implement dichopile and the other state-of-the-art algorithms, on a set of large automata from the VLTS benchmark suite [START_REF] Garavel | Very large transition systems[END_REF]). Finally we discuss in Section 4 the advantages and drawbacks of each of the algorithms, from both theoretical and experimental points of view. As shown in Table 4, it turns out that, compared to the known best alternatives, our algorithm offers an excellent compromise in terms of space and time requirements.

The Dichopile algorithm

At first, let us briefly recall the general principle of the classical recursive method for regular languages. Let us consider a deterministic finite automaton of L with q states {1, 2, . . . , q}. Obviously, there is a one-to-one correspondence between the words of L and the paths in A starting at the initial state and ending at any final state. For each state s, we write l s (n) for the number of paths of length n starting from s and ending at a terminal state. Such values can be computed with the following recurrences on n (where F denotes the set of final states in A):

       l s (0) = 1 if s ∈ F l s (0) = 0 if s ∈ F l s (i) = s→s l s (i -1) ∀i > 0 (1) Let L n = l 1 (n), l 2 (n), . . . , l q (n) .
The recursive method proceeds in two steps:

• Compute and store L k for all 1 ≤ k ≤ n. This calculation is done starting from L 0 and using Formula 1.

• Generate a path of length n by choosing each state according to a suitable probability to ensure uniformity among every path of length n. Thus, the probability of choosing the successor s i when the current state is s and the path has already n -m states is:

P(s i ) = l si (m -1) l s (m) . (2)
Note that in order to choose the successor of the initial state, we only need L n and L n-1 . Then, L n-1 and L n-2 allow to choose the next state and so on. Thus, if we have a method that compute efficiently L n , L n-1 , . . . , L 0 in descending order, we can store the two last vectors only and reduce space complexity compared to recursive method, which stores all L k 's in memory. This inverse approach constitutes the principle of Goldwurm (1995)'s method. In [START_REF] Oudinet | Random exploration of models[END_REF], the inverse recurrence is stated for regular languages, and it is shown that the algorithm is numerically instable, thus forbidding the use of floating-point arithmetics.

General principle

The idea of our dichopile algorithm is as follows. Compute the number of paths of length n from the number of paths of length 0 while saving in a stack a logarithmic number of intermediate steps: the number of paths of length n/2, of length 3n/4, of length 7n/8, etc. When we need to compute the number of paths of length n -i, we compute it again from the intermediate stage that is at the top of the stack. Figure 1 illustrates the principle of this algorithm. Recall that L j denotes the vector of q numbers of paths of length j, that is the l s (j)'s for all states s.

Algorithm 1 draws uniformly at random a path of length n by successively computing the numbers of paths of length i in descending order. This algorithm takes as inputs:

• a deterministic finite automaton A;

n/2 3n 4 n [ ] L 0 L n/2 L 3n/4 L n . . . Step 0: compute L n from L 0 j n -i n [ ] L 0 L j . . .
Step i: pop L j from the stack j

j+n-i 2 n -i n [ ] L 0 L j . . .

L k

Step i: compute L n-i from L j Algorithm 1 Draw a path of length n with the dichopile algorithm.

Require: an automaton A, a vector L 0 , a length n, and a function F that computes L j from L j-1 Ensure: returns a path σ of length n.

s ← s 0 {initialize s to the initial state} push (0, L 0 ) for i ← 0 to n do {Iteration for computing L n-i } (j, L cur) ← top of the stack if j > n -i then {useless value on the stack, get the next one} pop from the stack (j, L cur) ← top of the stack • a vector L 0 of q numbers defined as in the first two items of Formula 1;

end if while j < n -i -1 do {compute L n-i from L j } k ← j+n-i 2 for m ← j + 1 to k do L cur ← F (L cur) end for push (k, L cur) {push L k to the stack} j ← k end while if j = n -i -1 then L cur ← F (L cur) end if if i > 0 then {wait until L suc is defined to have L cur = L n-i and L suc = L n-i+1 } choose the next transition t i in A according to s, L cur and L suc σ ← σ.t i {concatenation of the transition} s ← the extremity of t i end if L suc ← L cur end for return σ L k L n lg(n -k) After step 0 L k L n-i lg(n -k -i) After step i
• the path length n;

• and a function F that computes L j from L j-1 , thus F computes the l s (j)'s for all s ∈ S using the third item of Formula 1.

It uses a stack and two local variables L cur and L suc, which are vectors of q numbers. The vector L suc saves the previous value of L cur.

Step i = 0 computes L n from L 0 , pushing to the stack the vectors L n/2 , L 3n/4 , L 7n/8 , etc. All the divisions are integer divisions.

Step i seeks to compute L n-i . For that, it starts by retrieving the top of the stack L j (if j > n -i then it takes the next item on the stack) and computes L n-i from L j , pushing to the stack a logarithmic number of intermediate vectors L k . At the end of each iteration of the main loop and just before updating L suc, we have L cur = L n-i and L suc = L n-i+1 . Then, using Formula 2, we can draw the successor state according to the current state and values contained in these two vectors.

Complexity analysis

Theorem 1. Using floating-point numbers with a mantissa of size O(log n), bit complexities of Algorithm 1 are O(q log 2 n) in space and O(dqn log 2 n) in time, where d stands for the maximal out-degree of the automaton.

Proof. Unlike the classical recursive method, there is no preprocessing phase. Values cannot be saved between two path generations because the content of the stack changes during the drawing.

The space complexity depends on the stack size. After the first iteration, there are lg n elements on the stack, and there will never be more elements on the stack in subsequent iterations: if L k is on top of the stack at the i-th iteration, there were lg(n -k) elements above L k after the first iteration and there will be lg(n -k -i) elements after the i-th iteration; Hence, fewer elements. This property is illustrated in Figure 2. Each stack element contains an integer (≤ n) and q path numbers represented by floating-point numbers with a mantissa of size O(log n), i.e. O(q log n) bits per item. Hence the size occupied by the stack is O(q log 2 n) bits.

The time complexity depends on the number of calls to function F and this number depends on the difference between the last stacked value (j) and the one to compute (n-i). To calculate this complexity, we refer to a diagram describing the successive calculations done by the algorithm using two functions, f and g, which call each other. The calculation scheme is shown in Figure 3. The function f counts the number of calls to function F. For f (n), it will first call F n times. Then, since the first saved value is n/2 , it will call f ( n/2 -1). The next saved value is 3n/4 , hence the call to f ( n/4 -1), and so on.

0 n/2 3n 4 7n 8 n f (n) f (n/2) f (n/4) f (n/8) f (1)
Thus, the number of calls to F to generate a path of length n is equal to f (n), which is defined as:

f (1) = 1 f (n) = n + lg n i=1 f ( n 2 i -1)
As the number of terms in the sum is not a constant, we cannot apply directly the Master theorem nor the Akra-Bazzi theorem [START_REF] Akra | On the solution of linear recurrence equations[END_REF]. However, their theorem 1 does apply to the generalized case where the number of terms in the sum, k, is a variable of n. This theorem allows us to ignore both the floor function and the minus 1 in the index of f since the same asymptotic behavior is kept. Thus, we study the asymptotic behavior of T defined as follows:

T (n) = n + lg n i=1 T ( n 2 i ) Lemma 2.
The function T has the following recurrence form:

T (n) = 2T (n/2) + n/2 ∀n ≥ 2 Proof. T (2) = 2 + T (1) = 3 = 2T (1) + 1 T (2n) = 2n + lg 2n i=1 T ( 2n 2 i ) = 2n + 1+lg n i=1 T ( n 2 i-1 ) = 2n + lg n i=0 T ( n 2 i ) = n +   n + lg n i=1 T ( n 2 i )   + T (n) = 2T (n) + n
Thanks to Lemma 2, we can apply the master theorem (Cormen et al., 2001, sect. 4.3) on T . As n/2 ∈ Θ(n), we are in the second case of the master theorem, which states:

T (n) ∈ Θ(n log n)
Hence a time complexity of f (n) in Θ(n log n) calls to function F . The cost of a call to the function F is in O(dq log n) because it corresponds to compute l s (i) from L i-1 for all s ∈ S, using a floating-point arithmetic with numbers of O(log n) bits. Thus, a bit complexity of O(dqn log 2 n) in time.

Experimental results

In this section, we present the experiments we did to measure and compare the relative efficiency of Dichopile with regards to the other algorithms for the uniform random generation of words in regular languages. We first describe the implementation and the methodology used to conduct those experiments in Section 3.1. Then, we present an explanation together with empirical evidence of the numerical instability of Goldwurm's method in Section 3.2. As written previously the proof is given in details in [START_REF] Oudinet | Random exploration of models[END_REF], but we also briefly present it here for the sake of completeness. Finally, we draw cross views of the time and the memory used for each algorithm (except Boltzmann) to generate paths in various automata, from 289 states to more than 10 7 states.

Implementation and methodology

All the algorithms for the uniform random generation of words in regular languages, but Boltzmann, are implemented and freely available in our free C++ library: Rukia1 . This library is based on several other libraries: the BGL (Boost Graph Library) to manipulate graphs [START_REF] Siek | Boost graph library[END_REF], the GMP (Gnu Multiple Precision) library [START_REF] Granlund | GNU MP: The GNU Multiple Precision Arithmetic Library[END_REF] and the Boost Random library [START_REF] Maurer | Boost random number library[END_REF] in order to generate random numbers. Note that Goldwurm's method should not use a floating-point arithmetic because of its numerical instability (see Section 3.2 for an in-depth explanation).

Unfortunately we could not implement the Boltzmann method: it needs to precompute a numerical parameter whose value depends on the generating series of the language. And computing this value for huge automata is still an unsolved problem, up to our knowledge.

In order to evaluate those algorithms, we have selected a pool of automata from the VLTS (Very Large Transition Systems [START_REF] Garavel | Very large transition systems[END_REF])) benchmark suite. These automata correspond to real industrial systems. Each name of an automaton writes vasy X Y, where X is the number of states divided by 1000, and Y is the number of transitions divided by 1000. The general properties of these automata can be found in Table 1. The number of paths, in those automata, grows exponentially according to the path length. Thus, there are for example more than 10 479 paths of length 1000 in vasy 1 4 automaton.

We did all our experiments on a virtual machine on a dedicated server whose hardware is composed of an Intel Xeon 3GHz processor with 32GiB memory. An experiment was to measure the amount of memory used and the elapsed time by an algorithm to generate a path of a certain length in an automaton. Each experiment was run 5 times and only the average value is presented (The standard deviation is always lower than 5% of the average value, which means it is not necessary to repeat the experiment more than 5 times). In case where the path generator was too fast to be measured by the timer (less than 10 milliseconds), 100 paths were generated and the average time was divided by 100.

Numerical instability of Goldwurm's method

Denise and Zimmermann (1999) proved that the operations used by the recursive method are numerically stable. We experimentally checked that it is also true for our implementation of the recursive method: Figure 4 shows that the maximal relative error is less than 10 -17 ; hence, our implementation is numerically stable. The relative error is computed as follows:

Err = max s∈S 0<i≤n |l s (i) -l s (i)| l s (i) (3)
where l s (i) (resp. l s (i)) is the number of paths of length i starting from the state s and computed with a floating-point (resp. exact) arithmetic. On the other hand, the extra operations performed by Goldwurm's method are numerically unstable as illustrated by the maximal relative error computed on our implementation of this method (see Figure 5). Such extra operations for rational languages correspond to the computation of the following formulae (see [START_REF] Oudinet | Random exploration of models[END_REF] for details):

l s (i) = s r s,s × l s (i + 1) with r s,s ∈ Q If ε s
i is the error of l s (i), then the error propagation is as follows:

ε s n-1 = t ε t n × |r st | ε n = max s ε s n ε n-1 ≤ max s t ε n |r st | ε n-1 ≤ max s t |r st | ε n
Thus, the error propagation follows a geometric distribution with parameter max s,s |r s,s |. In other word, Goldwurm's method is numerically unstable as soon as there is a |r s,s | greater than 1, which is the case in the automata studied.

Memory requirement

Here we present the memory consumption for the three methods that use floating-point arithmetic: the recursive method [START_REF] Denise | Uniform random generation of decomposable structures using floating-point arithmetic[END_REF], the Divide and Conquer method [START_REF] Bernardi | A linear algorithm for the random sampling from regular languages[END_REF] the recursive method (Rec) using floating-point arithmetic, the divide and conquer method (DaC), and the Dichopile method. The symbol means that the method needs more than 32GiB of memory, and the symbol -means the value is unavailable due to a timeout during the preprocessing

Table 2 shows the results of comparing the memory requirement for the three methods, when generating a path of a defined length in various automata. Even though the space complexity of Bernardi and Giménez's divideand-conquer (DaC) method is close to the space complexity of the recursive method in terms of path length, the difference of space consumption for those two methods is huge due to the size of automata considered here. In consequence, the divide-and-conquer method could not be used with automata of 5-thousand states and more. Dichopile performed much better than the two other methods and was the only one capable of generating paths in an automaton of more than 10 7 states.

Execution time

In addition to the measurement of space consumption, we also measured elapsed time for the three same methods in order to generate a path in various automata. Results for the other methods are provided in Appendix A.

Table 3 shows elapsed time for the three methods. The preprocessing time (Pre), which has to be done only once whatever the number of paths generated, is separated from the running time (Draw)2 . For example, the elapsed time necessary to generate 1000 paths with one method is equal to Pre + 1000 × Draw.

Figure 6 shows memory consumption and elapsed time to generate 1000 paths for each method. In those experiments, the recursive method was faster than the two other methods but Dichopile was capable of generating longer paths in larger automata. Again, even though the time complexity of divideand-conquer is the best with regards to the path length, the automata considered in our experiments are too large to see the advantage of divide-and-conquer.

Finally, the recursive method is fast but does not scale up very well (both in terms of long path or large automaton). The divide-and-conquer algorithm is well-suited for very long paths but for small automata only. Thus, the Dichopile algorithm is an efficient algorithm to generate long paths in large automata. Table 3: Comparison of the execution time to generate paths in various automata by: the recursive method (Rec) using floating-point arithmetic, the divide and conquer method (DaC), and the Dichopile method. The symbol ∞ means that the generation was not complete after 24 hours of execution, and the symbol -means the value is unavailable due to a memory overflow during the preprocessing. 2 andTable 3, assuming 1000 paths are generated. Thus, each time value corresponds to Pre + 1000 × Draw.

Conclusion

Table 4 summarizes the bit complexities in time and in space of the known algorithms to generate random words in regular languages. Since our goal is to be able to explore at random very large models, we are interested in the complexity in terms of both the path length n and the automaton size q. For the sake of clarity, we consider as constants the following values: the maximum degree d of the automaton and the mantissa size b chosen for the floating-point numbers. The inverse method using floating-point arithmetic is crossed out since it can not be used due to its numerical instability.

Table 4: Summary of the binary complexities in time and in space according to the method used. We consider the path length n and the number of states q in the automaton.

Time

Method

Arith Space Preprocessing Generation

recursive exact O(qn 2 ) O(qn 2 ) O(n 2 ) inverse exact O(qn) O(q 2 + qn 2 ) O(qn 2 ) recursive float O(qn log n) O(qn log n) O(n log n) inverse float O(q log n) O(q 2 + qn log n) O(qn log n) divide and conquer float O(q 2 log n log(qn)) O(q 3 log n log 2 (qn)) O(qn log(qn)) Boltzmann float O(q) O(q k log k n) O(n 2 ) dichopile float O(q log 2 n) O(1) O(qn log 2 n)
From the point of view of space complexity, obviously the best algorithm is the Boltzmann method with its O(q) complexity. The main limitation of the recursive method is the space needed to store the counting table. Even when using floating-point arithmetic, the space complexity is still O(qn log n), which becomes very problematic for large n and q. The inverse method has similar problems, with its O(qn) complexity. Both divide and conquer and dichopile perform well due to their polylogarithmic complexity in n, but dichopile uses more than q times less memory than divide and conquer (up to a constant factor). This is illustrated by the experimental results (Table 2) where dichopile can manage huge automata while divide and conquer cannot. On the other hand, for small automata and long paths, divide and conquer uses much less memory than the recursive method, as expected.

If for any reason, the tiny difference from the uniformity induced by the use of a floating-point arithmetic is not acceptable, the inverse method can be used as it offers the best space complexity, but it requires a long generation time (experimental results for this method can be seen in Appendix A.2). Exactly uniform generation can also be done by combining floating-point arithmetic and exact arithmetic [START_REF] Denise | Uniform random generation of decomposable structures using floating-point arithmetic[END_REF], but the space complexity becomes larger than for the quasi-uniform generation shown here. For example, the bit space complexity becomes O(qn 2 ) for the recursive method and O(q 2 n) for the divide-and-conquer algorithm.

Regarding the time complexity only, the best algorithm is the classical recursive scheme with its O(n log n) complexity in floating point arithmetic. As for Boltzmann, the rejection procedure necessary to obtain paths of length n raises the time complexity to O(n 2 ) in the general case. Both divide and conquer and dichopile are linear in q and almost linear in n, with an advantage for divide and conquer. Regarding the experiments (Table 3), this hierarchy is respected: when divide and conquer has enough memory for running, it runs faster than dichopile; and both are outperformed by the recursive method. Note that, as stated by [START_REF] Bernardi | A linear algorithm for the random sampling from regular languages[END_REF], faster versions of these three algorithms could be implemented, by using a bit-by-bit random number generator.

Altogether, the classical recursive method is fast (after the preprocessing stage), but is unusable for large n and q due to its huge space requirement. On the other hand, the Boltzmann method needs little memory but is slow3 according to n. Unfortunately, as written in Section 3, we could not experiment the Boltzmann method because computing the value of its parameter for large automata is still an unsolved problem, up to our knowledge. Both divide and conquer and dichopile are excellent compromises. From a theoretical point of view, the latter offers a better space complexity, and comparatively the increase of time complexity according to n is quite acceptable. Experiments confirm that divide and conquer is well fitted for small automata and long words: it needs much less memory than the recursive method and the time for generating a word remains very reasonable. On the other hand, dichopile is the only method that can generate long paths in large automata.

Figure 1 :

 1 Figure 1: Principle of the dichopile algorithm.

Figure 2 :

 2 Figure 2: The stack size is maximal after the first iteration of the dichopile algorithm.

Figure 3 :

 3 Figure 3: Recursive scheme of the number of calls to function F done by Algorithm 1. We omit floor and ceiling notations to clarify the figure.

Figure 4 :

 4 Figure 4: Measurement on two automata of the maximum relative error for the recursive method, defined by Equation (3).

Figure 5 :

 5 Figure 5: Measurement on two automata for the Goldwurm's method of the maximum relative error, defined by Equation (3).

Figure 6 :

 6 Figure 6: Cross view of the results presented in Table2 and Table 3, assuming 1000 paths are generated. Thus, each time value corresponds to Pre + 1000 × Draw.

Table 1 :

 1 Detailed description of the VLTS benchmark suite.

	Name	# states # transitions # labels Branching factor avg [min -max]
	vasy 0 1	289	1224	2	4.24 [4 -8]
	vasy 1 4	1183	4464	6	3.77 [2 -5]
	vasy 5 9	5486	9676	31	1.76 [0 -6]
	vasy 10 56	10849	56156	12	5.18 [4 -6]
	vasy 12323 27667 12323703	27667803	119	2.25 [0 -13]

Table 2 :

 2 Comparison of the memory requirement to generate paths in various automata by:

http://rukia.lri.fr

The stack used by Dichopile has to be computed again from scratch in order to generate a new path. Thus, there is no preprocessing time for this method.

Meanwhile, if an error margin is tolerated on the path length, that is if it suffices to generate paths whose length lay in the interval [(1 -ε)n, (1 + ε)n] for a fixed ε > 0, then the time complexity of the Boltzmann method is O(n).

Acknowledgement

We warmly thank Carine Pivoteau for helpful discussions on Boltzmann generation.

Appendix A. Experimental results using exact arithmetic

Beside results presented in Section 3, we also conducted experiments using exact arithmetic for the recursive method and the Goldwurm's method. Comparing performances of those methods with the ones using floating-point arithmetic is not useful because methods with floating-point arithmetic are faster and need less memory than their counterparts using exact arithmetic. For the sake of completeness, we present here results for the recursive method then for the Goldwurm's method using exact arithmetic.

Appendix A.1. The recursive method using exact arithmetic

In this section, we measured the performance of the original algorithm from [START_REF] Hickey | Uniform random generation of strings in a contextfree language[END_REF]. Its bit complexity is in O(qn 2 ) space and time for the preprocessing stage and O(n 2 ) for the generation, where n denotes the length of the path to be generated, and q denotes the number of states of a deterministic finite automaton.

Table A.5 shows the memory usage and the elapsed time to generate a path depending on the length of the path and the size of the automaton. We note that the generation is very fast but the space needed to store the counting table is very important and thus limits the length of paths that can be drawn. Appendix A.2. The inverse recursive method using exact arithmetic

In this section, the parsimonious version of Goldwurm is used to keep in memory a few coefficients only. As this version is numerically unstable, floatingpoint arithmetic cannot be used. The space complexity is in O(qn) bits. Generation of a path of length n is done in O(qn 2 ), after a pre-processing time in O(q 2 + qn 2 ).

Table A.6 shows memory usage and elapsed time to generate a path in the automata vasy 0 1, vasy 1 4 and vasy 5 9. Path length varies from 1000 to 128, 000. The computation of the values needed for reversing the recurrences was not possible for the other automata, due to insufficient memory resources.

Thanks to reducing space complexity, it is possible to generate much longer paths. Thus, a path of length 64, 000 was generated with this method in the automaton vasy 5 9. While with the recursive method using exact arithmetic, we could not get paths of length 16, 000 in the same automaton. However, the generation time is much longer: almost 10 minutes for a path of length 8, 000 while the same path can be generated in 24 milliseconds by the recursive method. The symbol ∞ means that the generation was not complete after 24 hours of execution.