S Zozor

P Ravier

O Buttelli

On Lempel-Ziv complexity for multidimensional data analysis

Keywords: Complexity measures, Lempel-Ziv complexity, Shannon Entropy, nonlinear deterministic multidimensional systems, random boolean network, minority game PACS: 87.10.+e, 02.50.Le, 64.60.Cn, 02.50.Sk

In this paper, a natural extension of the Lempel-Ziv complexity for several finitetime sequences, defined on finite size alphabets is proposed. Some results on the defined joint Lempel-Ziv complexity are given, as well as properties in connection with the Lempel-Ziv complexity of the individual sequences. Also, some links with Shannon entropies are exhibited and, by analogy, some derived quantities are proposed. Lastly, the potential use of the extended complexities for data analysis is illustrated on random boolean networks and on a proposed multidimensional extension of the minority game.

Introduction

Many physical, biological or financial signals result from the dynamics of great dimensional systems. As an example, in biology, the cells of the cardiac tissue exchange ions with the extracellular field by a nonlinear reaction, and with the connected cells by (linear) diffusion. Hence, electric waves propagate on the tissue, and the electrocardiogram (ECG) is the electric field produced by the propagation and measured by an electrode. Another example can be found in the collective actions of the genes, generating the production of proteins in certain quantities. In finance, one can take the example of the variation in the price of an asset, resulting from the collective action of the buyers and sellers. In all these examples, the challenge is to describe or extract meaningful information from the data [START_REF] Rajković | Extracting meaningful information from financial data[END_REF][START_REF] Rajković | Quantifying complexity in the minority game[END_REF], to characterize, to detect or to classify different pathologies [START_REF] Botteron | A technique for measurment of the extent of spatial organization of atrial activation during atrial fibrillation in the intact human heart[END_REF][START_REF] Cysarz | Entropies of short binary sequences in heart period dynamics[END_REF][START_REF] Owis | Study of features based on nonlinear dynamical modeling in ECG arrhythmia detection and classification[END_REF][START_REF] Quiroga | Kulback-Leibler and renormalized entropies: Applications to electroencephalograms of epilepsy patients[END_REF][START_REF]Task Force of the European Society of Cardiology and the North American Society of Pacing and Electrophysiology, Heart rate variability: Standards of measurement, physiological interpretation, and clinical use[END_REF][START_REF] Thakor | Ventricular tachycardia and fibrillation detection by a sequantial hypothesis testing algorithm[END_REF][START_REF] Wang | A short-time multifractal approach for arrhythmia detection based on fuzzy neural network[END_REF][START_REF] Zhang | Detecting ventricular tachycardia and fibrillation by complexity measure[END_REF], etc.

Due to the complex origin of the signals, researchers generally use tools from either information theory or from dynamic nonlinear processes. The goal is to characterize as well as possible the degree of organization of the measured signals. The first approach is statistical, and the tools generally employed are entropies [START_REF] Cysarz | Entropies of short binary sequences in heart period dynamics[END_REF][START_REF] Quiroga | Kulback-Leibler and renormalized entropies: Applications to electroencephalograms of epilepsy patients[END_REF][START_REF] Ebeling | Entropy, complexity, predictability and data analysis of time series and letter sequences[END_REF][START_REF] Pincus | Approximate entropy (ApEn) as a complexity measure[END_REF][START_REF] Torres | Relative complexity changes in time series using information measures[END_REF] (or "approximate entropies" or multi-resolution entropies), correlation measures [START_REF] Botteron | A technique for measurment of the extent of spatial organization of atrial activation during atrial fibrillation in the intact human heart[END_REF], spectral analysis [14], etc. The second approach is based on the fact that the underlying mechanisms producing the measured signals are generally deterministic and nonlinear. This is clear in the example of ECG, even if measurements may be corrupted by noise. The nonlinear tools generally employed to study the signals come from chaos theory such as Lyapunov exponents or dimensions (fractal, etc.) [START_REF] Rajković | Extracting meaningful information from financial data[END_REF][START_REF] Owis | Study of features based on nonlinear dynamical modeling in ECG arrhythmia detection and classification[END_REF][START_REF] Wang | A short-time multifractal approach for arrhythmia detection based on fuzzy neural network[END_REF]. Other tools often employed come from the concept of complexity in the sense of Kolmogorov, more particularly Lempel-Ziv complexity [START_REF] Rajković | Quantifying complexity in the minority game[END_REF][START_REF] Thakor | Ventricular tachycardia and fibrillation detection by a sequantial hypothesis testing algorithm[END_REF][START_REF] Zhang | Detecting ventricular tachycardia and fibrillation by complexity measure[END_REF][START_REF] Radhakrishnan | Estimating regularity in epileptic seizure time-series data -a complexity-measure approach[END_REF][START_REF] Radhakrishnan | Quantifying physiological data with Lempel -Ziv complexity -certain issues[END_REF][START_REF] Szczepański | Application of Lempel -Ziv complexity to the analysis of neural discharges[END_REF][START_REF] Zhang | EEG complexity as a measure of depth of anesthesia for patients[END_REF]. Several complexity measures for data analysis have already been proposed as presented in [START_REF] Ebeling | Entropy, complexity, predictability and data analysis of time series and letter sequences[END_REF] and in the references therein. However, in spite of the terminology, the approach is generally more likely statistical. In this paper, we will focus on Lempel-Ziv complexity. The motivation for this is that the data generally studied in biomedical engineering have a nonlinear and deterministic origin. Furthermore, tools such as Lyapunov exponents are difficult to evaluate 1 and require long time-computation. In contrast, the Lempel-Ziv complexity contains the notion of complexity in the deterministic sense (Kolmogorov sense) as well as in a statistical sense (Shannon sense), and can be computed with a low cost. In the literature, this tool is generally employed to analyze mono-dimensional signals. In this paper we show that contrary to what has been previously claimed [START_REF] Thakor | Ventricular tachycardia and fibrillation detection by a sequantial hypothesis testing algorithm[END_REF], this tool can be used for multidimensional signals. The goal is to use the natural extension of the Lempel-Ziv complexity to characterize the "regularity" of a system using several signals.

Section 2 recalls the basics on Lempel-Ziv complexity (LZC). Section 3 then shows how the LZC can naturally be extended for multidimensional sequences.

In this section we will give some properties of multidimensional LZC in connection with "marginal" LZCs. Some parallels and differences with Shannon entropy will be exhibited. In section 4 we will illustrate how the extended LZC can capture the spatiotemporal organization of the signal delivered by a multidimensional nonlinear and nonlinearly coupled system. The illustration will be done for two particular binary systems: random boolean networks and an extension of the minority game. Finally, section 5 will give some conclusions and perspectives.

Basics on the Lempel-Ziv complexity

Consider a sequence (or a word) S = s 1 . . . s n of length n where each letter s i is taken from an alphabet A of finite size α. Kolmogorov in 1965 defined the complexity of such a sequence as the size of the smallest binary program that can produce this sequence [START_REF] Cover | Elements of Information Theory[END_REF]. The complexity in the sense of Kolmogorov seems very general and computer-dependent. It was shown that even if a Kolmogorov complexity can be defined up to a constant, the evaluation of such a complexity in a finite time is not guaranteed [START_REF] Cover | Elements of Information Theory[END_REF]. Several years later, in their seminal work [START_REF] Lempel | On the complexity of finite sequences[END_REF], A. Lempel and J. Ziv proposed to define a complexity in the sense of Kolmogorov, but limiting their definition to programs based on two operations: recursive copy and paste operations. They defined two fundamental notions to evaluate the Lempel-Ziv complexity (LZC): production and reproduction operations. As in the original paper, we will denote SQ the concatenation of two words S and Q, (S) the length of a sequence S, S(i, j) the sub-sequence s i s i+1 . . . s j of S and π the operation suppressing the last letter of a sequence, i.e. Sπ = S(1, (S) -1).

Reproduction: An extension R = SQ of a sequence S is said reproducible from S if the sequence Q is in the vocabulary of SQπ, i.e. if Q is a sub-sequence of SQπ. As an example, for S = 101 and Q = 010, R = SQ = 101010 is a reproducible extension of S since q 1 = s 2 = r 2 , q 2 = s 3 = r 3 and q 3 = q 1 = r 4 previously copied, and Q is a sub-sequence of SQπ = 10101. In other words, SQ can be reproduced from S by recursive copy and paste operations. As in the initial paper, we will denote the reproduction by S -→ R and the index p ≤ (S) such that Q = R(p, (Q) + p -1) is called the pointer for the reproduction S -→ R.

Production: A non-empty sequence S is said to be producible from a prefix S(1, j) if S(1, j) -→ Sπ and j < (S). From the previous example, it can be seen that S = 101 produces SQ = 101010 since S can reproduce SQπ = 10101. But S also produces T = 101011 while S does not reproduce this sequence. The difference between reproduction and production is that in production the last letter can come from a supplementary copy-paste but can also be "new".

We will denote production S(1, j) =⇒ S and S(1, j) is called a basis of S.

To understand how a program can build a sequence using these two operations, consider a given word S. An index h 2 can be found such that s 1 =⇒ S(1, h 2 ). h 2 = 2 always works, and if s 2 = s 1 one can also choose h 2 = 3. If s 2 = s 1 and s 3 = s 2 , then h 2 = 4 also works; and so on. Then, there is an index h 3 so that S(1, h 2 ) =⇒ S(1, h 3 ). And so on. Thus, S can be built as [START_REF] Lempel | On the complexity of finite sequences[END_REF]. The process Hi(S) = S(1, h 1 )S(h 1 + 1, h 2 ) . . . S(h m-1 + 1, h m ) is called the history of the production and the sub-sequences S(h i-1 + 1, h i ) are called the components of the process. The size of the process is then defined as the number of components of the process, i.e. c Hi (S) = h m . Following the idea of Kolmogorov, Lempel and Ziv sought the shortest production possible process. Hence they defined a complexity in the sense of Lempel-Ziv as

s 1 = S(1, h 1 ) =⇒ S(1, h 2 ) =⇒ S(1, h 3 ) • • • =⇒ S(1, h m ) = S
c(S) = min Hi∈{histories of S} c Hi (S) ( 1 )
This quantity is not precisely a complexity since it does not give directly the program and its size that can produce S, but c(S) is directly linked to the size of such a program. As an illustration of LZC, for the sequence S = 0100100100 the process 0 =⇒ 01 =⇒ 0100 =⇒ 0100100 =⇒ 010010010 =⇒ 0100100100 is of size 6 while the minimal process 0 =⇒ 01 =⇒ 0100 =⇒ 0100100100 is only of size 4. In this example, it can be seen that in the minimal process each S(1, h i+1 ), except the last one, cannot be reproduced but only produced by S(1, h i ). In [START_REF] Lempel | On the complexity of finite sequences[END_REF], Lempel and Ziv defined as exhaustive components S(h i-1 + 1, h i ) all the components that can be produced but not reproduced by S(1, h i-1 ) and a history of a process is said to be exhaustive if all the components are exhaustive, except possibly the last one. Then, they showed that an exhaustive history is unique and that the number of components of the exhaustive history is precisely the LZC [START_REF] Lempel | On the complexity of finite sequences[END_REF]. This last remark led to algorithms to evaluate the LZC of a sequence (see e.g. [START_REF] Kaspar | Easily calculable measure for the complexity of spatiotemporal patterns[END_REF]). With the procedure to evaluate the LZC, it can be seen that if the length of each component, each pointer and the last letter of each component are stored, the complete sequence can easily be retrieved2 .

In conclusion, the LZC contains the notion of complexity in the sense of Kolmogorov. Furthermore, if the sequence S (n) = s 1 . . . s n is random and if the source is stationary and ergodic it can be shown that lim n→+∞ c(S (n) ) log(n) n = H(S) where H(S) is the entropy rate of the source [START_REF] Cover | Elements of Information Theory[END_REF][START_REF] Lempel | On the complexity of finite sequences[END_REF]: this result says that the LZC also contains a notion of "average information quantity" in the Shannon sense.

Since the LZC tries to capture a degree of redundancy, or patterns that are similar in a sequence, this tool seems interesting for the analysis of sequences that appear complex but that may hide some simple underlying behaviors. In this way, the LZC has been proposed to analyze chaotic sequences [START_REF] Kaspar | Easily calculable measure for the complexity of spatiotemporal patterns[END_REF]. Furthermore, since the LZC can also be viewed as an estimation of Shannon entropy, it seems that it can also be used for non deterministic signals. This tool seems to bridge the two above-mentioned approaches. However, although the LZC has already been employed for data analysis [START_REF] Thakor | Ventricular tachycardia and fibrillation detection by a sequantial hypothesis testing algorithm[END_REF][START_REF] Zhang | Detecting ventricular tachycardia and fibrillation by complexity measure[END_REF][START_REF] Radhakrishnan | Estimating regularity in epileptic seizure time-series data -a complexity-measure approach[END_REF][START_REF] Radhakrishnan | Quantifying physiological data with Lempel -Ziv complexity -certain issues[END_REF][START_REF] Szczepański | Application of Lempel -Ziv complexity to the analysis of neural discharges[END_REF][START_REF] Zhang | EEG complexity as a measure of depth of anesthesia for patients[END_REF], it was only used for scalar sequences. In the next section, we will show that the LZC is in fact also naturally defined for vectorial sequences and thus can be applied for multidimensional data analysis. Furthermore we will draw a parallel with Shannon entropy, pointing out some similarities, but also some differences.

Multidimensional Lempel-Ziv complexity

As far as we know, the first attempt to use the LZC for spatiotemporal data analysis was made by Kaspar and Schuster [START_REF] Kaspar | Easily calculable measure for the complexity of spatiotemporal patterns[END_REF]. In their paper, they proposed to analyze spatiotemporal signals defined on a finite size alphabet by computing a spatial LZC at each time t, i.e. for N sequences x i (t), i = 0, . . . , N -1, c(t) = c(x 0 (t) . . . x N -1 (t)). Then, with their approach, the time evolution of the spatial complexity c(t) can be analyzed to deduce a spatiotemporal behavior. As an example, if the spatial LZC decreases with time, it can be concluded that the signal/system tends to be more and more organized. However, such an approach remains only spatial, and the temporal links are omitted. To better understand this fact, if the signal has only two components (e.g. from two sensors), by definition the spatial LZC remains constant and equal to 2 at each time t. Hence, the evolution of the spatial LZC fails in describing the spatiotemporal behavior of the signal: in this example the spatial dimension is not large enough.

A more natural approach is to extend directly the LZC for vectorial data. This can be done naturally by extending the alphabet. To this aim, consider k sequences X i = x i,1 . . . x i,n for i = 0, . . . , k -1, where the letters are respectively in the alphabets A 0 , . . . , A k-1 of respective sizes α 0 , . . . , α k-1 . Consider now the sequence Z = z 1 . . . z n , defined on the extended alphabet

B = A 0 × • • • × A k-1 of size α 0 . . . α k-1
, where the components are the kuplets z j = (x 0,j , . . . , x k-1,j ). Notice that Z is a sequence of n k-uplets and not a sequence of k × n letters: sequence Z does not result from a letter mixing approach. Since the equality relation holds for k-uplets, the production and reproduction operations defined by Lempel and Ziv hold for k-uplets. As a conclusion, all the work of Lempel and Ziv remains valid for vectorial sequences, although not explicitly spelled out in their paper [START_REF] Lempel | On the complexity of finite sequences[END_REF]. But once again, it is important to notice that for vectorial sequences the alphabet has no scalar elements but only k-uplet elements. Hence, we can define what we will call the joint Lempel-Ziv complexity of sequences X 0 , . . . , X k-1 by

c(X 0 , . . . , X k-1 ) = c(Z) ( 2 )
Furthermore, if the alphabets are the same and are of the form A = {0, . . . , α-1}, we can also define a sequence Z = z 1 . . . z n considering that each z j has the x i,j as α-ary decomposition, i.e. z j = k-1 i=0 x i,j α i . Defining joint LZC of the X i as that of Z is exactly similar to the previous definition. The last assumption is not restrictive since bijections can be found from the A k to a subset of A of size α k to achieve such a case, provided that α ≥ max k α k . Moreover, using this approach, the algorithm proposed in [START_REF] Kaspar | Easily calculable measure for the complexity of spatiotemporal patterns[END_REF] to evaluate the LZC can still be used, comparing scalars. With this definition, the LZC of multidimensional signals can then be viewed as a joint LZC and analyzed regarding the LZC of its components. The first property of the joint LZC is its obvious symmetry by permutation of the X j :

Property 1 The joint LZC is invariant by any permutation σ of {0, . . . , k - 1}, i.e. c(X 0 , . . . , X k-1 ) = c(X σ(0) , . . . , X σ(k-1) ) ( 3 )
PROOF. Consider the exhaustive history of the vectorial sequence Z = z 1 . . . z n where z j = (x 0,j , . . . ,

x k-1,j ), i.e. E Z (Z) = Z(1, h 1 ) . . . Z(h m-1 + 1, h m ). Hence, whatever i = 2, . . . , m, we have both Z(1, h i-1 ) =⇒ Z(1, h i ) and Z(1, h i-1 ) -→Z(1, h i ).
Thus, for any i, let us denote p i the pointer

p i ≤ h i-1 such that both Z(h i-1 + 1, h i -1) = Z(p i , p i + h i -h i-1 -2) and z h i = z p i +h i -h i-1 -1 . Furthermore, for any q i ≤ h i-1 , Z(h i-1 + 1, h i ) = Z(q i , q i + h i - h i-1 -1). The direct consequence is for all j, X j (h i-1 + 1, h i -1) = X j (p i , p i + h i -h i-1 -2) and (x 0,h i , . . . , x k-1,h i ) = (x 0,p i +h i -h i-1 -1 , . . . , x k-1,p i +h i -h i-1 -1 ). Furthermore, whatever q i ≤ h i-1 there is at least one j q i ∈ {0, . . . , k -1} so that X jq i (h i-1 + 1, h i ) = X jq i (q i , q i + h i -h i-1 -1)
. The equality between the subsequences and the inequality between the k-uplets remain unchanged by any permutation σ in {0, . . . , k -1} and the inequality between the sequences holds with σ(j q i ), that finishes the proof. 2

Due to this property, the joint LZC is uniquely defined. Note that a definition using a sample mixing approach would destroy this symmetry (e.g. X 0 = 00001 and

X 1 = 10100 leads to c(X 0 , X 1 ) = c(X 1 , X 0 ) = 4 = c(x 0,1 x 1,1 . . . x 0,5 x 1,5 ) = c(0100010010) = 5 = c(x 1,1 x 0,1 . . . x 1,5 x 0,5 ) = c(1000100001) = 4
). Furthermore, by a sample mixing approach, one can feel that the temporal links of a sequence will be less well captured than by using this natural approach.

Since the joint LZC is simply an LZC as defined by Lempel and Ziv, it possesses all the properties of the LZC given in [START_REF] Lempel | On the complexity of finite sequences[END_REF]. The first of them is the asymptotic link with Shannon entropy:

Property 2 Consider k random sequences X 0 , . . . , X k-1 that are jointly stationary and ergodic. Hence

lim n→+∞ c(X 0 , . . . , X k-1 ) log(n) n = H(X 0 , . . . , X k-1 ) ( 4 )
where H(X 0 , . . . , X k-1 ) is the joint entropy rate of the sequences.

PROOF. Consider the sequence Z = z 1 . . . z n where z j = (x 0,j , . . . , x k-1,j ) Since the X i are jointly stationary and ergodic, Z is stationary and ergodic. On the extended alphabet, from [START_REF] Lempel | On the complexity of finite sequences[END_REF],

lim n c(Z) log(n) n = H(Z) = lim n E[-log(Pr[z 1 ,...,zn])] n . Since Pr[z 1 , . . . , z n ] = Pr[x 0,1 , . . . , x k-1,1 , . . . , x 0,n , . . . , x k-1,n ] we obtain that H(Z) = H(X 0 , . . . , X k-1 ) that finishes the proof. 2
This complexity clearly exhibits the fact that both spatial and temporal relations are captured by the joint LZC, at least statistically and asymptotically.

Notice now that in the LZC sense, a multidimensional signal is obviously more complex than each single component. This can be summarized as follows.

Property 3

The joint LZC of k sequences X 0 , . . . , X k-1 is greater than or equal to the LZC of each component, that reads c(X 0 , . . . , X k-1 ) ≥ max(c(X 0 ), . . . , c(X k-1 ))

(5)

PROOF. Consider sequence Z of the k-uplets (x 0,j , . . . , x k-1,j ) and its exhaustive history

E Z (Z) = Z(1, h 1 ) . . . Z(h m-1 + 1, h m ). Since for all i we have Z(1, h i-1 ) =⇒ Z(1, h i ), obviously for all j, X j (1, h i-1 ) =⇒ X j (1, h i ). As a consequence E Z (X j ) is also a production process of X j . Using the definition of the LZC, c(X j ) ≤ c E Z (X j ) = c(Z). 2
This property, as well as the previous one, are shared with the Shannon entropies.

Notice that Z(1, h i-1 ) -→Z(1, h i ) does not imply the same property for all X j ; this property holds at least for one of the X j . Hence, the exhaustive history of Z is not necessarily exhaustive for the X j and it can even be not exhaustive for any X j . The consequence is that the equality says nothing more in terms of links between the joint sequence and the marginal one. The following property give some links in particular cases.

Property 4 Consider X and Y two sequences defined in alphabets A x and A y respectively. If there is a bijection σ from A x to A y such that

y i = σ(x i ) for all i = 1, . . . , n, then c(X, Y ) = c(X) = c(Y ) ( 6 ) A particular case is c(X, X) = c(X) ( 7 )
The reciprocal property is false.

PROOF. Consider first the exhaustive history of X, E X (X) = X(1, h 1 ) . . . X(h m-1 + 1, h m ). Since X(1, h i-1 ) =⇒ X(1, h i ) and X(1, h i-1 ) -→X(1, h i ), there is a pointer p i ≤ h i-1 such that X(h i-1 + 1, h i -1) = X(p i , p i + h i -h i-1 -2), that x h i = x p i +h i -h i-1 -1 and whatever q i ≤ h i-1 , X(h i-1 + 1, h i ) = X(q i , q i + h i -h i-1 -1). Clearly (σ(x h i-1 +1 ), . . . , σ(x h i -1 )) = (σ(x p i ), . . . , σ(x p i +h i -h i-1 -2 )
) and the bijectivity of function σ yields that σ(

x h i ) = σ(x p i +h i -h i-1 -1 )
and that whatever

q i ≤ h i-1 , (σ(x h i-1 +1 ), . . . , σ(x h i )) = (σ(x q i ), . . . , σ(x q i +h i -h i-1 -1 )
). As a conclusion, E X (Y ) is an exhaustive history of Y and then clearly c(X) = c(Y ). Now, consider Z of components z j = (x j , y j ). We immediately obtain that

Z(h i-1 + 1, h i -1) = Z(p i , p i + h i -h i-1 -2), that z h i = z p i +h i -h i-1 -1 and that whatever q i ≤ h i-1 , Z(h i-1 + 1, h i ) = Z(q i , q i + h i -h i-1 -2)
. As a consequence, E X (Z) is also an exhaustive history of Z, that finishes the proof. Concerning non reciprocity, consider the example of two binary strings, X = 010010 and Y = 011011. We obtain c(X) = c(Y ) = c(X, Y ) but Y cannot be predicted through X with a bijection from {0, 1} to itself. 2

This property can naturally be extended for more than two sequences. Furthermore, this property is similar to that of Shannon entropy. Indeed, for two random sources X and Y , if and only if Pr[Y | X ] = 1 (i.e. Y completely determined by X), then H(X, Y ) = H(X). But contrary to Shannon entropy, c(X) = c(X, Y ) does not imply that X is completely determined by Y . Because of the non reciprocity of this property for the LZC, some care will have to be taken in the results of the joint LZC analysis, when an analogy is made with Shannon entropy.

To go farther with the analogy between Shannon entropy and LZC, one can define conditional LZC (CLZC) and informational LZC (ILZC) by

       c(Y |X) = c(X, Y ) -c(X) I c (X; Y ) = c(X) + c(Y ) -c(X, Y ) = c(X) -c(X|Y ) = c(Y ) -c(Y |X) ( 8
) and so on for higher order informational LZC (as for the Shannon entropies [START_REF] Matsuda | Physical nature of higher-order mutual information: Intrinsic correlations and frustration[END_REF]). Since the joint LZC and the marginal LZC tend to the joint and marginal Shannon entropies for stationary ergodic sequences, this holds for the conditional LZC and for informational LZC, respectively toward conditional entropy and mutual information. Furthermore, property 3 leads to the positivity of the conditional complexity as for the conditional Shannon entropy.

In many classification problems or source separation problems, the Kullback-Leibler divergence, and particularly mutual information, is used as a separation criterion. In a certain sense, mutual information can be understood as a distance between the joint density of two sequences and independent sequences sharing the same marginal densities. Using the LZC, by analogy with the Shannon entropy, I c can be expected to be understood as a divergence between two sequences. However, the ILZC I c can be negative, e.g. X = 00100100100 and Y = 01010100000 leads to I c (X, Y ) = -1. This nonproperty of non-negativity is another difference with mutual information in the Shannon sense. Furthermore, the non-reciprocal of property 4 and the non existence of triangular inequality does not permit to build a metric with ρ = c(X|Y ) + c(Y |X) as can be done, in a certain sense, with Shannon entropy (see [START_REF] Cover | Elements of Information Theory[END_REF]), e.g. X = 00110, Y = 01000 and Z = 01010 leads to

4 = ρ(X, Z) > ρ(X, Y ) + ρ(Y, Z) = 1 + 1 = 2.
The LZC tool seems interesting in itself for the analysis of complicated multidimensional data (biomedical, economic,. . . ), since it makes a kind of bridge between nonlinear tools for deterministic data analysis and information theory tools. However, this bridge is not complete because of the limited analogies with Shannon entropy. As a conclusion, even if this tool is interesting in itself to characterize or to classify signals, the results obtained need to be interpreted with great caution .

Data analysis using the joint LZC

In this section, the use of the joint LZC and the informational LZC (ILZC) for data analysis will be illustrated on two examples. The first example illustrates how the joint LZC can capture a spatiotemporal degree of organization of a random boolean network (RBN), while the second example will show the interest of using ILZC on a proposed multidimensional variant of the minority game (MG). The RBN as well as the MG produce boolean signals, and can be naturally studied by LZCs.

Illustration on a random boolean network

An RBN, also called Kauffman network, is given by N binary automata, where each element or cell i is spatially connected to K other cells; K is called the connectivity of the network [START_REF] Ballesteros | Random boolean networks response to external periodic signals[END_REF]. Then, the temporal evolution is given by maps f i : {0, 1} K → {0, 1} such that x i (t + 1) = f i (x i 1 (t), . . . , x i K (t)). The connections of the cells are directional, i.e. one of the i l can be equal to j even if j m = i whatever m. The initial state (x 0 (1), . . . , x N -1 (1)) is randomly chosen in {0, 1} N where 1 (resp. 0) is drawn with probability p (resp. 1p). The maps f i are also randomly drawn from all the 2 2 K possible functions {0, 1} K × {0, 1} with the same probability p to draw 1. In other words, the value of f i (0, . . . , 0) is drawn in {0, 1}, . . . , the value of f i (1, . . . , 1) is drawn in {0, 1}, for all the i. This is done to initialize the process, and once built, these functions remain unchanged during the time-evolution of the process. Hence, for a given choice of map, a specific network is considered. The probability p is called the bias of the network [START_REF] Ballesteros | Random boolean networks response to external periodic signals[END_REF]. It is shown in [START_REF] Ballesteros | Random boolean networks response to external periodic signals[END_REF][START_REF] Solé | Phase transitions and antichaos in generalized kauffman networks[END_REF] that, according to the values of (p, K), the behavior of the network can produce order or disorder. Furthermore, it is shown that the network can be controlled by freezing at each time t, F (t)γN cells, where γ is the maximum proportion of cells that can be frozen and where F (t) ∈ [0 ; 1] is a control function. The cells that can be frozen are randomly chosen, but this choice is time-independent. Here again, at the beginning of the process, the cells that will be potentially frozen are randomly chosen and this choice remains unchanged during the evolution of the process. The γNF (t) cells that are frozen are the first γNF (t) cells that were drawn. Furthermore, the frozen states are randomly chosen in the initialization of the process and remain the same during all the process. As a conclusion, even if there are random choices in the network, all the random variables are drawn in the initialization, thus the spatiotemporal evolution is clearly deterministic. This RBN has been proposed as a description of discrete genetic network models. In this model, the two possible states for the gene represent if the transcription process is active (e.g. state 1) or not (e.g. state 0). Hence, the control has been proposed for the modeling of external forcing (biological rhythms, etc.). More information about the RBN can be found in [START_REF] Ballesteros | Random boolean networks response to external periodic signals[END_REF][START_REF] Solé | Phase transitions and antichaos in generalized kauffman networks[END_REF] and in the references therein.

First of all, in [START_REF] Ballesteros | Random boolean networks response to external periodic signals[END_REF][START_REF] Solé | Phase transitions and antichaos in generalized kauffman networks[END_REF] it is shown that according to the values of K and p, the network can exhibit an ordered behavior, or a more complex one, akin to a sort of spatiotemporal chaos. More rigorously, since the number of states for the cells is finite and equals 2 N , the behavior must be periodic with a period less or equal to 2 N . The chaos-like behavior is possible, analyzing the RBN during a time sufficiently small compared to 2 N . In particular, it is shown that when K < K c = 1 2p(1-p) the system exhibits order and when K > K c it exhibits chaos-like behavior. Here, we have chosen to analyze such a coupled oscillator, in order to use the joint LZC. Thus, we have built an RBN, for different values of p and K, and for different realizations of maps (i.e. several specific networks), initial conditions, etc. For each vectorial signal (i.e. one realization for the set (p, K)), the joint LZC of the N oscillators is computed and is denoted c(p, K). Figure 1 represents, then, the average LZC c(p, K) over the realizations. This figure clearly shows that under the critical curve K = 1 2p(1-p) , the mean joint LZC is low, corresponding to a quite spatiotemporally regular signal, while the mean joint LZC is high over the critical curve. This result exhibits that the joint LZC is able to capture the degree of organization of such a coupled system, at least on average. Furthermore, since the signal cannot be rigorously chaotic, an approach using a Lyapunov exponent seems not rigorous for the analysis of an RBN, even if it is possible by analogy with continuous systems as shown in [START_REF] Ballesteros | Random boolean networks response to external periodic signals[END_REF][START_REF] Luque | Lyapunov exponents in random boolean networks[END_REF]. As is concluded for chaotic signals in [START_REF] Kaspar | Easily calculable measure for the complexity of spatiotemporal patterns[END_REF], the LZC can capture similar information to a Lyapunov exponent, but its evaluation has a lower cost than that of a Lyapunov exponent. Here, we illustrate that this holds for the joint LZC previously presented.

The second illustration concerns a controlled RBN. An RBN was generated, with a free evolution in the beginning, then the RBN is controlled during a finite time and with a periodic control function, and there is again a free evolution. To analyze the spatiotemporal complexity of this system, we have first evaluated at each time a spatial LZC as proposed by Kaspar and Schuster [START_REF] Kaspar | Easily calculable measure for the complexity of spatiotemporal patterns[END_REF] for spatiotemporal analysis. This spatial LZC is compared with the LZC evaluated on a sliding window of N w samples (sliding step of 1 sample) for the first and the middle cells, i.e. at each time t the LZCs c(x 1 (t -N w + 1) . . . x 1 (t)) and c(x (N +1)/2 (t -N w + 1) . . . x (N +1)/2 (t)) (e.g. for N odd) are computed. Lastly the joint LZC on the sliding window was evaluated, i.e. at each time t the joint LZC c(x(t -N w + 1) . . . x(t)) is computed, where x(t) is the N-uplet (x 1 (t), . . . , x N (t))). The results are depicted in figure 2. In these figures it can be seen that, for this example, the spatial LZC does not capture the spatiotemporally organized area in the RBN. One can also see that the 1-dimensional LZC is sometimes able to capture the more organized area, when the periodic behavior concerns the chosen cell for the analysis. However, the analyzed cell can also have a periodic behavior during the noncontrolled period, that does not reflect the spatiotemporally non-organized area. Lastly, the joint LZC is clearly able to detect the organized period and to differentiate it from the non-organized areas. In particular, the behavior of the joint LZC reflects the controlled function. Since this joint LZC takes simultaneously into account all the cells, during the non-controlled period, its value is large, reflecting the globally non-organized area. Figure 2 also depicts the averaged spatial LZC and the averaged joint LZC (sliding window) over 500 realizations of the controlled RBN. It can be seen in these figures that on average while the spatial LZC does not capture the controlled area, the joint LZC does. This is due to the fact that for each realization, the frozen cells are not the same and the spatial links vary from one realization to the other. Since the RBN is only spatially analyzed, the temporal periodicities cannot be captured by the spatial LZC. Some regularities can be captured only when the frozen cells present regularities in their spatial choices, but on average there are not regularities. Hence on average, the spatial LZC cannot exhibit low values when the RBN is more organized. In contrast, since the joint LZC takes into account simultaneously all the cells, on average it can clearly distinguish the organized period: in each realization, even if there is not necessarily spatial periodic behavior, the N-uplet has a periodically organized behavior in time. The mean marginal LZCs are not pictured here, but they exhibit the same kind of behavior as the joint LZC. Indeed, since during the control a given cell can exhibit periodic behavior for some realizations and not for others, by averaging over the realizations the non time-periodic LZC (or periodic but not linked to that of the control) are smoothed, resulting in a periodic mean marginal LZC. During the free evolution, the marginal LZC is also smoothed by averaging over the realizations. However, since for the joint LZC the RBN is globally taken into account, during the control the regularity of the RBN is better captured, the effect is more pronounced for the average joint LZC than for each marginal LZC.

Illustration on the minority game

The minority game consists of N agents where each agent i makes a binary decision a i (t) ∈ {0, 1} at time t. The goal for each agent is to be in the minority, i.e. agent i is in the minority if the number of agents taking the same decision a i (t) is less than N 2 [START_REF] Challet | Emergence of cooperation and organization in an evolutionary game[END_REF][START_REF] Zhang | Modeling market mechanism with evolutionary games[END_REF][START_REF] Savit | Adaptive competition, market efficiency, and phase transition[END_REF]. Let us denote by µ(t) the decision of the minority of agents at time t. To make a decision, each agent looks at the m previous minorities (µ(t-m), . . . , µ(t-1)) where m is the memory of the game. Each agent possesses a look up table of s functions f k : {0, 1} m -→ {0, 1} called strategies. There are 2 2 m possible strategies and each agent has the same number s of strategies f i k , i = 0, . . . , N -1, k = 0, . . . , s -1. The i k are generally different from one agent to another, but some strategies can possibly be shared by several agents. A gain G i k is attached to each strategy f i k . At each time, the gains of all the strategies that would have led to being in the minority are increased by 1, i.e. at time t, if f i k (µ(tm), . . . , µ(t -1)) = µ(t), G i k is increased by 1 point. Hence, at each time, to make its decision, agent i looks at all the strategies at his disposal and chooses the one which has the greatest gain, i.e. i l = Argmax i k G i k . The decision of agent i is then a i (t) = f i l (µ(tm), . . . , µ(t -1)). If for agent i several strategies shared the same maximal gain over the G i k , k = 0, . . . , s -1, the chosen strategy is randomly (uniformly) chosen over these optimal strategies. To initiate the game, one only has to randomly choose the s strategies of each agent, to randomly choose an action a i (0) for each agent, and an initial memory string (µ(-m), . . . , µ(-1)). This game is often used for social, biological, ecological or financial modeling, etc (see [START_REF] Rajković | Quantifying complexity in the minority game[END_REF][START_REF] Challet | Emergence of cooperation and organization in an evolutionary game[END_REF][START_REF] Zhang | Modeling market mechanism with evolutionary games[END_REF][START_REF] Savit | Adaptive competition, market efficiency, and phase transition[END_REF][START_REF] Marsili | Continuum time limit and stationary states of the minority game[END_REF][START_REF] Challet | Shedding light on el farol[END_REF] and ref. therein). The game can be easily understood in the case of financial modeling. Say that decision 0 corresponds to buying an asset, and 1 to selling an asset. It is clear that a player prefers to sell an asset while the majority buys this asset, and vice-versa. Several variants have been proposed. For example, a strategy may be rewarded by taking into account the "quality" of the minorities. It was shown that the behavior of such a game can be characterized by the number L 1 (t) of agents choosing 1 at each time. In particular the (time-averaged) variance σ 2 1 of L 1 , known as volatility in finance, has been shown to be a good characteristic of the behavior of the MG [START_REF] Challet | Emergence of cooperation and organization in an evolutionary game[END_REF][START_REF] Savit | Adaptive competition, market efficiency, and phase transition[END_REF][START_REF] Marsili | Continuum time limit and stationary states of the minority game[END_REF]. It characterizes the total waste of the game. This variance decreases as m increases until a critical size m c of the memory m, and increases to attain the same variance as that of a game with agents playing completely randomly (symmetric phase for m < m c , and then phase transition and symmetry breaking [START_REF] Marsili | Continuum time limit and stationary states of the minority game[END_REF]). However, it was emphasized that this characterization is insensitive to the introduction in the game of agents playing with a random history [START_REF] Rajković | Quantifying complexity in the minority game[END_REF]. In this variant of the game, where a proportion of agents play using the real history, and the other agents play using a random history, the behavior of σ 2 1 remains unchanged as the proportion of "random agents" increases. Rajkovic showed in [START_REF] Rajković | Quantifying complexity in the minority game[END_REF] that the Lempel-Ziv complexity of the minority sequence c(µ) shows the opposite behavior to volatility, with the same critical size of memory, for MG. But as the proportion of agents playing with a random history increases, c(µ) exhibits the modified nature of the game. Intuitively, this is due to the fact that if all the agents play with a random history, there is no time-organization in the game.

A multidimensional extension of the game has been proposed in [START_REF] Földy | Hierarchically organized minority games[END_REF] using a three-level structure. It consists of several MG in parallel, but where a connection is made through the payoffs of the strategies, by taking into account the minorities of the whole game. Here, to show how a multidimensional LZC can be used in the MG, we first propose the following simpler multidimensional variant of the game. We consider N players playing with N a assets. For asset j, agent i makes a decision using the m last minorities of asset j, but also taking into account the m a last minorities of K other assets. K will be called the connectivity of the multidimensional minority game (MMG). The K connections k i,j,0 , . . . k i,j,K-1 ∈ {0, . . . , N a -1} of agent i and variable j are randomly chosen in the initiation of the game and an asset can possibly be connected to itself (in this case there is an index l so that k i,j,l = j). Then at any time t, the decision of agent i for the asset j is made using the data (µ j (tm), . . . , µ j (t -1), µ k i,j,0 (tm a ), . . . , µ k i,j,0 (t -1), . . . , µ k i,j,K-1 (t -1)). For K = 0, it is as if N a MG were played independently, but for K = 0, there are clearly spatiotemporal links. As for the standard MG, for each asset the agents possess s strategies in the functions {0, 1} m+Kma -→ {0, 1}. Each strategy for each asset and of each agent possesses a gain that evolves in time (one point won for a good strategy, 0 otherwise).

In our illustrations, the MMG is composed of N = 101 agents, playing on N a = 5 assets, with s = 2 strategies for each agent and each asset. An additional memory m a = 1 and m a = 2 for each connected asset is considered, and the time duration of the game is 10000 samples. Figure 3 ; the normalized joint LZC of all the assets c(µ 0 , . . . , µ Na-1 ) log 2 (T ) T is plotted figure 6 for m a = 1 and m a = 2; Lastly, the behavior of the normalized informational LZC of the first two assets I c (µ 0 , µ 1 ) log 2 (T ) T and of all the assets I c (µ 0 , . . . , µ Na-1 ) log 2 (T ) T is given in figure 7. In these figures, the solid line with stars represents the result given by a totally random game. For K = 0, since the first asset comes from a standard MG, the behavior of σ 2 1 N is that described in [START_REF] Challet | Emergence of cooperation and organization in an evolutionary game[END_REF][START_REF] Savit | Adaptive competition, market efficiency, and phase transition[END_REF]. As shown in [START_REF] Rajković | Quantifying complexity in the minority game[END_REF] the LZC c(µ 0 ) shows the opposite behavior, leading to an equivalent characterization of the game. c(µ 0 , µ 1 ) has also the same behavior. In this case, m a does not change the behavior of the curves, since with no connection this additional memory is not taken into account. As K increases, the volatility exhibits a behavior that is increasingly similar to that of a totally random game. Since the volatility is similar for the five games, the curve would have the same shape on average over the assets. Although this average variable still represents the total waste of the game, it is not efficient to capture the self-organization hidden in this MMG. The LZC of the first asset, c(µ 0 ), leads to the same conclusion, as well as the joint LZC of the first two assets, c(µ 0 , µ 1 ), and the joint LZC of all the assets, c(µ 0 , . . . , µ Na-1 ). Furthermore these characteristics seem insensitive to the additional memory m a . However, some fluctuations of the LZCs as a function of m can be seen in these curves. Since the joint LZCs and the marginal ones have the same shape, this suggests that the combination of marginal and joint LZC can reveal the self-organized characteristic of the proposed MMG. Figure 7 depicts the two particular combinations represented by the informational LZC: I c (µ 0 , µ 1 ) = c(µ 0 ) + c(µ 1 )c(µ 0 , µ 1 ) and I c (µ 0 , . . . , µ Na-1 ) = Na l=1 (-1) l+1 i 1 <...<i l c(µ i 1 , . . . , µ i l ) (see [START_REF] Matsuda | Physical nature of higher-order mutual information: Intrinsic correlations and frustration[END_REF] for the analogy with Shannon entropy). In both cases, it is clear that the behavior of the ILZC is far from that of a random game. Moreover, some differences can clearly be seen for different values of connectivity K. In conclusion, the ILZC can clearly reveal the self-organized behavior of the proposed MMG. Since the purpose of the paper is not to study in detail this MMG, no further analyzis is given here. However, it will be interesting to investigate if some critical memory sizes in this MMG, taking into account both m and m a or not, would be exhibited through the ILZC, as is the case for the LZC of a standard MG with the memory m [START_REF] Marsili | Continuum time limit and stationary states of the minority game[END_REF].

Discussion

The Lempel-Ziv complexity has long been known and has been extensively used in the compression domain. However, this tool is not so extensively employed for data analysis. Recent studies present this tool as a potential tool for both analyzing biomedical sequences [START_REF] Zhang | Detecting ventricular tachycardia and fibrillation by complexity measure[END_REF][START_REF] Radhakrishnan | Quantifying physiological data with Lempel -Ziv complexity -certain issues[END_REF][START_REF] Szczepański | Application of Lempel -Ziv complexity to the analysis of neural discharges[END_REF][START_REF] Zhang | EEG complexity as a measure of depth of anesthesia for patients[END_REF] and signals from complex systems such as the minority game [START_REF] Rajković | Quantifying complexity in the minority game[END_REF]. In this paper we show that the LZC can naturally be understood for vectorial sequences, and then used to analyze multidimensional signals. We propose, then, derived quantities such as informational LZC, by analogy with Shannon entropy, and we have shown that these derived quantities can also be used for data analysis. Notice, however, that the use of the LZCs seems interesting for the analysis of signals which are naturally coded on a finite size alphabet. It is the case for random boolean networks or for the different variants of the minority game. For real world signals, one can envisage analyzing DNA sequences that are naturally coded on an alphabet of size 4. But it seems at present more difficult to use this tool for continuous-state signals such as ECG, even if some attempts have already been made. In the latter case, a crucial point is the quantification of the signals: static or dynamic quantization [START_REF] Ebeling | Entropy, complexity, predictability and data analysis of time series and letter sequences[END_REF]? How many levels of quantization [START_REF] Ebeling | Entropy, complexity, predictability and data analysis of time series and letter sequences[END_REF]? Very often, empirical rules have been chosen [START_REF] Zhang | Detecting ventricular tachycardia and fibrillation by complexity measure[END_REF][START_REF] Zhang | EEG complexity as a measure of depth of anesthesia for patients[END_REF], but the choice may have consequences on the LZC of the quantized signals. Working on these points seems very challenging. Another challenge would be also to investigate links between signals of different natures, such as electrocardiograms, electromyograms (etc.) via the joint LZC or derived quantities.

depicts the normal- ized volatility σ 2 1 N

 1 , as a function of the memory m, considering only the first asset of the MMG, in the case of an additional connecting memory m a = 1 and m a = 2; Figure4describes the normalized LZC of the minority sequence for the first asset c(µ 0 ) log 2 (T ) T ; Figure5represents the normalized LZC of the minority sequence for the first asset c(µ 0 ) log 2 (T ) T

Fig. 1 .Fig. 2 .

 12 Fig. 1. Mean joint LZC for RBNs, as a function of the bias p and the connectivity K. The RBNs contain N = 50 cells, the duration is T = 50 and c has been averaged using 100 realizations. The solid line represents the critical line K c = 1 2p(1-p) . The two figures on the right give two particular illustrations of an RBN, respectively for (K, p) = (8, .7) and (4, .9). The black (resp. white) points correspond to the value 1 (resp. 0).

reconstruction of the phase space, with several estimations to determine the embedding dimension and the optimal delay; then, estimation of the whole Lyapunov spectrum or just of some exponents (positive, max. . . )

Notice that it is the idea used in many well known compression algorithms, such as "gzip". Other parsing schemes also exist[START_REF] Ziv | Compression of individual sequences via variable-rate coding[END_REF], leading to several variants of compression algorithms.


David Monniaux

Laure Gonnord

Using Bounded Model Checking to Focus Fixpoint Iterations *

Two classical sources of imprecision in static analysis by abstract interpretation are widening and merge operations. Merge operations can be done away by distinguishing paths, as in trace partitioning, at the expense of enumerating an exponential number of paths.

In this article, we describe how to avoid such systematic exploration by focusing on a single path at a time, designated by SMT-solving. Our method combines well with acceleration techniques, thus doing away with widenings as well in some cases. We illustrate it over the well-known domain of convex polyhedra.

Introduction

Program analysis aims at automatically checking that programs fit their specifications, explicit or not -e.g. "the program does not crash" is implicit. Program analysis is impossible unless at least one of the following holds: it is unsound (some violations of the specification are not detected), incomplete (some correct programs are rejected because spurious violations are detected), or the state space is finite (and not too large, so as to be enumerated explicitly or implicitly). Abstract interpretation is sound, but incomplete: it over-approximates the set of behaviours of the analysed program; if the over-approximated set contains incorrect behaviours that do not exist in the concrete program, then false alarms are produced. A central question in abstract interpretation is to reduce the number of false alarms, while keeping memory and time costs reasonable [START_REF] Blanchet | A static analyzer for large safety-critical software[END_REF].

Our contribution is a method leveraging the improvements in SMT-solving to increase the precision of invariant generation by abstract fixpoint iterations. On practical examples from the literature and industry, it performs better than previous generic technique and is less "ad-hoc" than syntactic heuristics found in some pragmatic analyzers.

Listing 1: C implementation of y = sin(x)/x -1, with the -0.01 ≤ x ≤ 0.01 range implemented using a Taylor expansion around zero in order to avoid loss of precision and division by zero as sin(x) ≃ x → 0. The first source of imprecision in abstract interpretation is the choice of the set of properties represented inside the analyser (the abstract domain). Obviously, if the property to be proved cannot be reflected in the abstract domain (e.g. we wish to prove a numerical relation but our abstract domain only considers Boolean variables), then the analysis cannot prove it.

In order to prove that there cannot be a division by zero in the first branch of the second if-then-else of Listing 1, one would need the non-convex property that x ≥ 0.01 ∨ x ≤ -0.01. An analysis representing the invariant at that point in a domain of convex properties (intervals, polyhedra, etc.) will fail to prove the absence of division by zero (incompleteness).

Obviously, we could represent such properties using disjunctions of convex polyhedra, but this leads to combinatorial explosion as the number of polyhedra grows: at some point heuristics are needed for merging polyhedra in order to limit their number; it is also unclear how to obtain good widening operators on such domains. The same expressive power can alternatively be obtained by considering all program paths separately ("merge over all paths") and analysing them independently of each other. In order to avoid combinatorial explosion, the trace partitioning approach [START_REF] Rival | The trace partitioning abstract domain[END_REF] applies merging heuristics. In contrast, our method relies on the power of modern SMT-solving techniques.

The second source of imprecision is the use of widening operators [START_REF] Cousot | Abstract interpretation frameworks[END_REF]. When analysing loops, static analysis by abstract interpretation attempts to obtain an inductive invariant by computing an increasing sequence X 1 , X 2 , . . . of sets of states, which are supersets of the sets of states reachable in at most 1, 2, . . . iterations. In order to enforce convergence within finite time, the most common method is to use a widening operator, which extrapolates the first iterates of the sequence to a candidate limit. Optional narrowing iterations may regain some precision lost by widening.

Illustrating Example Consider Listing 2, a simplification of a fragment of an actual industrial reactive program: indexing of a circular buffer used only at certain iterations of the main loop of the program, chosen non-deterministically. If the non-deterministic choice nondet() is replaced by true, analysis with widening and narrowing finds [0, 99]. Unfortunately, the "narrowing" trick is brittle, and on 2 Listing 2: Circular buffer indexing

i n t x = 0 ; while ( t r u e ) { i f ( nondet ( ) ) { x = x + 1 ; i f ( x >= 1 0 0 ) x = 0 ; } }
Listing 2, widening yields [0, +∞), and this is not improved by narrowing! 1 In contrast, our semantically-based method would compute the [0, 99] invariant on this example by first focusing on the following path inside the loop: Listing 3: Example focus path assume ( nondet ( ) ) ; x = x + 1 ; assume ( x < 1 0 0 ) ; If we wrap this path inside a loop, then the least inductive invariant is [0, 99]. We then check that this invariant is inductive for the original loop. This is the basic idea of our method: it performs fixpoint iterations by focusing temporarily on certain paths in the program. In order to obtain the next path, it performs bounded model checking using SMT-solving. x := 0 We consider programs defined by a control flow graph: a set P of control points, for each control point p ∈ P a (possibly empty) set I p of initial values, a set E ⊆ P × P of directed edges, and the semantics τ e : P(Σ) → P(Σ) of each edge e ∈ E where P(Σ) is the set of possible values of the tuple of program variables. τ e thus maps a set of states before the transition expressed by edge e to the set of states after the transition.

Background and Notations in Abstract Interpretation

x := x + 1 x ≥ 100 x := 0 x < 100 (a) With original variables p 1 p 2 : x 2 = φ(x 1 , x 2 , x 3 , x 4 ) p 3 x 1 = 0 x 3 := x 2 + 1 x 3 ≥ 100 x 4 = 0 x 3 < 100 (b) SSA version. x = φ(e 1 ,
To each control point p ∈ P we attach a set X p ⊆ Σ of reachable values of the tuple of program variables at program point p. The concrete semantics of the program is the least solution of a system of semantic equations [START_REF] Cousot | Abstract interpretation frameworks[END_REF]:

X p = I p ∪ (p ′ ,p)∈E τ (p ′ ,p) (X p ′ ).
Abstract interpretation replaces the concrete sets of states in P(Σ) by elements of an abstract domain D. In lieu of applying exact operations τ to sets of concrete program states, we apply abstract counterparts τ ♯ . 2 An abstraction τ ♯ of a concrete operation τ is deemed to be correct if it never "forgets" states:

∀X ∈ D τ(X) ⊆ τ ♯ (X) (1)
We also assume an "abstract union" operation ⊔, such that X ∪ Y ⊆ X ⊔ Y. For instance, Σ can be Q n , D can be the set of convex polyhedra and ⊔ the convex hull operation [START_REF] Halbwachs | Détermination automatique de relations linéaires vérifiées par les variables d'un programme[END_REF][START_REF] Cousot | Automatic discovery of linear restraints among variables of a program[END_REF][START_REF] Bagnara | The Parma Polyhedra Library[END_REF].

In order to find an inductive invariant, one solves a system of abstract semantic inequalities:

       ∀p I p ⊆ X p ∀(p ′ , p) ∈ E τ ♯ (p ′ ,p) (X p ′ ) ⊆ X p . (2)
Since the τ ♯ e are correct abstractions, it follows that any solution of such a system defines an inductive invariant; one wishes to obtain one that is as strong as possible ("strong" meaning "small with respect to ⊆"), or at least sufficiently strong as to imply the desired properties.

Assuming that all functions τ ♯ e are monotonic with respect to ⊆, and that ⊔ is the least upper bound operation in D with respect to ⊆, one obtains a system of monotonic abstract equations:

X p = I p ⊔ (p ′ ,p)∈E τ ♯ (p ′ ,p) (X p ′ ). If (D, ⊆) has no infinite ascending sequences (d 1 d 2 . . . with d 1 , d 2 , • • • ∈ D)
, then one can solve such a system by iteratively replacing the contents of the variable on the left hand side by the value of the right hand side, until a fixed point is reached. The order in which equations are iterated does not change the final result.

Many interesting abstract domains, including that of convex polyhedra, have infinite ascending sequences. One then classically uses an extrapolation operator known as widening and denoted by ▽ in order to enforce convergence within finite time. The iterations then follow the "upward iteration scheme":

X p := X p ▽          X p ⊔ (p ′ ,p)∈E τ ♯ (p ′ ,p) (X p ′ )          (3)
where the contents of the left hand side gets replaced by the value of the right hand side. The convergence property is that any sequence u n of elements of D of the form u n+1 = u n ▽ v n , where v n is another sequence, is stationary [START_REF] Cousot | Abstract interpretation frameworks[END_REF]. It is sufficient to apply widening only at a set of program control nodes P W such that all cycles in the control flow graph are cut. Then, through a process of chaotic iterations [13, Def. 4.1.2.0.5, p. 127], one converges within finite time to an inductive invariant satisfying Rel. 2.

Once an inductive invariant is found, it is possible to improve it by iterating the ψ ♯ function defined as Y = ψ ♯ (X), noting X = (X p ) p∈P and Y = (Y p ) p∈P , with

Y p = I p ⊔ (p ′ ,p)∈E τ ♯ (p ′ ,p) (X p ′ ).
If X is an inductive invariant, then for any k, ψ ♯ k (X) is also an invariant. This technique is an instance of narrowing iterations, which may help recover some of the imprecision induced by widening [14, §4]. Choose p 1 ∈ A 7:

Algorithm 1 Classical Algorithm

A ← A \ {p 1 } 8:
for all outgoing edge (e) from p 1 do 9:

Let p 2 be the destination of e : A naive implementation of the upward iteration scheme described above is to maintain a work-list of program points p such that X p has recently been updated and replaced by a strictly larger value (with respect to ⊆), pick and remove the foremost member p, apply the corresponding rule X p := . . . , and insert into the work-list all p ′ such that (p, p ′ ) ∈ E (This algorithm is formally described in Algorithm 1). 'd) Figure 1(a) gives the control flow graph obtained by compilation of Listing 2. Node p 2 is the unique widening node.

Example of Section 1 (Cont

The classical algorithm (with the interval abstract domain) performs on this control flow graph of the following iterations :

• Initialisation : X p 1 ← (-∞, +∞), X p 2 ← X p 3 ← X p 4 ← ∅. • Step 1: X p 2 ← [0, 0], then the transition to p 3 is enabled, X p 3 ← [1, 1], then
the return edge to p 2 gives the new point x = 1 to X p 2 , the new polyhedron is then X p 2 = [0, 1] after performing the convex hull. Widening gives the polyhedron

X p 2 = [0, ∞).
(The widening operator on intervals is defined as

[x l , x r ] ▽[x ′ l , x ′ r ] = [x" l , x" r ] where x" l = x l if x l = x ′ l else -∞, and x" r = x r if x r = x ′ r else +∞.) • Step 2: X p 3 becomes [1, +∞).
The second transition from p 3 to p 2 is thus enabled, and the back edge to p 2 gives the point x = 0 to X p 2 . At the end of step 2 the convergence is reached.

• If we perform a narrowing sequence, there is no gain of precision because of the simple loop over the control point p 2 .

Our Method

We have seen two examples of programs where classical polyhedral analysis fails to compute good invariants. How could we improve on these results?

• In order to get rid of the imprecision in Listing 1, one could "explode" the control-flow graph: in lieu of a sequence of n if-then-else, with n merge nodes with 2 input edges, one could distinguish the 2 n program paths, and having a single merge node with 2 n input edges. As already pointed out, this would lead to exponential blowup in both time and space.

• One way to get rid of imprecision of classical analysis (Sec. 2) on the program from Fig. 1(a) would be to consider each path through the loop at a time and compute a local invariant for this path. Again, the number of such paths could be exponential in the number of tests inside the loop.

The contribution of our article is a generic method that addresses both of these difficulties.

Reduced Transition Multigraph and Path Focusing

Consider a control flow graph (P, E) with associated transitions (τ e ) e∈E , a set of widening points P W ⊆ P such that removing P W cuts all cycles in the graph, and a set P R of abstraction points, such that P W ⊆ P R ⊆ P (On the figures, the nodes in P R are in bold). We make no assumption regarding the choice of P W ; there are classical methods for choosing widening points [9, §3.6]. P R can be taken equal to P W , or may include other nodes; this makes sense only if these nodes have several incoming edges. Including other nodes will tend to reduce precision, but may improve scalability. We also make the simplifying assumption that the set of initial values I p is empty for all nodes in P \ P R -in other words, the set of possible control points at program start-up is included in P R .

We construct (virtually) the reduced control multigraph (P R , E R ), with edges E R consisting of the paths in (P, E) that start and finish on nodes in P R , with associated semantics the composition of the semantics of the original edges

τ e 1 →•••→e n = τ e n • • • • • τ e 1 .
There are only a finite number of such edges, because the original graph is finite and removing P R cuts all cycles. There may be several edges between two given nodes, because there may exist several control paths between these nodes in the original program. Equivalently, this multigraph can be obtained by starting from the original graph (P, E) and by removing all nodes p in P \ P R as follows: each couple of edges e 1 , from p 1 to p, and e 2 , from p to p 2 , is replaced by a single edge from p 1 to p 2 with semantics τ p 2 • τ p 1 .

Example of Section 1 (Cont'd)

The reduced control flow graph obtained for our running example is loop

x := 0 guard x ≥ 99 x := 0 guard x < 99 x := x + 1
Our analysis algorithm performs chaotic iterations over that reduced multigraph, without ever constructing it explicitly. We start from an iteration strategy, that is, a method for choosing which of the equations to apply next; one may for instance take a variant of the naive "breadth-first" algorithm from §2, but any iteration strategy [9, §3.7] befits us (see also Alg. 1). An iteration strategy maintains a set of "active nodes", which initially contains all nodes p such that I p ∅. It picks one edge e from an active node p 1 to a node p 2 , and applies X p 2 := X p 2 ⊔ τ ♯ e (X p 1 ) in the case of a node p 2 ∈ P R \P W , and applies X p 2 := X p 2 ▽(X p 2 ⊔τ ♯ e (X p 1 )) if p 2 ∈ P W ; then p 2 is added to the set of active nodes if the value of X p 2 has changed.

Our alteration to this algorithm is that we only pick edges e from p 1 to p 2 such that there exist x 1 ∈ X p 1 , x 2 ∈ τ e ({x 1 }) and x 2 X p 2 with the current values of X p 1 and X p 2 . In other words, going back to the original control flow graph, we only pick paths that add new reachable states to their end node, and we temporarily focus on such a path.

How do we find such edges e out of potentially exponentially many? We express them as the solution of a bounded reachability problem -how can we go from control state p 1 with variable state in X p 1 to control state p 2 with variable state in X p 2 -, which we solve using satisfiability modulo theory (SMT). (See Alg. 2)

Finding Focus Paths

We now make the assumption that both the program transition semantics τ e and the abstract elements x ♯ ∈ D can be expressed within a decidable theory T (this assumption may be relaxed by replacing the concrete semantics, including e.g. multiplicative arithmetic, by a more abstract one through e.g. linearization [START_REF] Miné | Symbolic methods to enhance the precision of numerical abstract domains[END_REF]).

Such is for instance the case if the program operates on rational values, so a program state is an element of Σ = Q n , all operations in the program, including guards and assignments, are linear arithmetic, and the abstract domain is the domain of convex polyhedra over Q n , in which case T can be the theory of linear real arithmetic (LRA). If program variables are integer, with program state space Σ = Z n , but still retaining the abstract domain of convex polyhedra over Q n , then we can take T to be the theory of linear integer arithmetic (LIA). Deciding the satisfiability of quantifier-free formulas in either LIA or LRA, with atoms consisting in propositional variables and in linear (in)equalities with integer coefficients, is NP-complete. There however exist efficient decision procedures for such formulas, known as SMT-solvers, as well as standardised theories and file formats [START_REF] Barrett | The satisfiability modulo theories library (SMT-LIB)[END_REF]; notable examples of SMT-solvers capable of dealing with LIA and LRA are Z3 and Yices. Kroening & Strichman [START_REF] Kroening | Decision procedures[END_REF] give a good introduction to the techniques and algorithms in SMT solvers.

We assume that the program is expressed in SSA form, with each program variable being assigned a value at only a single point within the program [START_REF] Cytron | An efficient method of computing static single assignment form[END_REF]; standard techniques exist for converting to SSA. Figure 1 gives both "normal" and SSA-form control-flow graphs for Listing 2.

We transform the original control flow graph (P, E) in SSA form by disconnecting the nodes in P R : each node p r in P R is split into a "source" node p s r with only outbound edges, and a "destination" node p d r with only inbound edges. We call the resulting graph (P ′ , E ′ ). Figure 2(a) gives the disconnected SSA form graph for Listing 2 where p 1 and p 2 have been split.

We consider execution traces starting from a p s r node and ending in a p d r node. We define them as for doing bounded model checking [START_REF] Armando | Bounded model checking of software using smt solvers instead of sat solvers[END_REF]. To each node p ∈ P ′ we attach a Boolean b p or reachability predicate, expressing that the trace goes through program point p. For nodes p ′ not of the form p s r , we have a constraint b p ′ = p e p,p ′ , for e p,p ′ ranging over all incoming edges. To each edge p → p ′ we attach a Boolean e p,p ′ , and a constraint e p,p ′ = b p ∧ τ p,p ′ . The conjunction ρ of all these constraints, expresses the transition relation between the p s r and p d r nodes (with implicit existential quantification).

If the transitions τ (p,p ′ ) are non-deterministic, a little care must be exercised for the path obtained from the b p to be unique. For instance, if from program point p 1 one can move non-deterministically to p 2 or p 3 through edges e 2 and e 3 an incorrect way of writing the formula would be (b 2 = e 2 )∧(b 3 = e 3 )∧(e 2 = b 1 )∧(e 3 = b 1 ), in which case b 2 and b 3 could be simultaneously true. Instead, we introduce special "choice" variables c i that model non-deterministic choices (Fig. 2).

In order to find a path from program point p 1 ∈ P R , with variable state x 1 , to program point p 2 ∈ P R , with variable state x 2 , we simply conjoin ρ with the formulas x 1 ∈ X p 1 and x 2 X p 2 , with x 1 , x 2 , x 1 ∈ X p 1 and x 2 X p 2 expressed in terms of the SSA variables. 3 For instance, if X p 1 and X p 2 are convex polyhedra defined by systems of linear inequalities, one simply writes these inequalities using the names of the SSA-variables at program points p 1 and p 2 .

p s 1 p s 2 p 3 p d 2 : x ′ 2 = φ(x 1 , x 4 , x 3 , x 2 )
e 3 x 3 = x 2 + 1 e 5 x 3 ≥ 100 x 4 = 0 e 4 x 3 < 100 e 2 e 1 x 1 = 0 (a) Disconnected (SSA) CFG p s 1 p s 2 p 3 p d 2 : x ′ 2 = φ(x 1 , x 4 , x 3 , x 2 )
e 3 x 3 = x 2 + 1
We apply SMT-solving over that formula. The result is either "unsatisfiable", in which case there is no path from p 1 , with variable values x 1 , to p 2 , with variable values x 2 , such that x 1 ∈ X p 1 and x 2 X p 2 , or "satisfiable", in which case SMTsolving also provides a model of the formula (a satisfying assignment of its free variables); from this model we easily obtain such a path, unique by construction of ρ.

Indeed, a model of this formula yields a trace of execution: those b p predicates that are true designate the program points through which the trace goes, and the other variables give the values of the program variables.

Example of Section 1 (Cont'd)

The SSA form of the control flow graph of Figure 1(a) is depicted in Figure 1(b). Fig. 2 shows the disconnected version of the SSA Graph (the node p 2 is now split), and the formula ρ expressing the semantics is shown beneath it.

Then, consider the problem of finding a path starting in control point 2 inside polyhedron x = 0 and ending at the same control point but outside of that polyhe-dron. Note that because there are two outgoing transitions from node p s 2 , which are chosen non-deterministically, we had to introduce a Boolean choice variable c s 2 . The focus path of Fig. 2

Algorithm

Algorithm 2 consists in the iteration of the path finding method of Sec. 3.2, coupled with forward abstract interpretation along the paths found and, optionally, path acceleration.

Correctness and Termination

We shall now prove that this algorithm terminates, and that the resulting X p define an inductive invariant that contains all initial states I p . The proof is a variant of the correctness proof of the chaotic iterations.

The invariant maintained by this algorithm is that all nodes p 1 ∈ P R \A are such that there is no execution trace starting at point p 1 in a state x 1 ∈ X p 1 and ending at point p 2 in a state x 2 X p 2 . Evidently, if A becomes empty, then this condition means that X p is an inductive invariant.

Termination is ensured by the classical argument of termination of chaotic iterations in the presence of widening: they always terminate if all cycles in the control flow graph are broken by widening points [13, Th. 4.1.2.0.6, p. 128]. In short, an infinite iteration sequence is bound to select at least one node p in P W an infinite amount of times, because P W breaks all cycles, but due to the properties of widening, X p should be stationary, which contradicts the infinite number of selections. Our comment at line 20 of Alg. 2 is important for termination: it ensures that for any widening node p, the sequence of values taken by X p when it is updated and reinserted into set A is strictly ascending, which ensures termination in finite time.

Self-Loops

The algorithm in the preceding subsection is merely a "clever" implementation of standard polyhedral analysis [START_REF] Cousot | Automatic discovery of linear restraints among variables of a program[END_REF][START_REF] Halbwachs | Détermination automatique de relations linéaires vérifiées par les variables d'un programme[END_REF] on the reduced control multigraph (P R , E R ); the difference with a naive implementation is that we do not have to explicitly enumerate an exponential number of paths and instead leave the choice of the focus path to the SMT-solver. We shall now describe an improvement in the case of selfloops, that is, single paths from one node to itself.

Algorithm 3 is a variant of Alg. 2 where self-loops are treated specially:

• The loopiter(τ ♯ , X) function returns the result of a widening / narrowing iteration sequence for abstract transformer τ ♯ starting in X; it returns X ′ such Choose p 1 ∈ A 10:

A ← A \ {p 1 } 11: repeat 12: res ← SmtSolve          ρ ∧ b p 1 ∧ x 1 ∈ X p 1 ∧ p 2 |(p 1 ,p 2 )∈E ′ b p 2 ∧ x 2 X p 2          13:
if res is not "unsat" then

X temp ← X p 2 ▽ X p 2 ⊔ Y ⊲ Final point p 2 is a widening point 18:
else 19:

X temp ← X p 2 ⊔ Y 20:
end if ⊲ at this point X temp X p 2 otherwise p 2 would not have been chosen until res="unsat" 25: end while ⊲ End of Iteration 26: Possibly narrow (see Sec. 4.1) 27: Compute X p i for p i P R 28: return all X p i that X ⊆ X ′ and τ ♯ (X ′ ) ⊆ X ′ .

• In order not to waste the precision gained by loopiter, the first time we consider a self-loop e ′ we apply a union operation instead of a widening; set U records the self-loops that have already been visited. This is a form of delayed widening [START_REF] Halbwachs | Delay analysis in synchronous programs[END_REF].

Termination is still guaranteed, because the inner loop cannot loop forever: it can visit any self-loop edge e ′ at most once before applying widening. Choose p 1 ∈ A 10:

A ← A \ {p 1 } 11: U = ∅ ⊲ U is a set of "
res ← SmtSolve          ρ ∧ b p 1 ∧ x 1 ∈ X p 1 ∧ p 2 |(p 1 ,p 2 )∈E ′ b p 2 ∧ x 2 X p 2          14:
if res is not "unsat" then

Compute e ′ ∈ E ′ from res 16:

if p 1 = p 2 then 17: Y ← loopiter(τ ♯ e ′ , X p 1 ) 18: else 19:
Y ← τ ♯ e ′ (X p 1 )

X p 2 ← X p 2 ⊔ Y 25: U ← U ∪ {e ′ } 26:
end if

27:

A ← A ∪ {p 2 }

28:

end if

29:

until res="unsat" 30: end while ⊲ End of Iteration 31: Compute X p i s for p i P R 32: return all X p i s

Example of Section 1 (Cont'd)

Let us perform our algorithm on our example :

• Step 1 : Is there a path from control point p 1 to control point p 2 feasible (without additional constraint) ? Yes. On Figure 2, the obtained model corresponds to the transition from p s 1 to p d 2 , and leads to the interval X p 2 = [0, 0]. • Step 2 : Is there a path from p 2 with x = 0 to p 2 with x 0 ? The answer to this query is depicted in Figure 2(b): there is such a path, on which we now focus. This path is considered as a loop and we therefore do a local iteration with widenings (loopiter). X p 2 becomes [0, 1], then after widening [0, ∞]. A narrowing step gives finally X p 2 = [0, 99], which is thus the result of loopiter.

• Step 3 : Is there a path from p 2 with x ∈ [0, 99] to p 2 with x ′ [0, 99] ? No.

The iteration thus ends with the desired invariant.

Extensions

Narrowing

Narrowing iterations can also be applied within our framework. Let us assume that some inductive invariant X p∈P R has been computed; it satisfies the relation ψ(X) ⊆ X component-wise, noting X = (X 1 , . . . , X |P| ), and ψ(X) denotes (Y 1 , . . . , Y |P| ) defined as

Y p 2 = I p 2 ∪ e∈E R e from p 1 to p 2 τ e X p 1 (4)
The abstract counterpart to this operator is ψ ♯ , defined similarly, replacing τ by τ ♯ and ∪ by ⊔. It satisfies the correctness condition (see Rel. 1) ∀X ∈ D ψ(X) ⊆ ψ ♯ (X).

As per the usual narrowing iterations, we compute a narrowing sequence X (k) = ψ ♯ k (X). It is often sufficient to stop at k = 1; otherwise one may stop when X (k+1) X (k) . Let us now see a practical algorithm for computing Y = ψ ♯ (X):

For all p ∈ P R , we initialise Y p := I p . For all p 2 ∈ P R , we consider all paths e ∈ E R from p 1 ∈ P R to p 2 such that there exist x 1 ∈ X p 1 , x 2 ∈ X p 2 , x 2 ∈ τ e ({x 1 }) as explained in §3.2. We then update Y p 2 := Y p 2 ⊔ τ ♯ e (X p 1 ).

Acceleration

In Sec. 3.5, we have described loopiter function that performs a classical widening / narrowing iteration over a single path. In fact, the only requirement over it is that loopiter(τ ♯ , X) returns X ′ such that X ⊆ X ′ and τ ♯ (X ′ ) ⊆ X ′ . In other words, X ′ is an over-approximation of τ ♯ * (X), noting R * the transitive closure of R.

In some cases, we can compute directly such an over-approximation, sometimes even obtaining τ ♯ * (X) exactly; this is known as acceleration of the loop. Examples of possible accelerations include the case where τ e is given by a difference bound matrix [START_REF] Comon | Multiple counters automata, safety analysis and Presburger arithmetic[END_REF], an octagon [START_REF] Bozga | Iterating octagons[END_REF], ultimately periodic integer relations [START_REF] Bozga | Fast acceleration of ultimately periodic relations[END_REF] or certain affine linear relations [START_REF] Gonnord | Combining widening and acceleration in linear relation analysis[END_REF][START_REF] Gonnord | Accelération abstraite pour l'amélioration de la précision en analyse des relations linéaires[END_REF][START_REF] Ancourt | A modular static analysis approach to affine loop invariants detection[END_REF].

For instance, the focus path of Fig. 2(b) consists in the operations and guards x = x + 1; x < 100; instead of iterating that path, we can compute its exact acceleration, yielding x ∈ [0, 99].

Partitioning

It is possible to partition the states at a given program point according to some predicate or a partial history of the computation [START_REF] Rival | The trace partitioning abstract domain[END_REF]. This amounts to introducing several graph nodes representing the same program point, and altering the transition relation.

Input-Output Relations

As with other analyses using relational domains, it is possible to obtain abstractions of the input-output relation of a program block or procedure instead of an abstraction of the set of states at the current point [START_REF] Ancourt | A modular static analysis approach to affine loop invariants detection[END_REF]; this also allows analyzing recursive procedures [START_REF] Halbwachs | Détermination automatique de relations linéaires vérifiées par les variables d'un programme[END_REF]Sec. 7.2]. It suffices to include in the set of variables copies of the variables at the beginning of the block or procedure; then the abstract value obtained at the end of the block or procedure is the desired abstraction.

Implementation and Preliminary Results

Our algorithm has been implemented as an option for Aspic, that computes invariants from counter automata with Linear Relation Analysis ( [START_REF] Feautrier | Accelerated invariant generation for C programs with Aspic and C2fsm[END_REF]). We wrote an Ocaml interface to the Yices SMT-solver ( [START_REF] Dutertre | A fast linear-arithmetic solver for DPLL(T)[END_REF]), and modified the fixpoint computation inside Aspic to deal with local iterations of paths. The implementation still needs some improvements, but the preliminary results are promising, and we describe some of them in Table 1. We provide no timing results since we were unable to detect any overcost due to the method. These two examples show that since we avoid (some) convex hulls, the precision of the whole analysis is improved.

The rate limiter example is particularly interesting, since, like the one in Listing 1 (which does not include a loop), it will be imprecisely analyzed by any method enforcing convex invariants at intermediate steps.

Related Work

Our algorithm may be understood as a form of chaotic iterations [13, §2.9.1, p. 53] over a certain system of semantic questions; we use SMT as an oracle to know which equations need propagating. The choice of widening points, and the order in which to solve the abstract equations, have an impact on the precision of the whole analysis, as well as its running time. Even though there exist few hard general results as to which strategy is best [13, §4.1.2, p. 125], some methods tend to experimentally behave better [START_REF] Bourdoncle | Sémantique des langages impératifs d'ordre supérieur et interprétation abstraite[END_REF].

"Lookahead widening" [START_REF] Gopan | Reps. Lookahead widening[END_REF] was our main source of inspiration: iterations and widenings are adapted according to the discovery of new feasible paths in the program. This approach avoids loss of precision due to widening in programs with

{ x = inp u t ( -1 0 0 0 , 1 0 0 0 ) ; i f ( x >= x_old +1) x = x_old + 1 ; i f ( x <= x_old -1) x = x_old -1; x_old = x ; } } ( 1 )
Source : [START_REF] Monniaux | Automatic modular abstractions for template numerical constraints[END_REF] In order to properly analyse such a program, Astrée distinguishes all four execution paths inside the loop through trace partitioning [START_REF] Rival | The trace partitioning abstract domain[END_REF], which is triggered by ad hoc syntactic criteria (e.g. two successive if-then-else). Our algorithm finds the invariant {-1000 ≤ x old ≤ 1000}, which is not found by classical analysis. multiple paths inside loops. It has proved its efficacy to suppress some gross overapproximations induced by naive widening. However, it does not solve the imprecisions introduced by convex hull (e.g. it produces false alarms on Listing 1).

Our method analyzes separately the paths between cut-nodes. We have pointed out that this is (almost) equivalent to considering finite unions of elements of the abstract domain, known as the finite powerset construction, between the cutnodes. 4 The finite powerset construction is however costly even for loop-free code, and it is not so easy to come up with widening operators to apply it to codes with loops or recursive functions [START_REF] Bagnara | Widening operators for powerset domains[END_REF]; for limiting the number of elements in the unions, some may be lumped together (thus generally introducing further over-approximation) according to affinity heuristics [START_REF] Sriram Sankaranarayanan | Static analysis in disjunctive numerical domains[END_REF][START_REF] Popeea | Inferring disjunctive postconditions[END_REF].

Still, in the recent years, much effort has been put into the discovery of disjunctive invariants, for instance in predicate abstraction [START_REF] Gulwani | Constraint-based invariant inference over predicate abstraction[END_REF]. Of particular note is the recent work by Gulwani and Zuleger on inferring disjunctive invariants [START_REF] Gulwani | The reachability-bound problem[END_REF] for finding bounds on the number of iterations of loops. We improve on their method on two points:

• In contrast to us, they assume that the transition relation is given in disjunctive normal form [26, Def. 5], which in general has exponential size in the number of tests inside the loop. By using SMT-solving, we keep the DNF implicit and thus avoid this blowup.

• By using acceleration, we may obtain more precise results than using widening, as they do for lattices that do not satisfy the ascending chain condition.

Nevertheless, their method allows expressing disjunctive invariants at loop heads, and not only at intermediate points, as we do. However, we think it is possible to get the best of both worlds and combine our method with theirs. In order to obtain a disjunctive invariant, they first choose a "convexity witness" (given that the number of possible witnesses is exponential, they choose it using heuristics) [26, p. 7], and then they compute a "transitive closure" [26, Fig. 6], which is a form of fixed point iteration of input-output relations (as in our Sec. 4.4) over an expanded control-flow graph. The choice of the convexity witness amounts to a partitioning of the nodes and transition (Sec. 4.3). Thus, it seems to possible to apply their technique, but replace their fixed point iteration [START_REF] Gulwani | The reachability-bound problem[END_REF]Fig. 6] by one based on SMT-solving and path focusing, using acceleration if possible.

In recent years, because of improvement in SMT-solving, techniques such as ours, distinguishing paths inside loops, have become tractable [START_REF] Monniaux | Automatic modular abstractions for linear constraints[END_REF][START_REF] Beyer | Path invariants[END_REF][START_REF] Monniaux | Automatic modular abstractions for template numerical constraints[END_REF][START_REF] Gawlitza | Improving strategies via SMT solving[END_REF]]. An alternative to using SMT-solving is to limit the number and length of traces to consider, as in trace partitioning [START_REF] Rival | The trace partitioning abstract domain[END_REF], used in the Astrée analyzer [START_REF] Cousot | Combination of abstractions in the astrée static analyzer[END_REF][START_REF] Cousot | Antoine Miné, David Monniaux, and Xavier Rival[END_REF][START_REF] Blanchet | A static analyzer for large safety-critical software[END_REF], but the criteria for limitation tend to be ad hoc. In addition, methods for abstracting the sets of paths inside a loop, weeding out infeasible paths, have been introduced [START_REF] Balakrishnan | Refining the control structure of loops using static analysis[END_REF].

With respect to optimality of the results, our method will generate the strongest inductive invariant inside the abstract domain if the domain satisfies the ascending chain condition and no widening is used; for other domains, like all methods using widenings, it may or may not generate it. In contrast, some recent works [START_REF] Gawlitza | Improving strategies via SMT solving[END_REF] guarantee to obtain the strongest invariant for the same analysis problem, at the expense of restriction to template linear domains and linear constructions inside the code.

Conclusion and future work

We have described a technique which leverages the bounded model checking capacities of current SMT solvers for guiding the iterations of an abstract interpreter. Instead of normal iterations, which "push" abstract values along controlflow edges, including control-flow splits and merges, we consider individual paths. This enables us, for instance, to use acceleration techniques that are not available when the program fragment being considered contains control-flow merges. This technique computes exact least invariants on some examples on which more conventional static analyzers incur gross imprecision or have to resort to syntactic heuristics in order to conserve precision.

We have focused on numerical abstractions. Yet, one would like to use similar techniques for heap abstractions, for instance. The challenge will then be to use a decidable logic and an abstract domain such that both the semantics of the program statements and the abstract values can be expressed in this logic. This is one direction to explore. With respect to the partitioning technique, 4.3, we currently express the partition as multiple explicit control nodes, but it seems desirable, for large partitions (e.g. according to Boolean values, as in B. Jeannet's BDD-Apron library) to represent them succinctly; this seems to fit nicely with our succinct encoding of the transition relation as a SMT-formula.

Another direction is to evaluate the scalability of these methods on larger programs. The implementation needs to be tested more to evaluate the precision of our method on middle-sized programs, the main advantage is that Aspic implements some of the acceleration techniques. Analyzers such as Astrée scale up to programs running a control loop several hundreds of thousands of lines long; translating such a loop to a SMT formula and solving for this formula and additional constraints does not seem tractable. It is possible that semantic slicing techniques [START_REF] Rival | Understanding the origin of alarms in Astrée[END_REF] could help in reducing the size of the generated SMT problems.

  i f ( x >= 0 ) { xabs = x ; } e l s e { xabs = -x ; } i f ( xabs >= 0 . 0 1 ) { y = s i n ( x ) / x -1 ; } e l s e { xsq = x * x ; y = xsq * ( -1/6. + xsq / 1 2 0 . ) ; }

  e 2 , . . . ) denotes a SSA φ-node: x takes value e 1 if control flows from the first incoming edge, e 2 from the second. . .

Figure 1 :

 1 Figure 1: Control flow graph corresponding to listing 2.

  With a focus path (solid edges) from x 2 = 0 at program point 2 to x ′ 2 = 1 at the same program point(e 1 = (x 1 = 0) ∧ b s 1 ) ∧ (e 3 = (x 3 = x 2 + 1) ∧ b s 2 ∧ c s 2 ) ∧ (e 2 = b s 2 ∧ ¬c s 2 ) ∧ (e 5 = b 3 ∧x 3 ≥ 100∧x 4 = 0)∧(e 4 = b 3 ∧x 3 < 100)∧(b 3 = e 3 )∧(b d 2 = e 1 ∨e 4 ∨e 5 ∨e 2 ) ∧ (x ′ 2 = ite(e 1 ,x 1 , ite(e 5 , x 4 , ite(e 4 , x 3 , x 2 )))) Figure2: Disconnected version of the SSA control flow graph of Fig.1(b), and the corresponding SMT formula. ite(b, e 1 , e 2 ) is a SMT construct whose value is "if b then the value of e 1 else the value of e 2 ". To each node p x corresponds a Boolean b x and an optional choice variable c x ; to each edge, a Boolean e y .

  (b) was obtained by solving the formula ρ ∧ b s 1 = false ∧ b s 2 = true ∧ b d 2 = true ∧ (x 2 = 0) ∧ ¬(x ′ 2 = 0): we impose that the path starts at point p s 2 (thus forcing b s 1 = false ∧ b s 2 = true) in the polyhedron x = 0 (thus x 2 = 0) and ends at point p d 2 (thus forcing b p 2 = true) outside of that polyhedron (thus ¬(x 2 = 0)).

18: end if 19: end for; 20: end while;

	10: 11:	if p 2 ∈ P W then X temp ← X p 2 ▽ X p 2 ⊔ τ ♯ e (X p 1 )	⊲ Widening node;
	12: 13:	else X temp ← X p 2 ⊔ τ ♯ e (X p 1 ) ;
	14:	end if
	15:	if X temp X p 2 then	⊲ The value must be updated
	16:	X p 2 ← X temp ;
	17:	A ← A ∪ {p 2 };
			⊲ End of Iteration
	21: Possibly narrow
	22: return all X p i s;

Algorithm 3

 3 Path-focused Algorithm with Self-Loops.

			marks changes from
	Alg. 2.
	1: Compute SSA-form of the control flow graph.
	2: Choose P R , compute the disconnected graph (P ′ , E ′ ) accordingly.
	3: ρ ← computeFormula(P ′ , E ′ )	⊲ Precomputations
	4: A ← ∅;
	5: for all p ∈ P R such that I p ∅ do
	6:	A ← A ∪ {p}
	7: end for;
	8: while A is not empty do	⊲ Fixpoint Iteration on the reduced graph
	9:

Table 1 :

 1 Invariant generation on two simple challenging programs

	Program	Automaton	Result and notes
			The compilation of the
			program gives an ex-
			panded control structure
			where some paths are
	Listing 4: Boustrophedon		"clearly"	unfeasible
	void boustrophedon ( ) { i n t x ; i n t d ; x = 0 ; d = 1 ; while ( 1 ) { i f ( x == 0 ) d = 1 ; i f ( x == 1 0 0 0 ) d= -1; x += d ; } }		(e.g. x < 0 and x > 1000), imposing both thus the only feasi-ble ones are guarded by x < 0, x = 0, 0 < x < 1000, x = 1000 and x > 1000. The tool finds the invariant {0 ≤ x ≤ 1000, -1 ≤ d ≤ 1} Classical Analysis with widening "upto" gives
			{d ≤ 1, d + 1999 ≥ 2x}
			and Gopan and Reps'
			improvement is not able
			to find x ≥ 0.
	Listing 5: Rate limiter
	void main ( ) {
	f l o a t x_old , x ;
	x_old = 0 ;
	while

On this example, it is possible to compute the [0, 99] invariant by so called "widening up-to" [28, Sec. 3.2], or with "thresholds"[START_REF] Blanchet | A static analyzer for large safety-critical software[END_REF]: essentially, the analyser notices syntactically the comparison x < 100 and concludes that 99 is a "good value" for x, so instead of widening directly to +∞, it first tries 99. This method only works if the interesting value is a syntactic constant.

Many presentations of abstract interpretation distinguish the abstract element x ♯ ∈ D from the set of states γ(x ♯ ) it represents. We opted not to, for the sake of brevity.

The formula defining the set of values represented by an abstract element X has sometimes been denoted by γ[START_REF] Reps | Symbolic implementation of the best transformer[END_REF].

It is equivalent if the only source of disjunctions are the splits in the control flow, and not atomic operations. For instance, if the test |x| ≥ 1 is considered an atomic operation, then we could take the disjunction x ≥ 1 ∨ x ≤ -1 as output. We can rephrase that as a control flow problem by adding a test x ≥ 0, otherwise said to express |x| as a piecewise linear function with explicit tests for splits between the pieces.

* This research was partially funded by ANR project "ASOPT".