Freddy Bouchet

Tobias Grafke

Tomás Tangarife

Eric Vanden-Eijnden

Large Deviations in Fast-Slow Systems

Keywords:

The incidence of rare events in fast-slow systems is investigated via analysis of the large deviation principle (LDP) that characterizes the likelihood and pathway of large fluctuations of the slow variables away from their mean behavior -such fluctuations are rare on short time-scales but become ubiquitous eventually.

This LDP involves an Hamilton-Jacobi equation whose Hamiltonian is related to the leading eigenvalue of the generator of the fast process, and is typically non-quadratic in the momenta -in other words, the LDP for the slow variables in fast-slow systems is different in general from that of any stochastic differential equation (SDE) one would write for the slow variables alone. It is shown here that the eigenvalue problem for the Hamiltonian can be reduced to a simpler algebraic equation for this Hamiltonian for a specific class of systems in which the fast variables satisfy a linear equation whose coefficients depend nonlinearly on the slow variables, and the fast variables enter quadratically the equation for the slow variables. These results are illustrated via examples, inspired by kinetic theories of turbulent flows and plasma, in which the quasipotential characterizing the long time behavior of the system is calculated and shown again to be different from that of an SDE.

I. INTRODUCTION

The evolution of many dynamical systems of interest involve the interplay of fast and slow variables. Examples include planetary motion, geophysical flows, climate-weather interaction models, macromolecules, etc. In such systems one is typically interested in the behavior of the slow variables on time-scales that are much longer than that over which the fast variables evolve.

Under suitable conditions, the fast variables are adiabatically slaved to the slow ones, and the latter only feel the average effects of the former. When this is the case the evolution of the slow variables on their natural time-scale can be captured by a closed limiting equation for these variables alone that is obtained by averaging out the effect of the fast variables on the slow motions. This equation is valid in the limit when the scale separation between the fast and slow variables is infinitely wide, and it is an instance of the Law of Large Numbers (LLN) in the present context. Of course the scale separation is never infinite in reality and the slow variables also experience fluctuations above their average motion. Small fluctuations are captured by the Central Limit Theorem (CLT) which provides a linear stochastic differential equation (SDE) for the difference between the slow variables and their mean. Large deviations from the LLN away from the CLT scaling, on the other hand, can be characterized via a large deviation principle (LDP). On time scales that are of order one with respect to the clock of the slow variables, these large fluctuations are rare events. However, the LDP also captures the long time behavior of these variables when the effect of fluctuations is no longer negligible and large deviations from the LLN are no longer rare. This is the case, for example, if the limiting equation given by the LLN possesses multiple stable fixed points (or, more generally, multiple attractors). In such situations, fluctuations may eventually push the system from the vicinity of one such attractor to another, and the way this occurs is typically not captured by the CLT, but rather by the LDP. The aim of the present paper is to analyze in detail the structure of the LDP for a specific class of fast-slow systems and thereby provide concrete tools to characterize rare events as well as noise induced transitions in such systems. We note that this question has been investigated by many authors (see e.g. [START_REF] Iosifovich | The averaging principle and theorems on large deviations[END_REF][START_REF] Kifer | Averaging in dynamical systems and large deviations[END_REF][START_REF] Kifer | Averaging principle for fully coupled dynamical systems and large deviations[END_REF][START_REF] Liptser | Large deviations for two scaled diffusions[END_REF][START_REF] Veretennikov | On large deviations for sdes with small diffusion and averaging[END_REF]), but few concrete results exist that permit to actually compute the Hamiltonian associated with the LDP: our goal here is to perform such calculations explicitly in simple examples or indicate how they could be performed numerically in more complicated ones. The examples we deal with have the peculiarity that, on the one hand the fast variables evolve linearly once the slow variable is held fixed, and on the other hand the fast variables act on the slow variable through a quadratic nonlinearity. Such a structure is very common in many examples related to the kinetic theories of both turbulent flows and plasma physics. It appears naturally for instance for the kinetic theory of the 2D Navier-Stokes or the family of quasigeostrophic models (see e.g. [START_REF] Bouchet | Kinetic theory of jet dynamics in the stochastic barotropic and 2d navier-stokes equations[END_REF] and references therein) or the kinetic theory of plasma physics leading to either the Lenard-Balescu or the Vlasov equation (see e.g. [START_REF] Nicholson | Introduction to plasma theory[END_REF]), or more generally the kinetic theory of systems with long range interactions (see e.g. [START_REF] Campa | Physics of long-range interacting systems[END_REF]).

A. Set-up

We will consider fast-slow systems of the type

       Ẋ = f (X, Y ) dY = 1 α b(X, Y )dt + 1 √ α σ(X, Y )dW (t) (1)
Here f : R m×n → R m and b : R m×n → R n are vector fields, W (t) ∈ R p is a standard pdimensional Wiener process, σ : R m×n → R n×p , and α > 0 is a parameter whose smallness measures the separation of time scale between the slow X ∈ R m and the fast Y ∈ R n . To analyze the behavior of the slow variables X when α 1, let us introduce the virtual fast process

d Ỹx = b(x, Ỹx )dτ + σ(x, Ỹx )dW (τ ) (2)
where x is fixed. This equation is obtained from the equation for Y in (1) by setting X to the fixed value x and rescaling time to the natural time scale of the fast Y , τ = t/α. Assume that the virtual fast process is ergodic at every x with respect to the invariant measure µ x (dy) (which may depend parametrically on x) and that the following expectation exists

F (x) = R n f (x, y)µ x (dy) ≡ (P f )(x) = lim T →∞ 1 T T 0 f (x, Ỹ x (τ ))dτ (3)
where the second equality follows from ergodicity. Then for any ε > 0 and any fixed T < ∞ we have

lim α→0 P x sup 0≤t≤T |X(t) -X(t)| < ε = 1 (4)
where P x denotes the probability conditional on X(0) = X(0) = x and X satisfies the limiting equation

Ẋ = F ( X) (5)
Property ( 4) is an instance of the LLN in the present context. For the reader's convenience, and also because the derivation of (5) involve formal asymptotic tools that we will need below, we recall this derivation in Appendix A.

The limiting equation in ( 5) is deterministic because we assumed that there is no explicit noise acting on the slow variables X in (1) and we let α → 0. In this limit the effect of the fast Y on the slow X completely averages out. For small but finite α, however, the slow variables are subject to fluctuations above their mean. To leading order, these fluctuations can be captured by the Central Limit Theorem (CLT). More precisely, denote by C f (x, τ ) the time-correlation matrix of f (x, y) ≡ f (x, y) -F (x) along the virtual fast process Ỹx ,

C f (x, τ ) = R n E y f (x, Ỹx (τ )) f T (x, y) + f (x, y) f T (x, Ỹx (τ )) µ x (dy) (6)
where E y denotes expectation over the virtual fast process conditional on Y x (0) = y. Assume that the following integral exists

A(x) = ∞ 0 C f (x, τ )dτ, = lim T →∞ 1 T T 0 dτ T 0 dτ f (x, Ỹx (τ )) f T (x, Ỹx (τ )) (7)
and the following expectation exists

B(x) = R n ∂ x f (x, y)µ x (dy) + ∞ 0 dτ R n ∂ y E y f (x, Ỹ x (τ )) ∂ x b(x, y)µ x (dy) (8)
where (∂ y f ) i,j = ∂f i /∂y j and (∂ x b) i,j = ∂b i /∂x j . Then on any interval t ∈ [0, T ] with T < ∞,

the process ξ = X - X √ α (9)
converges in distribution towards the Gaussian process solution of

dξ = B( X)ξdt + η( X)dW (t) (10)
where η(x) is a m × m matrix such that (ηη T )(x) = A(x) and W (t) is a m-dimensional Wiener process -the derivation of (10) via formal asymptotic expansion techniques is recalled in Appendix B.

While the CLT indicates that typical fluctuations of the slow variables around their mean are of order O( √ α) on time-scales that are O(1) in α, it does not permit to estimate the probability of deviations of order one away from this mean. On these time-scales, large deviations are expected to be rare, and their probability can be estimated by a LDP which takes the following form. Suppose that the following expectation over the virtual fast process exists

H(x, θ) = lim T →∞ 1 T log E y exp θ • T 0 f (x, Ỹx (τ ))dτ (11)
and define the Lagrangian associated with this Hamiltonian

L(x, y) = sup θ (y • x -H(x, θ)) (12)
as well as the action

I T (x) = T 0 L(x(t), ẋ(t))dt (13)
Then this action permits to estimate the probability that the slow process wanders away from its mean behavior in the sense that for any

Γ ⊂ {γ ∈ C([0, T ], R m ) : γ(0) = x} the following LDP holds -inf γ∈Γ • I T (γ) ≤ lim inf α→0 α log P(X ∈ Γ) ≤ lim sup α→0 α log P(X ∈ Γ) ≤ -inf γ∈ Γ I T (γ) (14)
The action (13) is also useful if one is interested in the slow motions on longer time-scales that can be O(α -1 ) or even O(exp(C/α)) for C > 0. On these time-scales, large deviations stop being rare and may, for example, lead to random transitions between the different attractors of the limiting equation in [START_REF] Iosifovich | The averaging principle and theorems on large deviations[END_REF] if there are more than one. For example suppose that D ⊂ R m is an open set that contains a single stable attracting point x D of the limiting equation ( 5) and let us consider the first exit time from D

T D = inf{t > 0 : X(t) ∈ D} (15)
Then for any x ∈ D we have

lim α→0 α log E x T D = inf y∈∂D V (x D , y) (16)
where V (x, y) is the quasipotential

V (x, y) = inf T >0 inf γ(0)=x γ(T )=y I T (γ) (17)
While the results above are well-known, in particular expression [START_REF] Nicholson | Introduction to plasma theory[END_REF] for the Hamiltonian in the LDP (see [START_REF] Iosifovich | The averaging principle and theorems on large deviations[END_REF][START_REF] Kifer | Averaging in dynamical systems and large deviations[END_REF][START_REF] Kifer | Averaging principle for fully coupled dynamical systems and large deviations[END_REF][START_REF] Veretennikov | On large deviations for sdes with small diffusion and averaging[END_REF]), little attention has been given to the explicit form this LDP takes via calculation of the Hamiltonian. As stated above, our main aim here is to provide tools to perform such calculations, either analytically or numerically. Even if it goes without saying, we should stress that the Hamiltonian (11) cannot be deduced simply from the limiting equations (5) from the LLN and (10) from the CLT. Indeed, we may naively try to somehow recombine these two limiting equation and write down a nonlinear SDE whose LLN and CLT are precisely ( 5) and [START_REF] Liptser | Large deviations for two scaled diffusions[END_REF]:

since the noise in this SDE would be scaled by a factor √ α one may be tempted to think that its associated LDP would also yield the LDP for the slow variables solution of (1). This program, however, is not achievable in general because there is no SDE for the slow variable alone that admits (5) as LLN and ( 9) as CLT, except in special cases. Indeed the drift term in this equation would have to be F (X) to leading order in α, and setting X = X + √ α ξ and expanding in α would give a linear drift in the CLT equation equivalent to [START_REF] Liptser | Large deviations for two scaled diffusions[END_REF] involving ∂ x F ( X). And there lies the problem:

B(x) = ∂ x F (x) in general.
As we will see below, this problem is consistent with the fact that the LDP for the slow variables X in ( 1) is different in form to that of an SDE for X alone -in particular, the Hamiltonian [START_REF] Nicholson | Introduction to plasma theory[END_REF] involved in this LDP is typically non-quadratic in the momenta conjugate to X, unlike that of an SDE with small noise. This essentially means that the effect of the noise induced by the fast variables on the slow ones cannot be modeled as a Gaussian white-noise in general. Below, we will re-derive the LDP stated above using formal asymptotic expansions tools. These results are complementary to the rigorous results proven e.g. in [START_REF] Iosifovich | The averaging principle and theorems on large deviations[END_REF][START_REF] Kifer | Averaging in dynamical systems and large deviations[END_REF][START_REF] Kifer | Averaging principle for fully coupled dynamical systems and large deviations[END_REF][START_REF] Liptser | Large deviations for two scaled diffusions[END_REF][START_REF] Veretennikov | On large deviations for sdes with small diffusion and averaging[END_REF]. In particular we obtain an equation for the Hamiltonian which, to the best of our knowledge, is new and can be solved explicitly in some nontrivial examples.

B. Organization

The remainder of this paper is organized as follows. In Sec. II we derive the LDP for the fastslow system (1). This derivation is formal, and to connect it with the results stated above we put it within the context of Donsker-Varadhan theory (Sec. II A) of large deviation and Gartner-Ellis theorem (Sec. II B). We also discuss in Sec. II C the link between the LLN and the LDP on the one hand, and the CLT and the LDP on the other. In Sec. III, we specialize the LPD to a class of fast-slow system for the which the equation for the Hamiltonian can be simplified. In Sec. IV we use these results to study a test-case example. Finally, some concluding remarks are given in Sec. V.

The slow-fast system in (1) defines a Markov process with generator L = L 0 + α -1 L 1 where

L 0 = f (x, y) • ∂ x , L 1 = b(x, y) • ∂ y + 1 2 a(x, y) : ∂ y ∂ y (18)
with a(x, y) = (σσ T )(x, y), and a(x, y) : ∂ y ∂ y is the contraction of the two tensors a and ∂ y ∂ y (i.e. the trace of the product a∂ y ∂ y ). Given any suitable test function φ : R m → R the expectation

u(t, x, y) = E x,y φ(X(t)) (19)
therefore satisfies the backward Kolmogorov equation

∂ t u = L 0 u + 1 α L 1 u, u(0) = φ (20)
Our derivation of the LDP for the slow variables X in ( 1) is based on formal asymptotic analysis of this equation. This approach is also at the core of the formal derivations of the LLN equation in ( 5) presented in Appendix A and the CLT equation in [START_REF] Liptser | Large deviations for two scaled diffusions[END_REF] presented in Appendix B.

Since the LDP for the slow variables X in ( 1) is concerned with estimating the probability of having a fluctuation of X away from X that is O(1) in α, and since such probability is expected to be exponentially small in α -1 on time-scales that are O(1) in α, let us consider the expectation

u(t, x, y) = E x,y exp 1 α h(X(t))) , (21)
so that u satisfies (20) for the initial condition

u(0, x, y) = exp 1 α h(x) , (22)
and look for a solution of the type

u(t, x, y) = w(t, x, y) + O(α) exp 1 α S(t, x, y) (23)
where both S and w are assumed to be independent of α. Inserting this ansatz in (20) and collecting terms of increasing power in α, we obtain at leading order, O(α -2 ):

1 2 a : ∂ y S∂ y S = 0 (24)
This equation indicates that S depends on x but not y, S(t, x, y) ≡ S(t, x), or, equivalently,

P S = S ( 25
)
where P is the expectation operator with respect to the invariant measure of the fast virtual process (2), see [START_REF] Dembo | Large deviations techniques and applications[END_REF]. The function S plays the role of the action of the LDP, since combining ( 21) and ( 23) with S(t, x, y) ≡ S(t, x) implies that

S(x, t) = lim α→0 α log E x,y exp 1 α h(X(t))) (26)
which is a variant of Varadhan's Lemma [START_REF] Dembo | Large deviations techniques and applications[END_REF]. The question which remains to be addressed is how to estimate S(x, t)?

To this end, go to next order, O(α -1 ), where the ansatz (23) used in (20) gives

w∂ t S = wL 0 S + L 1 w. ( 27
)
where we used the properties that L 0 is a first order operator, i.e. L 0 (f g) = f L 0 g + gL 0 f , and

∂ y S = L 1 S =
∂ t S = L 0 S + w -1 L 1 w = P L 0 P S + P (w -1 L 1 w) (28)
Inserting (28) back in (27) gives

L 1 w + L 0 Sw = P L 0 P S + P (w -1 L 1 w) w (29)
This equation can be written explicitly as

b(x, y) • ∂ y w + 1 2 a(x, y) : ∂ y ∂ y w + f (x, y) • ∂ x S w = H(x, ∂ x S)w (30)
where we defined

H(x, ∂ x S) = F (x) • ∂ x S + 1 2 R n µ x (dy) a(x, y) : ∂ y log w ∂ y log w (31)
Here we used w -1 L 1 w = L 1 log w + 1 2 a(x, y) : ∂ y log w ∂ y log w and P L 1 = 0, to express the last term at the right hand side of (28) as

P (w -1 L 1 w) = 1 2 P (a(x, y) : ∂ y log w ∂ y log w) (32)
Since x and ∂ x S only enter as parameters in (30), this equation can be viewed an eigenvalue problem for the operator L 1 , where H(x, ∂ x S)-f (x, y)•∂ x S plays the role of (leading) eigenvalue.

The function H(x, ∂ x S) also is the Hamiltonian of the LDP we are trying to derive. Indeed, if we go back to equation (28) for S, we see that it can be expressed in terms of H(x, ∂ x S) as the following Hamilton-Jacobi equation:

∂ t S = H(x, ∂ x S), S(x, 0) = h(x) (33)
This also means that if we introduce the Lagrangian associated with H via

L(x, y) = sup θ (y • θ -H(x, θ)) (34)
then the process X will satisfy a large deviation principle with respect to the action

I T (x) = T 0 L(x(t), ẋ(t))dt (35)
and the solution to (33) can be expressed as

S(t, x) = inf{h + I t (x)} (36)
where the infimum is taken over all paths x : [0, t) → R m such that x(t) = x. Since S(t, x) permits via (26) to evaluate the limit of α time the logarithm of the expectation in (26), this is again Varadhan's Lemma.

To make these results concrete, it remains to see whether we can analyze (30) and get a more convenient equation for H(x, θ). This will be done in Sec. III for a specific class of systems. Before going there, however, we show that our results are consistent with those of Donsker-Varadhan large deviations theory (Sec. II A) as well as with Ellis-Gartner theorem (Sec. II B). We also discuss the link between the LDP we just derived and the LLN on the one hand and the CLT on the other (Sec. II C).

A. Connection with Donsker-Varadhan Theory of Large Deviations

Denoting the solution of (30) for x and ∂ x S = θ fixed by w(y, x, θ), this solution can be expressed as the expectation w(y, x, θ) = lim

T →∞ E y exp -T H(x, θ) + θ • T 0 f (x, Ỹx (τ ))dτ = lim T →∞ e -T H(x,θ) E y exp θ • T 0 f (x, Ỹx (τ ))dτ (37)
where Ỹx (t) is the virtual fast process, solution to (2), and H(x, θ) needs to be adjusted to make the limit converge to a finite, nonzero value. The requirement that such an adjustment be possible also gives the condition for the solution to (30) to exist and be uniquely given by (37). By taking the logarithm of the factor under the limit in (37), it is easy to see that this requirement imposes as a necessary condition that the following limit exists

H(x, θ) = lim T →∞ 1 T log E y exp θ • T 0 f (x, Ỹx (τ ))dτ (38)
This is the alternative expression for H already given in [START_REF] Nicholson | Introduction to plasma theory[END_REF] that is consistent with the one in Donsker-Varadhan large deviations theory.

B. Connection with Gartner-Ellis Theorem

Gartner-Ellis Theorem tells us that the Hamiltonian in (31) can also be defined by taking the limit as α → 0 of

α log E exp 1 α T 0 θ(t) • f (x(t), Y x (t))dt (39)
where Y x (t) is the solution of the second equation in (1) for a given x(•), i.e.

dY x = 1 α b(x(t), Y x )dt + 1 √ α σ(x(t), Y x )dW (t) (40)
Letting

û(t, y) = E y exp 1 α t 0 θ(t ) • F (x(t ), Y x (t ))dt (41)
this function satisfies

∂ t û = 1 α L 1 û + 1 α θ(t) • f (x(t), y) û (42)
Look for a solution of the type

û(t, y) = ŵ(t, y) exp 1 α φ(t) (43)
Proceeding similarly as above then lead to the following system of equations for φ and ŵ (compare (30) and (33)):

b(x, y) • ∂ y ŵ + 1 2 a(x, y) : ∂ y ∂ y ŵ + f (x, y) • θ ŵ = H(x, θ) ŵ (44)
and

∂ t φ = H(x, θ) (45)
for the Hamiltonian H defined in (31). This also means that

lim α→0 α log E exp 1 α T 0 θ(t) • f (x(t), Y x (t))dt = T 0 H(x(t), θ(t))dt (46)
which leads again to the action in (35).

C. Link between the LDP, the LLN, and the CLT Suppose that we expand the exponential in (38) to second order in θ, take the expectation, then expand the logarithm to second order in θ as well, and finally take the limit as T → ∞. This sequence of operations corresponds to making a cumulant expansion of the variable θ • T 0 f (x, Ỹx (τ ))dτ truncated to second order, and it gives the following quadratic approximation for H:

H quad (x, θ) = θ • F (x) + 1 2 θ T A(x)θ (47)
where F (x) is defined in (3) and A(x) in [START_REF] Heymann | Pathways of maximum likelihood for rare events in nonequilibrium systems: Application to nucleation in the presence of shear[END_REF]. This is the Hamiltonian for the LDP associated with the SDE

dX = F (X)dt + √ α η(X)dW (t) (48)
in the limit as α → 0. The process defined by (48) satisfies a LLN principle with limiting equation [START_REF] Iosifovich | The averaging principle and theorems on large deviations[END_REF], meaning that the LDP contains the information about the LLN. Yet, it also highlights the subtle (and well-know) differences between the CLT and the LDP. Indeed, the process defined by (48) satisfies a CLT with respect to

dξ = ∂ x F ( X)ξdt + η( X)dW (t) (49)
This equation is not identical with (10) -their drift terms are different. In fact, it is easy to see that we would have to add deterministic terms of order O( √ α) in (48) in order that the CLT associated with this modified equation coincide with (10) -as already mentioned in the introduction, the addition of these terms would render the resulting equation unclosed since they depend on ξ rather than X. These additional terms do not affect the Hamiltonian of the LDP which would still be (47).

This goes to show that the CLT cannot be recovered from the LDP: in the range of values for X

where it applies, it contains finer information than that in the LDP.

Conversely, the LDP cannot be deduced from the CLT, as the actual H in (38) is different in general from its quadratic approximation H quad . More precisely, H = H quad iff θ• T 0 f (x, Ỹx (τ ))dτ is Gaussian, and Marcinkiewicz's theorem states that in all other cases H is not a polynomial of θ of any order (i.e. its expansion in θ involves infinitely many terms). As we will see below in Sec. III, H quad = H if the fast Y enter linearly the equation for the slow X. When this is not the case and H quad = H, H quad can be used to describe moderate fluctuations of order O(α ν ) with 1 2 ≤ ν < 1, that is, outside the range of validity of the CLT but only moderately so. To describe large fluctuations of order O(1), however, we need to use the actual H in (38).

III. A SPECIFIC CLASS OF SYSTEMS

Next let us specialize (1) to systems in which the dynamics of the fast Y is linear and they enter quadratically the equation for the slow X

     Ẋ = r(X) + Y T s(X) + Y T M (X)Y dY = - 1 α L(X)Y + 1 √ α σ(X)dW (50)
where X ∈ R (the generalization to the vectorial case is straightforward but it makes notations more cumbersome, so we will stick to the scalar case here -see however Sec. IV C for an illustration with X ∈ R 2 ), Y ∈ R n , r : R → R (e.g. r(x) = -νx for some ν > 0), s : R → R n , L : R → R n×n is a positive-definite matrix, M : R → R n×n is a symmetric matrix, and σ : R → R n×n . For such systems, the equation (30) leading to the Hamiltonian of the LDP can be written down more explicitly.

To see how, notice first that the virtual fast process defined in ( 2) is given explicitly by

Ỹ x (τ ) = e -L(x)τ y + τ 0 e -L(x)(τ -τ ) σ(x)dW (τ ) (51)
As a result

µ x (dy) = (2π) -n/2 (det C(x)) -1/2 exp - 1 2 y T C -1 (x)y dy ( 52
)
where C(x) is the equilibrium covariance matrix of Ỹ x (τ ) satisfying the Lyapunov equation

L(x)C(x) + C(x)L T (x) = a(x) (53)
where a(x) = (σσ T )(x). This means that the limiting equation ( 5) from the LLN reads

Ẋ = r( X) + tr C( X)M ( X) (54)
and the equation ( 10) from the CLT reads

dξ = r ( X)ξdt + tr C( X)M ( X) ξdt + g( X)ξdt + η( X)dW (t) (55)
Here

g(x) = - ∞ 0 tr C(x)(L (x) + [L (x)] T )e -L T (x)τ M (x)e -L(x)τ dτ (56)
and

η 2 (x) = ∞ 0 s T (x) C(x)e -L T (x)τ + e -L(x)τ C(x) s(x)dτ + 4 ∞ 0 tr C(x)e -L(x)τ M (x)e -L T (x)τ C(x)M (x) dτ = s T (x)L -1 (x)a(x)L -T (x)s(x) + 4 ∞ 0 tr C(x)e -L(x)τ M (x)e -L T (x)τ C(x)M (x) dτ (57)
where we used the Lyapunov equation ( 53).

Turning ourselves to the LDP next, (30) is explicitly

-L(x)y • ∂ y w + 1 2 a(x) : ∂ y ∂ y w + r(x) + y T • s(x) + y T M (x)y θ w = H(x, θ)w (58)
where θ = ∂ x S and

H(x, θ) = (r(x) + tr(C(x)M (x))) θ + 1 2 µ x (dy) a(x) : ∂ y log w ∂ y log w (59)
Look for a solution of the form

w = exp(y T m(x, θ) + y T N (x, θ)y) (60)
for some unknown m(x, θ) and N (x, θ). Then (dropping the dependencies in x for simplicity of notation)

H(θ) = (r + tr(CM )) θ + 1 2 m T (θ)am(θ) + 2 tr (CN (θ)aN (θ)) (61)
and (58) becomes

-y T N (θ)Ly -y T L T N (θ)y -y T L T m(θ) + tr(aN (θ)) + 2y T N (θ)am(θ) + 2y T N (θ)aN (θ)y + (y T s + y T M y)θ = tr(CM )θ + 2 tr (aN (θ)CN (θ)) (62)
Collecting the terms that are of order 0, 1, and 2 in y, respectively, gives the equations tr(aN (θ)) = tr(CM )θ + 2 tr (CN (θ)aN (θ)) (63)

(L T -2N (θ)a)m(θ) = s θ (64) N (θ)L + L T N (θ) = 2N (θ)aN (θ) + M θ (65)
By right multiplying (65) by C, taking the trace, and using the Lyapunov equation ( 53), it is easy to see that the result is (63), i.e. if (65) is satisfied then (63) automatically holds. We can also solve (64) in m(θ) to get

m(θ) = (L T -2N (θ)a) -1 s θ (66)
and be left with solving (65) in N (θ) -note that this equation may have more than one solution, and we should take the one such that N (0) = 0, so that H(0) = 0. Inserting ( 63) and ( 66) in ( 61)

we then obtain the Hamiltonian of the LDP in terms of N (θ) alone

H(θ) = rθ + tr (aN (θ)) + 1 2 s T (L T -2N (θ)a) -T a(L T -2N (θ)a) -1 s θ 2 (67)
The solution to (65) is not available explicitly in general. There is one trivial case, however, namely when M = 0. In this case it is easy to see that N (θ) = 0, m(θ) = L -T s θ and the Hamiltonian is quadratic

H(θ) = rθ + 1 2 s T L -1 aL -T s θ 2 (68)
This Hamiltonian is that of the LDP associated with SDE

dX = r(X)dt + √ αη 0 (X)dW (t) (69)
where η 0 (x) is the factor defined in (57) evaluated at M = 0. It is easy to see that this the limiting equation from the LLN for this equation is (54) (with M = 0) and the equation from the CLT is (55) (again with M = 0). Thus, if the fast variables are Gaussian, and their action on the slow one is linear, the LDP contains all the information about the CLT, and the SDE (69) can be used to investigate large deviations. Notice that this includes nontrivial situations with metastability, when ẋ = r(x) has more than one stable fixed point and one is interested in the rate and mechanism of transitions between these points.

Another case where the Hamiltonian can be computed explicitly is the following one: assume that a is invertible, and that the following conditions hold:

La = aL T and L T M a = M aL T . (70)
Then, it is straightforward to check that whenever B(θ) = L T 2 -2θM a admits a square root, the matrix

N (θ) = 1 2 L T -B(θ) a -1 (71)
satisfies N (θ)L = L T N (θ), and is a solution of (65). Using (70), we also have that aB(θ) is symmetric, so inverting (71) and using (65), we prove that N (θ) is symmetric, which is consistent with its definition (60). It is then also easy to prove that a B(θ) -1 = B T (θ) -1 a, so the Hamiltonian (67) reads

H(θ) = rθ + 1 2 tr L - √ L 2 -2θaM + 1 2 s T L 2 -2θaM -1 as θ 2 , (72)
whenever the square roots and inverses appearing in this equation exist (which is the case for θ = 0). Note that the square root in the trace should be chosen such that H(θ = 0) = 0.

IV. A CASE STUDY

In this section, we illustrate our results on the following test case example:

         Ẋ = 1 K K k=1 Y 2 k -νX dY k = - 1 α γ(X)Y k dt + σ √ α dW k , k = 1, . . . , K (73)
where W k are independent Wiener processes, α > 0, ν > 0 and σ are parameters and γ(X) > 0 is a function to be specified later. (73) consists of a scalar X (the vectorial case is discussed below in Sec. IV C) coupled to K Ornstein-Uhlenbeck processes Y k which each feel an independent noise and whose common decay rate depends on X. (73) is in the class of (50) with

r(x) = -νx, s(x) = 0, M (x) = K -1 Id, L(x) = γ(x) Id, σ(x) = σ Id (74)
We will be interested in studying (73) in the limits α → 0 and K → ∞ -the former limit is in the realm of the formalism developed here, whereas the latter can be estimated by direct calculation.

As we will below these two limits commute.

A. The limit as α → 0

Using the formulas given in section III, it is easy to see that the LLN equation (54) becomes

Ẋ = σ 2 2γ( X) -ν X (75)
and the CLT equation ( 55) is

dξ = -ν + σ 2 γ ( X) 2γ 2 ( X) ξdt + σ 2 2Kγ 3 ( X) dW (t) (76)
As far as the LDP is concerned, note that (73) is of the form

Ẋ = 1 K K k=1 f 0 (X, Y k ) + f 1 (X) (77)
where Y k are i.i.d. random processes. From (38), we see that the Hamiltonian reads

H(x, θ) = θf 1 (x) + KH 0 x, θ K (78)
with H 0 the one-particle Hamiltonian

H 0 (x, θ) = lim T →∞ 1 T log E y 1 exp θ • T 0 f 0 (x, Ỹ1,x (τ ))dτ . (79)
To get the Lagrangian associated with this Hamiltonian, we solve

ẋ = ∂H ∂θ = f 1 (x) + ∂H 0 ∂θ x, θ K (80)
and get θ/K as a function of (x, ẋ). Denoting this solution as θ = Kϑ(x, ẋ) and using ( 78) and

(80), we then deduce that the Lagrangian is proportional to K and reads

L(x, ẋ) = ẋθ -H(x, θ) = K ϑ(x, ẋ) ∂H 0 ∂θ (x, ϑ(x, ẋ)) -H 0 (x, ϑ(x, ẋ)) . (81)
As a result, the path that minimizes the action (35) (instanton) does not depend on K, but the probabilistic weight of this path decreases exponentially with K. The quasi-potential V (x), which by definition is the solution of

H(x, ∂ x V ) = 0 (82)
with H given by (78), is then also proportional to K.

Going back to the particular case (73), the one-particle Hamiltonian (79) can be computed using (72). Indeed, we are now solving a one-dimensional problem, so the conditions (70) are fulfilled and ( 72) is simply

H 0 (x, θ) = 1 2 γ(x) -γ 2 (x) -2σ 2 θ , (83)
defined whenever θ γ 2 (x)/2σ 2 . To get the Lagrangian associated with this Hamiltonian, we solve

ẋ + νx = ∂H 0 ∂θ x, θ K = 1 2 σ 2 γ 2 (x) -2σ 2 θ/K . ( 84
) ϑ(x, ẋ) = θ K = γ 2 (x) 2σ 2 - σ 2 8( ẋ + νx) 2 (85)
and get

L(x, ẋ) = K 8σ 2 |2γ(x)( ẋ + νx) -σ 2 | 2 ẋ + νx (86)
whenever ẋ + νx > 0.

The quasi-potential is given by (82), using (78) and (83) it is easy to see that this implies that either

∂ x V ≡ V = 0, or V (x) K = νxγ(x) -1 2 σ 2 ν 2 x 2 (87)
This result should be compared with the one obtained from the quadratic approximation to H,

H quad = -νxθ + K σ 2 θ/K 2γ(x) + σ 4 (θ/K) 2 4γ 3 (x) , (88)
for which we deduce

V quad (x) K = 4γ 3 (x) σ 4 νx - σ 2 2γ(x) (89)
The potentials V (x) and V quad (x) are different in general. To give a concrete example, consider the case γ(x) = x 4 /10x 2 + 3, ν = 1, σ = √ 3, which leads to bistability of the slow process

x(t). The potentials V and V quad are represented in figure 1. The extrema of V quad are also extrema of V , and the second derivatives of these two potentials are the same at these extrema, as should be expected. Figure 2 illustrates this last point for the potential minima. However, we see that the global shape of the potentials are very different, and in particular the energy barrier between the attractors of V quad is almost twice the one between the attractors of V , which means that the probability of rare transitions obtained from the quadratic approximation in this case will be much lower than the actual one. it is easy to see that (73) can be rewritten as

     Ẋ = E -νX dE = - 2 α γ(X)Edt + 1 α σ 2 dt + 2 σ √ αK √ E dW (91)
where we used the identity

1 K K k=1 Y k dW k = K -1/2 √ E dW (in law) (92)
Since the noise term in (91) is small when K is large, we can use large deviation theory to analyze the behavior of the system in the limit as K → ∞. The Freidlin-Wentzell action associated with (91) reads

I α T (E, x) = α 8σ 2 T 0 | Ė + 2γ(x)E/α -σ 2 /α| 2 E dt (93)
if E = ẋ + νx and I T (E, x) = ∞ otherwise. Letting α → 0 and keeping only the leading order term gives

I α T (E, x) ∼ 1 8σ 2 α T 0 |2γ(x)( ẋ + νx) -σ 2 | 2 ẋ + νx dt (94)
which is consistent with (86) since they both imply that the probability weight on paths is roughly

exp - K 8σ 2 α T 0 |2γ(x)( ẋ + νx) -σ 2 | 2 ẋ + νx dt ( 95
)
when α is small and/or K is large.

C. Two-dimensional generalization

To illustrate the impact that the non-quadratic nature of the Hamiltonian has on the pathway of the transition, let us now consider the following generalization of (73):

     Ẋi = -β i X i + κ 2 j=1 D ij X j + Y 2 i , i = 1, 2 dY i = - 1 α γ(X i )Y i dt + 1 √ α σdW i , i = 1, 2 (96)
with

D 11 = D 22 = -1, D 12 = D 21 = 1 (97)
The LLN equations for the system (96) are given by

Ẋi = σ 2 2γ( Xi ) -β i Xi + κ 2 j=1 D ij X j , i = 1, 2. (98)
For the specific choice γ(x) = (x -5) 2 + 1, and β 1 = 0.6, β 2 = 0.3 and σ = √ 10, the flow field associated with (98) is shown in figure 3: it has the two stable fixed points (shown as red circles in the figure) with one unstable critical point (shown as a red square) in between. The Hamiltonian associated with (96) can be written as

H(x, θ) = 2 i=1 h(x i , θ i ) - 2 i=1 β i x i θ i + κ 2 i,j=1 θ j D ij x j (99) with h(x, θ) = 1 2 γ(x) -γ 2 (x) -2σ 2 θ . (100)
A numerical computation of the transition trajectories between the two stable fixed points was performed using the geometric minimum action method GMAM [START_REF] Heymann | The geometric minimum action method: A least action principle on the space of curves[END_REF][START_REF] Heymann | Pathways of maximum likelihood for rare events in nonequilibrium systems: Application to nucleation in the presence of shear[END_REF] (building on the method introduced in [START_REF] Weinan | Minimum action method for the study of rare events[END_REF]) and is shown in figure 3: The red solid line depicts the minimizing trajectory for the full Hamiltonian (99), while the dashed line represents the minimizing trajectory for a quadratic approximation of H(x, θ), with clear differences between the two. The respective probability of the minimizers can be seen by evaluating the action along the trajectories, as shown in Fig. 4: The action evaluated along its minimizer is lower than evaluated along the minimizer of the quadratic action. Also showed in dashed lines are the relaxation paths from the unstable to the stable equilibrium points from the limiting equation ( 98): These paths are followed by the mini- mizing trajectories on the way down from the unstable critical points (as they should: no noise is necessary for this part of the transition paths) but not on the way uphill.

V. CONCLUDING REMARKS

We have investigated how large deviations affect the slow variables in fast-slow systems via analysis of the large deviation principle (LDP) that characterize their likelihood and pathways.

For a specific class of systems, we derived an algebraic equation for the Hamiltonian involved in this LDP, and we discussed several situations in which this equation can be solved explicitly.

These results show that the way rare events or infrequent transitions arise in fast-slow systems is intrinsically different from the way they would arise if the dynamics of the slow variables was approximated by an SDE -these difference stem from the fact that the Hamiltonian is non-quadratic in the momenta in general. The examples treated in the present paper were simple enough to allow for analytic treatment. However, we believe that our results will be useful in more complicated where B 1 (x) is

B 1 (x) = R n ∂ x f (x, y)µ x (dy) (B10)
The second term at the right hand side of (B8) is

P L 1 2 ∞ 0 dτ e τ L 1 L1 2 P v 0 = B 2 (x)ξ • ∂ ξ v 0 + A(x) : ∂ ξ ∂ ξ v 0 (B11)
where A(x) is the matrix defined in [START_REF] Heymann | Pathways of maximum likelihood for rare events in nonequilibrium systems: Application to nucleation in the presence of shear[END_REF] and

B 2 (x) = ∞ 0 dτ R n ∂ y E y f (x, Ỹ x (τ )) ∂ x b(x, y)µ x (dy) (B12)
Inserting (B9) and (B11) in (B8) shows that this equation is indeed the backward Kolmogorov equation of the joint process governed by ( 5) and [START_REF] Liptser | Large deviations for two scaled diffusions[END_REF].

B.FIG. 1 : 3 .FIG. 2 :

 132 FIG.1:The potential V (x)/K and the one obtained from a quadratic approximation of the HamiltonianV quad (x)/K, see (89), for γ(x) = x 4 /10x 2 + 3, ν = 1, σ = √ 3.The quadratic potential obtained from the quadratic approximation is quite different from the actual potential.

FIG. 3 :

 3 FIG. 3: Dynamics of coupled slow-fast system ODE model for κ = 0.2. The arrows denote the direction of the deterministic flow, the color its magnitude. The red solid line depicts the minimizer of the actual action associated with the Hamiltonian in (99), the red dashed line the minimizer of the quadratic approximation of this Hamiltonian, and the red dotted line the relaxation paths from the saddle via the limiting equation (98). Red markers are located at the fixed points (circle: stable; square: saddle).

2 FIG. 4 :

 24 FIG.4: Action density dI (variation of the action per unit length) versus path length s, for the coupled slowfast system ODE model in (96) for κ = 0.2, for paths up to the saddle point. The action density is computed with respect to the full Hamiltonian (99) for the three trajectories depicted in figure3, i.e. the minimizer of the actual action, the minimizer of the quadratic approximation of this action, and the relaxation pathway from the limiting equation (98).

Acknowledgments E. V.-E. thank David Kelly for interesting discussions. The research leading to these results has received funding from the European Research Council under the European Union's seventh Framework Programme (FP7/2007-2013 Grant Agreement no. 616811) (F. Bouchet, and T. Tangarife).

Appendix A: Derivation of the limiting equation [START_REF] Iosifovich | The averaging principle and theorems on large deviations[END_REF] from the LLN Here we derive the limiting equation ( 5) of the LLN by formally taking the limit as α → 0 on the backward Kolmogorov equation (20). To this end expand u as u = u 0 + αu 1 + O(α 2 ), insert this ansatz in (20), and collect term of increasing power in α. This gives the hierarchy

The first implies that u 0 is a only a function of x and not of y, or equivalently

Since L 1 is not invertible (P L 1 = 0), the second equation requires a solvability condition, which reads

It is easy to see that P L 0 P = F (x) • ∂ x , i.e. (A3) is the backward Kolmogorov equation of the limiting ODE [START_REF] Iosifovich | The averaging principle and theorems on large deviations[END_REF].

Appendix B: Derivation of the CLT equation [START_REF] Liptser | Large deviations for two scaled diffusions[END_REF] To derive the linear SDE (10) of the CLT, notice that, using ( 5), (1) can be rewritten as

This means that the joint process ( X, ξ, Y ) is Markov with generator

where L 1 is defined in (18) and

this function satisfies the backward Kolmogorov equation

), insert this ansatz in (B4), and collect term of increasing power in α:

The first equation implies that v 0 = P v 0 , i.e. v 0 is a function of x and ξ only. The solvability condition for the second equation is automatically satisfied since P f = 0 implies that P L 1 2 P = 0.

Therefore, the solution to this equation is

where L -1 1 denotes the pseudo-inverse of L 1 . Alternatively, this solution can also be expressed as

Using this expression in the solvability condition for the third equation in (B5) finally gives the evolution equation for v 0 :

The first term at the right hand side is explicitly


Jean-Philippe Chancelier

Ramine Nikoukhah

A novel code generation methodology for block diagram modeler and simulators Scicos and VSS

Block operations during simulation in Scicos and VSS environments can naturally be described as Nsp functions. But the direct use of Nsp functions for simulation leads to poor performance since the Nsp language is interpreted, not compiled. The methodology presented in this paper is used to develop a tool for generating efficient compilable code, such as C and ADA, for Scicos and VSS models from these block Nsp functions. Operator overloading and partial evaluation are the key elements of this novel approach. This methodology may be used in other simulation environments such as Matlab/Simulink.

Introduction

In this paper we present a new methodology of code generation for block diagram simulation tools such as Scicos [START_REF] Campbell | Modeling and Simulation in Scilab/Scicos with ScicosLab 4[END_REF], VisSim SIMULATE1 (VSS) and Simulink. This methodology takes advantage of the operator overloading facilities available in matrix based languages such as Matlab, Scilab, Nsp [START_REF] Nsp | [END_REF] or Octave. The implementation is done in Nsp/Scicos and VSS environments [START_REF] Campbell | Modeling and Simulation in Scilab/Scicos with ScicosLab 4[END_REF] but the code may easily be ported to Matlab/Simulink and other similar environments. In our implementation, the final C code can be either generated directly using a C pretty printer or through the use of Gnat Model Compiler [START_REF] Adacore | GNAT Pro MOdel Compiler User Guide[END_REF]. In the latter case, an .xmi file is generated and then gmc is invoked to generate the C code. Either way, the C code may be reimported automatically in the Scicos and VSS environments as a CBlock for validation.

Not only the block behaviors may easily and naturally be implemented in the Nsp language but they can also be simulated directly in Scicos, if implemented as such. The implementation is natural because most Scicos blocks support matrices and Nsp data types. But also because Nsp is the working language in the Nsp/Scicos environment. In particular Scicos model and block parameters are defined via Nsp scripts and expressions, and evaluated by the Nsp interpreter.

The main reason for which only parameter definitions in Scicos blocks use Nsp programs and run-time block code is invariably expressed in C is performance. Parameters are evaluated once at compile time, so peformance is not an issue, but if the run time block behavior is expressed in Nsp, the simulation performance would be significantly reduced.

The methodology proposed in this paper provides the possibility of defining Scicos/VSS blocks entirely in Nsp. These blocks, which can be developed by Scicos developers, toolbox developers or end-users, may be tested and validated directly by (low performance because interpreted) simulations in Scicos. The Nsp code of the block can then be used for code generation by the methodology presented in this paper.

In this methodology, an Nsp script is generated from the Scicos/VSS model. The generation of this script is based on the result of the compilation of the model by Scicos or VSS compilers. The compiler result includes in particular the order in which the blocks in the model must be executed (their functions called) and the types and sizes of all the block inputs and outputs. This generated script, which contains calls to the Nsp functions of the blocks, evaluates block outputs, new discrete states, state derivatives, zero-crossing surfaces, etc., i.e., all the information needed for both simulation and code generation.

The main obstacle in generating compilable code (such as C or ADA) from Nsp, and more generally from all Matlab-like languages, is the inability, in the general case, to determine the types and sizes of variables before run-time. But this is not an issue for the script generated from a Scicos or VSS model since in this case, all the sizes and types of the variables representing block inputs and outputs are known (determined by the compiler). The main component of the code generation tool presented here is the Nsp code generation facility that can produce compilable code from such Nsp scripts. It should be emphasized that this facility is not a general purpose code generator for arbitray Nsp scripts. It is developed especially for Scicos/VSS usage where in the generated code, data types involved are limited to matrices of doubles and various integer data types, and all variable sizes can be statically determined. This Nsp code generation facility is entirely developed in Nsp using its abilities to define new data types and performing operator overloading for these new data types. In particular we define a new data type called bvar representing a matrix and containing:

• type: numerics or symbolics

• value: a matrix • name: string If a variable of type bvar is of type numerics, then its value represents its actual value and it is similar to a regular Nsp matrix with the same value. However using the type bvar, it is possible to associate a name to the variable. If the variable is of type symbolics, then its value is just a "nominal value", often zero, but its data type (double, int8, int16,...) and its size represent the actual data type and size of the variable.

Operations involving numeric variables (of type bvar) and regular Nsp variables never produce symbolic variables (of type bvar). Operations involving symbolic variables produce in general symbolic variables except in special cases. In particular if X is a symbolic variable, then size(X) is not; neither is datatype(X). The overloading of basic operators and language primitives for these data types automatically overloads all Nsp functions that use these operators and primitives thanks to dynamic scoping property of the Nsp language.

The main idea of the approach used in this paper is to overload basic Nsp operators and functions for variables of type bvar so that if we have an Nsp script that performs operations on regular matrices, if these matrices are replaced by variables of type bvar, the script runs as before but in addition, a pseudo-code (a sequence of "instructions") is generated and stored in a global list variable named code. Every time an operation involves a variable of type bvar, the overloaded operation performs the original operation returning the corresponding value as if the arguments were not of type bvar but in addition it adds one or several "instructions" to the list code. This mechanism may be seen as spying the script as it is executed: when an operation involving a type bvar variable is performed, it is recorded in the code list. The record includes not only the nature of the operation (for example matrix multiplication) but also certain information on its operands. The information in the generated pseudo-code is rich enough so that it can be used to generate code in C, ADA or other languages, by simple pretty printing operations.

The execution of the script can be seen as partial evaluation. When the value is unknown, only the type and size information is propagated; when the value is known, the actual value is propagated. As long as the script can run to the end by propagating this partial information, the code generation can be performed.

Unlike most code generators where the objective is to transform a program written in a source language into a program in the target language preserving its functionality, in our approach the program in the source language contains not only the code to be translated but also specifications on how this transformation should be performed. This is done by providing a set of functions for use as code generation directives so that a code generation task may be expressed as a part of the Nsp script. The Scicos/VSS block Nsp codes however will not make use of these functions. Scicos users need not know about these directives; as far as they are concerned, the block codes are written in standard Nsp language with regular Nsp data. During code generation, the script will be calling these block codes with the variables expressing blocks input/outputs and states defined as symbolic bvar variables unbeknownst to the Scicos user.

There are many advantages in using the technique presented here as opposed to developing an independent code generator: this code generator does not use a different parser for the scripting language; the script is run using the tool itself. Similarly, the operator and functions used in the process of partial evaluation are that of the scripting language itself. Finally the operator overloading is being performed in the scripting language, so it can be customized easily for different targets.

The method presented here is implemented in Nsp but it may be developed in other languages providing similar overloading facilities, in particular Matlab, Scilab and Octave.

The tool based on this methodology provides a powerful and easily customizable code generator for Scicos and VSS. Despite the existence of code generators for Simulink (Matlab and Simulink coder in particular) a tool based on the approach presented in this paper would still be of interest since it provides more control over the generated code.

In this paper we will focus on code generation for discrete-time subsystems of Scicos and VSS models. This is what is needed in the majority of embedded code generation applications, in particular for the implementation of embedded controllers. We consider double, boolean and signed and unsigned 8, 16 and 32 bit integers data types. We do not consider the complex data type (seldom used in embedded code), nor do we consider fixed-point data types, which will be considered in a subsequent paper and for which the methodology presented here is very promising.

Simple Example

We start with a very simple example to illustrate the basic idea of our approach. Consider the Scicos model depicted in Fig. 1. The objective is to generate C code for the Super Block the contents of which is depicted in Fig. 2. Scicos code generator produces an Nsp script, based on the result of the Scicos compiler, corresponding to the subsystem inside the Super Block. This script calls the Nsp functions associated with the blocks in the subsystem. The block functions have an imposed calling sequence and provide specific callback functionalities. Here we review the functions associated with the basic blocks used in this subsystem to illustrate how block functions may be developed for code generation.

A simplified version of the Nsp code of the delay block 1/z is given below: The arguments of the block functions are an Nsp structure and an integer. The structure, named block in this code, contains the fields The value -1 is used for initialization. The function is called with this flag once at the beginning so that it can initialize the internal states of the block, if any. Flag values 1 and 2 are used respectively to code the output update and state update phases of the block behavior. Here, we consider only code generation for discrete-time subsystems. For more general blocks, other flag values are used for example for the computation of state derivative, evaluation of zero-crossing surfaces, etc. For more information on the usage of flags see [START_REF] Campbell | Modeling and Simulation in Scilab/Scicos with ScicosLab 4[END_REF], Chapter 9.

In the script generated from the Scicos or VSS model, the Nsp code of the blocks are called several times with different flag values. Unlike block parameters, which are known at the time of code generation, block input/output and state values are not known. Their values are defined as bvar variables when the block function is called. Even though the function has been developed to receive a structure containting regular matrix arguments, thanks to operator overloading, the function accepts this new argument.

The structures io and state are also new data types for which insertion and extraction operations have been overloaded. But once again, the developer of the block function code needs not know about it.

When this function is called with flag -1, normally we would simply copy the parameter value in the state, but the Delay block accepts that the parameter be defined as a scalar, even if the state is not scalar, in which case all the entries of the state matrix take the value of the parameter. Moreover, even if the parameter is defined as a double (default data type in Nsp), the block may be of a different data type so the parameter should be converted into correct data type before being copied in the state. Note that the expressions datatype(block.io(1)) and size(block.io( 1)) in this code can be evaluated even if block.io( 1) is symbolic (their values depend only on the data type and size of block.io(1)). Needless to say that the execution of this part of the code does not perform any code generation, and thus no trace of this code will be seen in the final generated code.

For this block, the state update and output update operations (Flag values 1 and 2) consist simply of copying respectively the input to the state, and the state to the output. The corresponding pseudo-code is generated thanks to the overloading of the insertion operation in the io and state fields.

The simplified version of the code for the Gain block is as follows: This block does not have a state, so its function contains no initialization phase or state update phase. As in the case of the previous block, the block parameter, the gain value in this case, may be defined as a double even if the block data type is something else, for example any type of integer.

The Gain block in Scicos may perform different types of operations. The gain parameter may be a matrix and the block output is the product of this matrix with the input vector (which must have compatible sizes). But the gain parameter may also be a scalar, in which case the input and output have identical sizes and the product multiplies all the elements of the input matrix with this parameter to compute the output. But the parameter may also be a matrix when the input is a scalar. But all of this complexity is captured in the Nsp multiplication operation used in the code. Indeed the Nsp multiplication operation has exactly the same properties. This similarity is often seen between basic Nsp operators and primitives, and basic Scicos and VSS block behaviors. This is the fundamental reason why using Nsp coding is a convenient way of describing the behavior of Scicos and VSS blocks, often leading to a compact, readable and maintainable Nsp code.

The function put_annotation is used to place comments in the generated code. We now give the complete code of the SUMMATION block function because it shows some advanced features of block definitions that can directly be handled by our approach: This block performs two very different summing operations. When the block has only one input, the output of the block is the sum of the entries of the input matrix. When it has more than one inputs, the block performs the addition (or substraction depending on the value of parameter p2) of the input matrices. The Nsp code contains both algorithms and the selection is made by testing the number of inputs. This test and other tests in particular for deciding if the operation to perform is an addition or a subtraction are "partial evaluated" and will not appear in the generated code.

Note that the conditional expressions of the if statements and the range of for loops must be known (through partial evaluation) at the time of code generation. For example in this case, nin, which represents the number of inputs, is known because the length of block.io is known. The generated code will of course contain no trace of these conditional statements or loops. Clearly if the conditional expression of an if statement is not known, the execution of the script cannot continue. This is a limitation of this approach.

Note also that the error functions generate errors at the time of code generation if a parameter value is incorrect. But if the Nsp code is used as run time code for simulation, they generate error messages at the beginning of the simulation if a parameter does not have correct value.

The Nsp code for the Scicos MUX block is also straightforward since the operation of this block corresponds to row concatenation in Nsp. The block may have more than two inputs in which case repeated concatenations are required within a for loop: The generated C code for the Super Block in Fig. 2 The calling sequence of the state update and output update C routines of the block contain pointers to the input and output data, in that order. The "state" of the block, which consists of the subsystem blocks' input/outputs (link values), and the state of the Delay block are defined as global variables so that their values are preserved from the output update phase to state update phase, and from one time step to the next. Not all link values need to be treated as global, further code optimization may be used to identify link values that need not be preserved and thus simplify the code by removing these links out of the global variables. In this case additional optimization will also remove the unnecessary temporary variables tmp_1 and tmp_2, but the C compiler does this anyway.

Even though there is no scalar type in Nsp (a scalar is just a one by one matrix) the generated C code uses, for example for real values, the type double for scalar variables and not double *, which is used for non-scalars only. But the arguments of function calls are all passed as pointers.

Following the code generation process, the Super Block is automatically replaced in the model with a CBlock containing the generated C code; see Fig. 3. The simulation of this new model can be used to compare simulation results before and after code generation. This is a fast and efficient method for testing the generated code.

Operator overloading and pseudo-code generation

Overloading of operators and primitives is a powerful method in Matlab like environments used for various applications, such as automatic differentiation [START_REF] Weinstein | A source transformation via operator overloading method for automatic differentiation in matlab[END_REF][START_REF] Patterson | An efficient overloaded method for computing derivatives of mathematical functions in matlab[END_REF]. Overloading can be used to alter the behavior of an existing code, in particular by adding additional operations such as the evaluation of the derivatives for automatic differentiation, or pseudo-code generation in our case.

New data types and overloading

In Nsp [START_REF] Nsp | [END_REF], defining new datatypes for which operator overloading is possible is done by writing the new datatype definition in a C code. Some Nsp internal tools are available to ease the writing of the code necessary to implement a new datatype definition. As an example, the C-code needed to implement the bvar datatype was generated from the following simple definition of the bvar object.

(define-object Bvar (in-module "Bvar") (parent "Object") (c-name "NspBvar") (fields '("gboolean" "sym" "hidden" "FALSE"); '("NspObject*" "value" "hidden" "NULL" ); '("char*" "varname" "hidden" ); ) (gtype-id "Bvar") )

When the new datatype is defined, the definition of the overloaded operators can be performed in Nsp or in C. Examples on how to define new data types in Nsp are given in the src/types-test directory of Nsp source code.

The main new data type used in this application is called bvar (block variable). An object of type bvar is a wrapper around a standard Nsp object decorated with new attributes: an attribute containing the name of the variable and an attribute to specify whether the object is numeric or symbolic. The following lines of code explain how to create a bvar variable, which encapsulates a 1x1 symbolic matrix of type double named x having nominal value 67.

-nsp->A=bvar(varname="x",value=67,symbolic=%t) A = "x",%t,Mat The fields in the A variable can be obtained or changed using methods, and the type of object A is obtained by using the type function. The short type name of a variable is the string which is used to define functions used for overloading. Suppose we define the function f_bvar in a function library as follows : Nsp operators and functions used for writing Scicos and VSS blocks are overloaded for bvar variables. As an example we give the code for the function inv_bvar, which implements matrix inversion. When the matrix to be inverted is 2x2 we use the explicit inverse formula. When used with a 2x2 symbolic matrix, the evaluation of the above function executes the previous code, which emits corresponding pseudo-code directives stored in a global variable called code. This variable can then be used to produce code for different targets. Here is the C code obtained after calling inv on a symbolic 2x2 matrix:

-nsp->codegen_init(); -nsp->A=symbolics(rand(2,2)); -nsp->inv(A); -nsp->global('code'); -nsp->[L,code,declarations]=code_optimize(list(),code,list(),list(),opt=[%t,%t,%f]); -nsp->txt=code_printer_c(code,list()) txt = s (10x1)

tmp_2[0]=(tmp_1[3]); tmp_2[3]=(tmp_1[0]); tmp_2[2]=(-(tmp_1[2])); tmp_2[1]=(-(tmp_1[1])); tmp_19=(((tmp_1[0])*(tmp_1[3]))-((tmp_1[2])*(tmp_1[1]))); tmp_20[0]=((tmp_2[0])/ tmp_19); tmp_20[2]=((tmp_2[2])/ tmp_19); tmp_20[1]=((tmp_2[1])/ tmp_19); tmp_20[3]=((tmp_2[3])/ tmp_19);
By looking at the code of function inv_bvar, it may seem surprizing that its execution emits pseudo-code. But, note that affectations, minus operation and multiplications are also overloaded. It is precisely the role of the overloaded functions to produce the required pseudo-code.

As an example of a function which contains explicit code directives, we give the overloaded code for the overloaded unary minus function: When the function argument is not symbolic then the function just performs ( -valueof(in)) and no code is emitted. When the variable is symbolic and of size 1x1 the function returns a new symbolics variable (out=symbolics(-valueof(in),getunique()) and the code is emitted by the call to the expression function (expression("-",list(in),overflow_option)) and by the call to the function gen_def(out,rhs). When the variable is a general matrix, the same function is applied to every entry of the matrix.

function out = minus_bvar(in) //

Code generation directives

A number of specific Nsp functions are used as directives to control code generation during the execution of the script. Some are easy to understand, for example directives to initiate and terminate the code generation process: codegen_init() textout=codegen_finalize() or to start and end a function definition in the target language:

StartFunction(f,_io) EndFunction()
More complex directives are used to specify properties of the generated code that cannot be naturally expressed in Nsp due to semantic differences between the Nsp language and the target language. The most important directives are described below.

Creating persistent variables

Persistent or static variables do not exist in Nsp but they exist in the target languages considered for code generation (in particular C) and are needed. So special directives are introduced for creating them. A pool of persistent variables may be created with the functions persistent_create() and persistent_insert(...). Each variable in the pool of persistent variables is a bvar object of type symbolics. When code is emitted, a top level declaration (of static type) will be generated for each persistent variable used in the generated code. The function persistent_extract() can be used to obtain the persistent variable.

When the function persistent_insert is used with a variable name already present in the pool of persistent variables then a local variable declaration is emited and the local variable value is copied into the persistent variable.

In the generated code a special function called initialize will contain code which resets all the persistent variables to their default values.

-nsp->// demo for creation of a pool of persistent variables -nsp->codegen_init(); -nsp->states=persistent_create(); -nsp->states=persistent_insert(states,'x1',1:3); -nsp->states=persistent_insert(states,'x2',7); -nsp->// note that not used persistent declarations are removed from the -nsp->// generated code (x3 is unsued).

-nsp->states=persistent_insert(states,'x3',1:6); -nsp->io=inouts(); -nsp->StartFunction("foo",io); -nsp->code_insert('annotation','copy [4:6]

Creating function arguments

Function arguments are created with the use of the functions inouts and inouts_insert. Arguments are bvar objects of type symbolics. When the inouts_insert function is used with a variable name already present in a sequence of function arguments, then a local variable declaration is emited and the local variable value is copied into the function argument.

In the following example a function foo with one argument is created. When calling the function foo, the persistent variable x1 is filled with the value of the function foo argument and then the function foo argument is filled with the constant value [START_REF] Weinstein | A source transformation via operator overloading method for automatic differentiation in matlab[END_REF]8,9].

-nsp->codegen_init(); -nsp->states=persistent_create(); -nsp->states=persistent_insert(states,'x1',1:3); -nsp->// in and out arguments of the function -nsp->io=inouts(); -nsp->io=inouts_insert(io,'inouts1', [START_REF] Nsp | [END_REF]8,9]); -nsp->// -nsp->StartFunction("foo",io); -nsp->code_insert('annotation','copy inouts1 to x1 with memcpy'); -nsp->states=persistent_insert(states,'x1',io.inouts1); -nsp->code_insert('annotation','copy [START_REF] Weinstein | A source transformation via operator overloading method for automatic differentiation in matlab[END_REF]8,9] into inouts1'); -nsp->io = inouts_insert(io,'inouts1', [START_REF] Weinstein | A source transformation via operator overloading method for automatic differentiation in matlab[END_REF]8,9]); -nsp->EndFunction("foo",io); -nsp->txt=codegen_finalize();

Function constant

The declaration and initialization of local variables (automatic variables) may be performed with the constant function.

-nsp->// declaration for a constant -nsp->codegen_init(); -nsp->out=constant([5,0;7,8],'x1'); // declaration + set -nsp->global('declarations','code','top_declarations','text'); -nsp->txt=code_printer_c(code,declarations) txt = s (1x1)

double x1[]={ 5, 7, 0, 8 };

Function expand

The function expand returns in a new symbolic variable. In particular, expand(in,m,n) returns a new symbolic bvar variable with the same type as the variable in but with size mxn. It also inserts a declaration for the new returned variable in the generated code.

-nsp->// demo for expand -nsp->codegen_init(); -nsp->in=numerics(%t); -nsp->var=expand(in,2,3) var = "tmp_2",%f,BMat -nsp->global('declarations','code','top_declarations','text'); -nsp->txt=code_printer_c(code,declarations) txt = s (1x1) int tmp_2[]={ TRUE, TRUE, TRUE, TRUE, TRUE, TRUE }; -nsp-> -nsp->

Functions bvarcopy and bvarempty

The functions bvarcopy and bvarempty are used to create new symbolic bvar variables with same type and size as their argument. In the case of bvarcopy, a code for copying the value of the argument in to the value of the newly created variable is also emitted.

-nsp->// demo for bvarcopy and bvarempty -nsp->codegen_init(); -nsp->in=numerics(rand(2,3)); -nsp->var=bvarcopy(in); // var is declared and set or mcopy is generated -nsp->global('declarations','code','top_declarations', ' 5 Application to code generation for Scicos

Code generation script

Code generation script is an Nsp script generated automatically from a Scicos/VSS model, and in particular from the content of the Super Block for which code generation is requested. Scicos compiler provides all the information needed to generate this script. The information includes in particular the list of blocks and their order of execution for the initialization, state update and output update phases. The compilation result contains also the data types and sizes of all the link variables inside the Super Block and those of its input and output ports.

To generate the script, the model variables are first determined. These variables contain the signals on various links and the internal block states. They lead to static C variables in the final C code because they have to retain their values from one call to the next. Links that transfer constant values are partial evaluated and do not appear in the generated code. This is done by propagating constants through the diagram, when possible.

The script starts with a section to determine the variables' initial values, which are then explicitly deefined in the generated code. This section calls the block routines with flag -1 and does not use bvar type variables. Subsequently the generated script calls in proper order the block routines with flags 1 and 2. The input and output values in these calls are of type bvar and the C code is generated thanks to operator overloading.

Unlike the code generation script, which is generated automatically and deals explicitly with variables of type bvar, the developer of the Nsp Scicos/VSS block functions writes the function codes in standard Nsp as if these functions were only used for simulation. In fact, he does not even need to know that the arguments of his functions may be at any time be of any type other than regular Nsp data types. There are however certain limitations on what he can use in his code. In particular, no conditional tests may be based on values that cannot be constant propagated (partial evaluated to a numeric value) and the compiler should be able to determine the number of iterations in for and while loops, assuming the types and sizes of all the input arguments to his functions are known. Note also that not all Nsp primitives may be used; only those that have been overloaded for the bvar type may be used (most basic Nsp primitives have already been overloaded).

The constraint on the usage of conditional statements such as the if then else construct may seem severe, but it turns out that for describing Scicos/VSS block behaviors, conditional statements based on non-constant propagated conditions are rarely encountered. In most cases such conditions may be expressed as expressional if then else and select case statements, which are implemented as If_exp and Select_exp functions.

An expressional if then else statement expressed as If_exp function is used as follows:

out = If_exp(cond, exp1, exp2)
If cond is true, out takes the value of the expression exp1, if not, it takes the value of the expression exp2. Note however that in either case both expressions exp1 and exp2 are evaluated, unlike for a real if then else statement where only one of the expressions is evaluated. These structurel conditions used in particular for subsampling in Scicos and VSS are realized using the special blocks IfThenElse and ESelect.

Conditional subsampling

We have seen so far that in order for the code generation script to function, the expressions of all conditional statements (if-then-else, when,...) must be partial evaluated, which means that it cannot be dependent on symbolic variables. But Scicos and VSS models may include conditional subsampling using in particular the IfThenElse block. Subsampling is an important mechanism available in Scicos, so a special method is introduced to handle such models for efficient code generation.

Figure 4: Model for testing the coding scheme.

An example of Scicos model with conditional subsampling is given in Fig. 4. The content of the Super Block for which C code is generated is illustrated in Fig. 5. This model implements a simple coding scheme where the Super Block receives a sequence of 0 and 1 integers and generates a sequence providing the number of consecutive 0's or 1's. The Super Block contains an IfThenElse block for testing whether or not two successive integers have the the same value, and depending on the result, a counter is either incremented or reset to 1.

Figure 5: The Super Block modeling the coding scheme using subsampling. Code is generated for this Super Block. as pointers, so that these external routines can only modify these inputs/outputs (in general only outputs) and global variables.

Example

The methodology presented in this paper may be used to generate code for models using existing Scicos blocks but it can also be used to generate code for models containing sciblks, i.e. blocks where the user has expressed the block's behavior in Nsp. The following example is adapted from an extended Kalman filter application developed in [3] for Simulink with an embedded Matlab block. Note that this is a very general Nsp code; it uses a number of matrix manipulation operators and functions. All these operators and functions are overloaded. The constraints on the usage of conditional statements and loops are also satisfied since there are no such statements in the code.

Conclusion

This paper has provided a short presentation of a new methodology based on partial evaluation and overloading used for the implementation of a code generation tool developed for Scicos and VSS. The code is freely available in Nsp 1.0 (See [START_REF] Nsp | [END_REF]).

A C code generated for the extended Kalman filter example

; static double link10024[]={ 0, 0, [0]); tmp_5=(z_10021 [START_REF] Campbell | Modeling and Simulation in Scilab/Scicos with ScicosLab 4[END_REF]); tmp_8=sqrt(((tmp_3*tmp_3)+(tmp_5*tmp_5))); tmp_11=atan2((z_10021 [START_REF] Campbell | Modeling and Simulation in Scilab/Scicos with ScicosLab 4[END_REF]),(z_10021 [0])); tmp_12 [0]=tmp_8; tmp_12 [START_REF] Adacore | GNAT Pro MOdel Compiler User Guide[END_REF]=tmp_11; /* Begin concatr of tmp_17 with unknown*/

Figure 1 :

 1 Figure 1: Main diagram.

Figure 2 :

 2 Figure 2: Content of Super Block for which C code is generated.

•

  io: a list of matrices of possibly different types and sizes representing block input and output values • state: a list of matrices of possibly different types and sizes representing block states • params: an Nsp structure coding the names and values of block parameters. The second argument of the function is a flag that indicates what job the function should perform.

Figure 3 :

 3 Figure 3: Super Block is automatically replaced with a C block for testing.

  function y=f_bvar(A) y=size(A.get_value[]);endfunction; Then calling f(A) will return the size of the value stored in the bvar variable A: -nsp->f(A) ans = r (1x2) | 3 4 |

Figure 6 :

 6 Figure 6: Model of an extended Kalman filter.

/

  * Scicos Computational function * Generated by Code_Generation toolbox of Scicos with scicos4.4.1 * date : 06 mars 2015 */ #include <scicos/scicos_block4.h> #include <string.h> #include <stdio.h> #include <stdlib.h> #include <stdint.h> #include <math.h> /* Start1002*/ void quote(double *res, double *a, double *dm,double *dn) { int i,j, m1=(int) (*dm), n1 = (int) (*dn) ; for (i = 0 ; i < (m1); i++) for (j = 0 ; j < (n1); j++) { res[j+(n1)*i]= a[i+(m1)*j]; } } void mult(double *res, double *a, double *b,double *md1,double *nd1,double *md2,double *nd2) { int i,j,k,m1=(int) (*md1),n1= (int) (*nd1),m2= (int) (*md2),n2=(int) (*nd2); for (i = 0 ; i < m1; i++) for (j = 0 ; j < n2; j++) { res[i+m1*j]=0; for (k = 0 ; k < n1; k++) res[i+(m1)*j] += a[i+(m1)*k]*b[k+(m2)*j];

  contains two functions for state update and output update:

	tmp_5[0]=tmp_1;
	tmp_5[1]=*inouts1;
	/* MUX block ends.*/
	inouts2[0]=tmp_5[0];
	inouts2[1]=tmp_5[1];
	}
	void updateState10001(double *inouts1,double *inouts2){
	z_10001=link10004;
	z_10002=*inouts1;
	}
	/* End1000*/
	void toto1000(scicos_block *block,int flag)
	{
	if (flag == 1) {
	updateOutput10001((GetRealInPortPtrs(block,1)),(GetRealOutPortPtrs(block,1)));
	}
	else if (flag == 2) {
	updateState10001((GetRealInPortPtrs(block,1)),(GetRealOutPortPtrs(block,1)));
	#include <scicos/scicos_block4.h> }
	#include <string.h> else if (flag == 4) {
	#include <stdio.h> initialize1000();
	#include <stdlib.h> }
	#include <stdint.h>
	#include <math.h>
	typedef int boolean;
	/* Start1000*/
	static double z_10001=0;
	static double z_10002=0;
	static double link10004=0;
	void initialize1000(){
	static double tmp_13=0;
	static double tmp_14=0;
	static double tmp_15=0;
	z_10001=tmp_13;
	z_10002=tmp_14;
	link10004=tmp_15;
	}
	void updateOutput10001(double *inouts1,double *inouts2){
	double tmp_1;
	double tmp_5[2];
	/* Gain block begins.*/
	/* Gain block ends.*/
	tmp_1=z_10001;
	/* Sum block begins with 2 inputs.*/
	link10004=(tmp_1-z_10002);
	/* MUX block begins with 2 inputs.*/

VisSim SIMULATE is a commercial modeling and simulation product based on Scicos.

* This code generation tool was developed during a three-years research projects funded within the French FUI 2011

The generated code shows that the test implemented by the IfThenElse block has produced an if statement in C: The code geneation script uses the if_cos function, which uses the If_exp function. The if_cos function is defined as follows:

function if_cos(in,f1,f2) // f1 and f2 are expressions of type call(..) if is_sym(in) then code_insert("if_expr", in > (0),f1,f2); elseif valueof(in)>0 then code_insert("ident",f1) else code_insert("ident",f2) end endfunction Note that the if_cos function does not return any value. This is consistent with the fact that the IfThenElse block has no regular output. The second and third arguments of this function are calls to external (generated) C routines. The arguments are the Super Block inputs and outputs, passed