Elodie L Plan
email: elodie.plan@farmbio.uu.se

Alan Maloney

France Mentré

Mats O Karlsson

Julie Bertrand

Performance Comparison of various Maximum Likelihood Nonlinear Mixed-effects Estimation Methods for Dose-Response Models AUTHORS

Keywords: MAXIMUM LIKELIHOOD ESTIMATION, FOCE, LAPLACE, ADAPTIVE GAUSSIAN

Estimation methods for nonlinear mixed-effects modelling have considerably improved over the last decades. Nowadays several algorithms implemented in different softwares are used. The present study aimed at comparing their performance for dose-response models.

Eight scenarios were considered using a sigmoid E max model, with varying sigmoidicity factors and residual error models. 100 simulated datasets for each scenario were generated. 100 individuals with observations at 4 doses constituted the rich design and at 2 doses for the sparse design. Nine parametric approaches for maximum likelihood estimation were studied: FOCE in NONMEM and R, LAPLACE in NONMEM and SAS, adaptive Gaussian quadrature (AGQ) in SAS, and SAEM in NONMEM and MONOLIX (both SAEM approaches with default and modified settings). All approaches started first from initial estimates set to the true values, and second using altered values. Results were examined through relative root mean squared error (RRMSE) of the estimates.

With true initial conditions, full completion rate was obtained with all approaches except FOCE in R. Runtimes were shortest with FOCE and LAPLACE, and longest with AGQ.

Under the rich design with true initial conditions, all approaches performed well except FOCE in R. When starting from altered initial conditions, AGQ, and then FOCE in NONMEM, LAPLACE in SAS, and SAEM in NONMEM and MONOLIX with tuned settings, consistently displayed lower RRMSE than the other approaches.

For standard dose-response models analyzed through mixed-effects models, differences could be identified in the performance of estimation methods available in current software.

INTRODUCTION

Non-linear mixed-effects models (NLMEM) were introduced to the biomedical field about 30 years ago [START_REF] Sheiner | Evaluation of methods for estimating population pharmacokinetic parameters. II. Biexponential model and experimental pharmacokinetic data[END_REF][START_REF] Sheiner | Evaluation of methods for estimating population pharmacokinetics parameters. I. Michaelis-Menten model: routine clinical pharmacokinetic data[END_REF][START_REF] Sheiner | Evaluation of methods for estimating population pharmacokinetic parameters. III. Monoexponential model: routine clinical pharmacokinetic data[END_REF] and have substantially improved the information learned from preclinical and clinical trials. Within drug development, NLMEM were initially used for pharmacokinetic (PK) analyses [START_REF] Ette | Population pharmacokinetics I: background, concepts, and models[END_REF], before being extended to pharmacokineticpharmacodynamic (PKPD) analyses [START_REF] Ette | Population Pharmacokinetics II: Estimation Methods[END_REF], along with dose-response analyses. On top of the structural mathematical model fit to PK or/and PD observations, the statistical model components enable the modeller to characterize results obtained in a set of individuals with the same parametric model and, in addition, to estimate the interindividual variability [START_REF] Sheiner | Estimation of population characteristics of pharmacokinetic parameters from routine clinical data[END_REF], and to quantify the unexplained variability [START_REF] Beal | Heteroscedastic Nonlinear Regression[END_REF].

The estimation of the fixed effect and random effect parameters involve complex estimation methods. Maximum likelihood estimation (MLE) approaches constitute a large family of methods commonly used in NLMEM analyses [START_REF] Dartois | Overview of model-building strategies in population PK/PD analyses: 2002-2004 literature survey[END_REF]. The non-linearity of the regression function in the random effects prevents a closed form solution to the integration over the random effects of the likelihood function [START_REF] Wang | Derivation of various NONMEM estimation methods[END_REF], thus several algorithms have been developed for MLE. Gaussian assumptions for the distribution of the random effects are common among MLE methods, and form the group of parametric approaches [START_REF] Pillai | Non-linear mixed effects modeling -from methodology and software development to driving implementation in drug development science[END_REF].

Along with methodological developments, different software have emerged, the most commonly used one [START_REF] Dartois | Overview of model-building strategies in population PK/PD analyses: 2002-2004 literature survey[END_REF] being NONMEM [START_REF] Beal | The NONMEM System[END_REF]. Estimation algorithms available were first restricted to First-Order (FO) and then First-Order Conditional Estimation (FOCE), which were subsequently implemented in Splus, R and WinNonMix. LAPLACE [START_REF] Wolfinger | Laplace's approximation for nonlinear mixed models[END_REF] then appeared in NONMEM, while SAS witnessed the addition of two macros MIXLIN and NLINMIX. A later procedure in SAS that represented a considerable improvement was NLMIXED, with FO and adaptive Gaussian quadrature (AGQ). Alternatives followed with stochastic expectation maximisation (EM) algorithms, and especially the SAEM algorithm [START_REF] Kuhn | Coupling a stochastic approximation version of EM with an MCMC procedure[END_REF] implemented in the MONOLIX [START_REF] Lavielle | MOdèles NOn LInéaires à effets miXtes)[END_REF] and the NONMEM [START_REF] Beal | NONMEM User's Guides[END_REF] software.

Whilst the estimation algorithms use different statistical methods, all aim at producing reliable estimates of the model parameters. The complexity of the model and the approximations embedded in the algorithm could potentially lead to poor estimation

METHODS

Statistical model

Let d = d 1 , …, d K be a set of ordered dose levels selected in a dose-response study and y ik be the response of subject i = 1, …, N to the dose d k . The dose-response is assumed to be adequately described by a function f such as:

ik k i ik y f(d , ) ε    ( 1
)
where i is the p dimensional vector of the model individual parameters for subject i and  ik is the measurement error.  ik given  i are assumed to be independent and normally distributed with a zero mean and a variance  ik 2 which can be additive ( ik 2 =  2 ) or proportional ( k 2 = f(d k  i   × 2 ). f is a function than can be nonlinear with respect to the parameters   i depend on the fixed effect p-dimensional vector  and the random effect q-dimensional vector  i in the following manner when considering an exponential model to ensure

positivity: i B η i μe     (2)
with the random effects following a Gaussian distribution with a zero mean and a variance matrix  of size (q×q), whose diagonal elements are variances   . The (p×q)-matrix B allows some components of not to have a random part. Also, the exponential random effect model ensures the positivity of the model parameter.

Finally, let define the vector of all the model parameters as = (',Vech()',) where the operator Vech(.) creates a column vector from the matrix by stacking its lower diagonal elements below one another.

Likelihood function

The log-likelihood L(y;) is the sum over the N subjects of the individual likelihoods, L(y i ;):

     N 1 i i i ) ;
(y L ) L(y; [START_REF] Sheiner | Evaluation of methods for estimating population pharmacokinetic parameters. III. Monoexponential model: routine clinical pharmacokinetic data[END_REF] where the individual log-likelihood L i (y i ; ) is defined as follows:

i i i i i i i
L (y ; ) log p(y , ; )d log p( y ; )p( ;

)d                         (4)
with p(y i | i ; ) the conditional density of the observations given the individual random effects, p( i ;) the density of the individual random effects, and p(y i ,  i ;) the likelihood of the 'complete' data which correspond to the observations plus the random effects,  i .

Estimation algorithms

Estimation methods are briefly described here. More details may be obtained in the original articles.

First-Order Conditional Estimation (FOCE)

As initially described by Lindstrom and Bates [START_REF] Lindstrom | Nonlinear mixed effects models for repeated measures data[END_REF], the algorithm approximates (4) by the log-likelihood of a linear mixed effect model. The  i and updated estimates of are obtained by minimizing a penalized nonlinear least square (PNLS) objective function using the current estimates of  and  Then, the model function f is linearized using a first-order Taylor expansion around the current estimates of  and the conditional mode of the  i so that (4) can be approximated by the log-likelihood of a linear mixed effect (LME) model to estimate  and . The maximization is realized through a hybrid approach starting with a moderate number of EM iterations before switching to Newton-Raphson iterations. The approach alternates between PNLS and LME until a convergence criterion is met. They implemented their method in the nlme function of the R software [START_REF]R: A Language and Environment for Statistical Computing[END_REF].

In the NONMEM software, the conditional modes of the  i are obtained by maximizing the empirical Bayes posterior density of  i , p( i |y i ;), using the current estimates of vector

:                )d ; )p( ; y p( ) ; )p( ; y p( ; y p i i i i i i i i (5)
Also, (4) is approximated by a second order Taylor expansion of the integrand (also called Laplacian approximation) around the  i ; however the Hessian is approximated by a function of the gradient vector to avoid the direct computation of second-order derivatives. For an additive residual error model, both the approximation by the linearization of the function f and the Laplacian approximation using an approximated Hessian have been shown to be equivalent asymptotically [START_REF] Wang | Derivation of various NONMEM estimation methods[END_REF]. However, this equivalence no longer holds in case of interaction between the  i and the  ik , as in the proportional error model. A derivative-free quasi-Newton type minimization algorithm is used.

Laplacian approximation (LAPLACE)

The principle of this algorithm is to approximate (4) by a second order Taylor expansion of the integrand around the conditional mode of the  i , which are obtained by maximizing the empirical Bayes posterior density of the  i using the current estimates of vector .

In the NLMIXED procedure of the SAS software [START_REF]9.1.2[END_REF], this algorithm is implemented as a special case of the adaptive Gaussian quadrature algorithm (see below) where only one abscissa is defined at the conditional modes of the  i with a corresponding weight equal to 1. Also, the  i are also obtained by maximizing p( i |y i ;) with a default dual quasi-Newton optimisation method.

Adaptive Gaussian Quadrature (AGQ)

The principle of this algorithm is to numerically compute (4) by a weighted average of p(y i | i ;) p( i ;) at predetermined abscissa for the random effects using a Gaussian kernel.

Pinheiro and Bates (31) suggested using standard Gauss-Hermite abscissa and weights [START_REF] Abramowitz | Handbook of Mathematical Functions: with Formulas, Graphs, and Mathematical Tables[END_REF], with the abscissa centred around the conditional mode of the  i and scaled by the Hessian matrix from the conditional mode estimation [START_REF] Pinheiro | Mixed-Effects Models in S and S-PLUS[END_REF]. The adaptive Gaussian approximation can be made arbitrarily accurate by increasing the number of abscissa.

Stochastic Approximation Expectation Maximization (SAEM)

SAEM is an extension of the EM algorithm where individual random effects are considered as missing data [START_REF] Lindstrom | Newton-Raphson and EM Algorithms for Linear Mixed-Effects Models for Repeated-Measures Data[END_REF]. It converges to maximum likelihood estimates by repeatedly alternating between the E and M steps. As the E step is often analytically intractable for nonlinear models, the E step in SAEM is replaced by a simulation step where the  i are drawn by running several iterations of a Hastings-Metropolis algorithm using three different kernels successively [START_REF] Lavielle | Estimation of population pharmacokinetic parameters of saquinavir in HIV patients with the MONOLIX software[END_REF]. Then the expectation of the complete log-likelihood Q() = E(log(p(y, η ))) is computed according to a stochastic approximation:

          m m 1 m m m 1 QQ γ log p(y, ; ) Q           (6)
where  m is a decreasing sequence of positive numbers over the m = 1, …, M algorithm iterations with  1 = 1. The SAEM algorithm has been shown to converge to a maximum (local or global) of the likelihood of the observations under very general conditions [START_REF] Delyon | Convergence of a stochastic approximation version of the EM algorithm[END_REF].

Simulation and estimation study

This simulation study consisted, for each studied scenario, of 100 stochastic simulated datasets generated in NONMEM and subsequently analysed with the different studied approaches (i.e. implementation of the estimation algorithms in the various software).

Simulations

Design

The dataset structure mimicked a clinical trial including 100 individuals and investigating four dose levels: 0, 100, 300 and 1000 mg. A continuous PD outcome was recorded for each individual following two simulation designs: (i) the rich design counted four observations per individual, one at each dose level, whereas (ii) in the sparse design each individual was randomly allocated to only two of the four dose levels.

Base model

A dose-response model based on a sigmoid E max function with a baseline (E 0 ) was constructed as in [START_REF] Beal | Heteroscedastic Nonlinear Regression[END_REF]. The Hill factor (γ) is responsible for the sigmoidicity, i.e. the degree of non-linearity of the function shape.

Gaussian random components with normal zero-mean distribution were assumed for all individual parameters except for γ. A correlation in the variances of the random effects for E max and ED 50 was assumed. The residual error model was assumed to be additive or proportional (see 2.1). Selected parameters values are reported Table 1.

Scenarios

Eight simulation scenarios (s = 8) were derived, exploring (i) the two previously described simulation designs: rich (R) and sparse (S), (ii) three values of γ: 1, 2, and 3, and (iii) two error models: additive (A) and proportional (P). They were referred to as: R1A, R2A, R3A, R1P, R2P, R3P, S3A, and S3P, and corresponded to eight sets of 100 simulated datasets to be analysed. Note that for the sparse design only sets with γ = 3, the most non-linear model, were evaluated.

Estimations

Initial conditions

The same model from which the simulated datasets were generated was used for estimation. Each dataset was analysed twice: (i) with true initial conditions, i.e. starting estimate values set to the original parameter values on which simulations were based, and (ii)

with altered initial conditions: γ set to 1, the other fixed effects to two fold of their true value, and random effects to low numbers (Table 1). This procedure explored the robustness of the approaches.

Software settings

Estimation algorithms were mostly utilised with the default settings with which they are available in the different studied software. Changes from these defaults were listed Concerning the convergence, it was stopped in the same manner as SAEM_NM, but instead of every 9999 iterations being submitted to the convergence test system, only every 25 were.. SAEM_MLX was run with setting the maximal number of iterations for the stochastic (K1≤500) and the cooling phase (K2≤200) using the following stopping rules: i) the stochastic phase is ended before K1 is reached if an iteration m is met where p(y, η m ; Ψ m ) < p(y, η m- 1K1 ; Ψ m-1K1 ) with lK1=100 and ii) the cooling phase is ended before K2 is reached if an iteration m is met where the variances of the parameters, computed over a window of lK2 iterations, is reduced by a factor rK2 compared to their values at the end of the stochastic phase, with lK2=50 and rK2=0.1. SAEM_MLX_tun was tuned in the way that it had a number of iterations for the stochastic phase, K1=500 (i.e. not using the stopping rule for this phase), and increased individual samples, nmc=5.

Hence nine approaches (a = 9) were explored through the estimation of the simulated datasets: FOCE_NM, FOCE_R, LAP_NM LAP_SAS, AGQ_SAS, SAEM_NM, SAEM_NM_tun, SAEM_MLX, and SAEM_MLX_tun.

Computer power

FOCE, LAPLACE and SAEM run in NONMEM 7.1.0 were assisted with PsN 3.2.5 [START_REF] Lindbom | PsN-Toolkit--a collection of computer intensive statistical methods for non-linear mixed effect modeling using NONMEM[END_REF] on a Linux cluster node of 3.59 GHz with a G77 Fortran compiler. Estimations with FOCE in R were done on a 2.49 GHz CPU as well as some with SAEM in MONOLIX (others were on a 1.83 GHz), assisted by a Matlab version R2009b. All SAS runs (LAPLACE and AGQ) were performed on a 2.66 GHz computer using SAS 9.2 for Windows.

Performance comparison

Completion rates

The proportion of completed estimations, i.e. the number K of the 100 analysed datasets that produced parameter estimates with each approach was reported. Other computations were executed with these Z sets of results; however when less than 50 of the runs completed, statistical measures were not produced. Z, thereafter expressed as a percentage, was therefore assessing the stability of the different approaches, whereas results were given only when K ≥ 50.

Runtimes

Runtimes were recorded as the CPU time needed to estimate each of the 100 copies of a simulated scenario. Then the average was calculated. A correction was done with the clock rate of the processor in the computer on which runs were performed as in [START_REF] Dartois | Overview of model-building strategies in population PK/PD analyses: 2002-2004 literature survey[END_REF]. Parallelization was not possible with the investigated approaches, so did not have to be accounted for.

K s ,a s ,a ,k .,a k1

1 N I C P U t C P U f K    (8)
where NI s,a is the calculated number of instructions in billions for scenario s with approach a, CPUt s,a,k the real time in seconds recorded on a CPU to perform the corresponding k th estimation, and CPUf .,a the frequency in GHz (equivalent to billion instructions per second) of the clock in the utilized CPU.

Accuracy and precision

Relative estimation errors (RER), relative bias (RBias), and root mean squared error (RMSE) were computed such that the accuracy and the precision of the estimation algorithms were evaluated for each of the 9 components (p) of the vector Ψ. The RER (%) are evaluated for each estimate and box-plot of RER(%) show both bias (mean) and imprecision (width).

The RBias (%) describes the deviation of the mean over the estimated parameters from their true value. The relative RMSE (RRMSE %) summarize both the bias and the variability in estimates. The Standardized RRMSE (%) was constructed for each parameter and each approach as the RRMSE divided by the lowest RRMSE value obtained across all approaches for that parameter in [START_REF] Beal | The NONMEM System[END_REF].

k* p ,s ,a p ,s ,. k p ,s ,a * p ,s ,. R E R ( ) 100             (9) K p ,s ,a k p ,s ,a k1 1 R Bias( ) R E R ( ) K      (10) K 2 p ,s ,a k p ,s ,a k1 1 R R M S E ( ) R E R ( ) K      (11)   p ,s ,a p ,s ,a a p ,s ,a R R M S E ( ) S tandardized R R M S E ( ) m in R R M S E ( )    (12)
where k p   is the estimated p component for the k th data set and Ψ p * the true value.

For each scenario and each approach, mean standardized RRMSE across the 9 components of  was computed as a global measure of the performance.

Computations were conducted in R 2.11.1.

RESULTS

Completion rates

100 % of the analyses started from true initial conditions completed with final estimates for all the approaches except FOCE_R (99, 62, 5, 69, 32, 2, 16, and 33 % for the R1A, R2A, R3A, R1P, R2P, R3P, S3A, and S3P scenarios, respectively) (Figure 2). The same simulated datasets estimated with altered starting values gave completion rates of the same order with FOCE_R (98, 76, 16, 68, 8, 3, 5, and 10 % for the R1A, R2A, R3A, R1P, R2P, R3P, S3A, and S3P scenarios, respectively), decreased ones with SAEM_NM (97, 91, 16, 74, 81, and 75 % for the R1A, R3A, R1P, R2P, R3P, and S3P scenarios, respectively) and SAEM_NM_tun (91 and 67 % for the R3A and S3A scenarios), and maximum completion (100 %) for all the other approaches. Therefore 133 sets of estimates were considered for further comparison statistics, 11 failing to meet the 50 % completion criterion.

Runtimes

Runtimes

Accuracy and precision

Boxplots of RER for ED 50 and ω 2 (ED 50 ) estimates are displayed on Figures 3a and3b as they often are the main parameters of interest in dose-response studies. Standardized RRMSE star-plots with 9 radii for each of the elements of  are represented in Figure 4; on a given radius, the closer to 1, the closer is the performance relative to the approach with the smallest RRMSE for the parameter of interest. For a global assessment across parameters, mean standardized RRMSE are illustrated in Figure 5.

True initial conditions

As displayed in Figure 3a, the parameter ED 50 was globally accurately estimated under true conditions, but presented a lower precision for scenarios with γ = 1. The highest and most consistent biases were observed with FOCE_R, on the few scenarios for which metrics were produced due to poor completion rates. ED 50 was better estimated with AGQ_SAS, LAP_NM and FOCE_NM on the sparse design than with the other tested approaches, which produced some bias, especially LAP_SAS (interquartile range excluding zero), and exhibited imprecision (wide interquartile range and longer whiskers), especially the SAEM approaches (except SAEM_NM). For the parameter ω 2 (ED 50 ) (Figure 3b), estimates were slightly more biased, but essentially more imprecise, especially with γ = 1, and the additive error model. For the sparse design, most approaches exhibited a bias, except the four SAEM approaches, which appeared to provide more accurate but less precise estimates than the other approaches.

SAEM_NM obtained the lowest RRMSEs whatever the scenario and parameter (values available in appendix); as illustrated in Figure 4, when γ > 1 and the error model was additive, all approaches but SAEM_NM estimated large E max , and when γ = 3 and the error model was proportional, all approaches but SAEM_NM estimated large ED 50 .

Globally on the rich design, as represented Figure 5, all approaches had a mean standardized RRMSE below 1.5 for most of the scenarios with the exception of FOCE_R.

Nevertheless, for scenario R3A, FOCE_NM and SAEM_MLX had it slightly above 1.5. On the sparse design, the LAPLACE methods, AGQ_SAS, and SAEM_NM had mean standardized RRMSEs below 1.5, whereas SAEM_MLX had it above 1.5 for both error models and SAEM_NM_tun and SAEM_MLX_tun for only the S3A scenario.

Altered initial conditions

On the rich design, most of the approaches estimated ED 50 similarly as when starting from true values, as illustrated in Figure 3a. However, the results of the SAEM approaches changed compared to the previous initial conditions case and sometimes drastically for the versions with the default settings, even failing to reach 50 % of completion for SAEM_NM with scenario R1P. On the sparse design, most of the methods obtained biased estimates, with the exceptions of AGQ_SAS, SAEM_NM, and FOCE_NM, which gave the distributions of 100 estimated ED 50 the most centred on the true value and tight. FOCE_R results could not be assessed, but the other approaches presented tailed distributions of estimated ED 50 , with quartiles not including the true value for LAP_SAS with both scenarios models and for SAEM_MLX with S3A. As shown in Figure 3b, the bias and imprecision in the ω 2 (ED 50 )

estimates were increased by starting from altered initial conditions particularly for SAEM_NM, whereas SAEM_MLX_tun yielded the boxplot most centred on zero.

It can be observed in Figure 4 that FOCE_NM and AGQ_SAS obtained standardized RRMSEs below 1.5 on most scenarios and parameters. When the sparse design was adopted the SAEM approaches and the LAPLACE approaches obtained standardized RRMSEs above 1.5 on most parameters, but for the proportional error model scenario they were below 1.5

with SAEM_NM_tun. FOCE_R estimated most parameters with poor standardized RRMSE, but especially γ and σ.

On Figure 5, FOCE_NM, and AGQ_SAS are shown to have lowest mean standardized RRMSE whatever the scenario, with LAP_SAS and SAEM_MLX_tun having mean standardized RRMSE below 1.5 for all but one scenario (S3P and S3A respectively).

FOCE_R obtained mean standardized RRMSE above 1.5 on all scenarios where its performance could be evaluated, whereas SAEM_NM, SAEM_MLX and LAP_NM also obtained elevated mean standardized RRMSE on at least half of the scenarios.

DISCUSSION

The present work provides a comparison in terms of speed, robustness, bias and precision of the most commonly used likelihood-based estimation approaches in nonlinear mixed effect modelling for the fitting of a dose-response model.

FOCE_R was shown to be the least robust approach with less than 50 % completion rate on 9 of the 16 combinations of scenarios and initial conditions settings investigated. All other approaches could be evaluated as they completed at least half of the data sets, with the exception of SAEM_NM in one situation. However the convergence criteria differed across estimation methods. In NONMEM, convergence of classical methods (FOCE and LAPLACE) is based only on the parameter estimation gradient, whereas it was set to be based on objective function, thetas, sigmas, and all omega elements for the SAEM methods. In MONOLIX, the automatic stopping rule for the stochastic phase is based on the complete log-likelihood. In SAS, convergence is primarily based on 6 key criteria, relating to the absolute and relative changes in the likelihood, gradients, and parameter values. The difficulty in defining convergence complicates these comparisons.

The convergence criteria used will affect runtimes, with less strict convergence criteria yielding shorter runtimes. However it is believed that the trends would remain the same, with the classical methods FOCE and LAPLACE being the fastest, and AGQ being the slowest.

AGQ slow runtimes were due to the high number of quadrature points chosen (9 quadrature points across 3 random effects imply 729 (9 3 ) likelihood evaluations for each subject at each iteration). Reducing this (e.g. to 3 quadrature points) would have significantly shortened the runtimes, and may have led to similar results (not inspected). Unsurprisingly, the estimation process speed was driven by the extent of the likelihood function simplification, with firstorder linearization-based algorithms achieving the shortest run times. Within each iteration, the SAEM approaches are faster than the Gaussian quadrature-based method because they sample the integrand rather than fully integrating it, but many more iterations are needed with SAEM than with AGQ. Increasing the number of chains to the SAEM algorithm was additionally time-consuming in MONOLIX, whereas SAEM_NM_tun was overall faster than SAEM_NM due to the number of iterations being decreased.

Globally, the approximation based on a linearization of the model, but for FOCE_R, gave good results for the fixed effects (relative biases typically less than 3 %) when starting from the true conditions, with ω 2 (ED 50 ) and Cov(E max ,ED 50 ) being least well estimated. As for their precision, it was decreasing in a similar extent using altered conditions and/or on a sparse design. The performance of adaptive Gaussian quadrature was high on all cases. The conclusions were less straightforward for the SAEM approaches. Indeed, SAEM_NM lacks a global search first step in order to refine the initial estimates; this could be appreciated with the results of the scenarios starting from altered values compared to SAEM_MLX. . However increasing the number of individual samples and linearly mu-referencing the parameters substantially improved the results. Mu-referencing appeared to yield more efficient behaviour of SAEM_NM_tun according to the implementation of the algorithm in NONMEM.

SAEM_MLX performance with altered initial conditions comes from the fact that it is coupled with a simulated annealing algorithm slowing up the decrease in variance estimates during the first iterations allowing escape from the local maxima of the likelihood and convergence to a neighbourhood of the global maximum. However, the more reduced the information is in the data, the more iterations and the more chains are needed to be provided in order to improve the convergence. Of note, on the S3A scenario with altered initial conditions, which is a particularly challenging combination of error model, Hill parameter value and design, the SAEM_NM_tun performance was improved using a user-supplied Omega shrinking algorithm for fixed effects parameters without interindividual variability instead of the default gradient process (results not shown). A similar Omega shrinking approach is implemented in MONOLIX.

One noticeable aspect about the investigated approaches is the possibility for user-defined options. The main advantage is the opportunity for the modeller to adapt the search to their specific problem. This makes it necessary for the user to be educated to the different alternatives, and their need might change during the model building, or worst, their nonutilization might influence the model selection. Nevertheless, an implementation always entails default settings, chosen by the developer and enlightened by common usage. Hence the same estimation algorithm existing in distinct software represents a dissimilar approach not only because of the implementation, but also because of the defaults settings. For that reason, explored approaches were primarily run with the options set to the defaults and secondarily with settings changed or tuned, when possible.

As estimation approaches in NLMEM require the user to provide initial values for the parameters to estimate, it was decided to assess the impact of these values on their performances. For the sake of simplicity, only two scenarios were considered, with initial guesses respectively correct and reasonably altered. The real case scenario would probably lie in between both situations as the user would first explore the data at hand, as well as use prior knowledge on the compound to come up with reasonable guesses. Of note, low initial values for the variances may provide less power to the EM-like algorithms for exploring the parameter search space, however in MONOLIX the simulated annealing inflates initial values for the variances.

Models investigated in the present study were dose-response models, based on the most commonly used structure in the field, a sigmoid E max . This model is fairly simple and contains a low number of parameters. The degree of nonlinearity is linked to the value of the Hill factor, which was varied across scenarios. Non-linearity is the major difficulty for ML estimation methods, for the reason mentioned earlier of no closed form solution for the integrand, whether the algorithm performs a linear approximation, a numerical integration or a stochastic approximation of the likelihood. Decreasing performance could hence be observed along the γ-increase with the additive error models, but not with the proportional error models, revealing other factors to take into account, such as the design. Models defined by ordinary differential equations represent also a challenge for estimation methods, and would perhaps result in conclusions of a different nature, but were not investigated in the present study.

Random effects are keys in the analysis of repeated data, allowing the modeller to quantify interindividual variability. The number of random effects that can be included in a model primarily depends on the amount of information generated under the chosen design, but also on the capacity of the algorithm to estimate them in addition to the fixed effects. The structure plays likewise a role, with considerations about the size of the variance-covariance matrices; therefore the studied structure included random effects on all parameters except one, plus one correlation.

Studies performing comparisons are bound to be limited by their tools. In the present work we used RMSE to sum-up information on both accuracy and precision which is a metric known to be sensitive to outliers. Yet, these choices provided us with the opportunity to present a readable comparison of 9 different estimation approaches across several combinations of true parameter values, error models and designs.

Drawbacks of FOCE_R experienced in this study had been described before [START_REF] Girard | A comparison of estimation methods in nonlinear mixed effects models using a blind analysis[END_REF].

Nevertheless, previously reported [START_REF] Plan | Performance in population models for count data, part I: maximum likelihood approximations[END_REF][START_REF] Jonsson | Estimating bias in population parameters for some models for repeated measures ordinal data using NONMEM and NLMIXED[END_REF] poor performance of LAP_NM for skewed distributions was not as evident in this study, where LAP_NM mean standardized RRMSE was low for all scenarios. However parameters on which performance was the poorest were variance of random effects, which was the case here also. These studies and additional ones [START_REF] Savic | Performance in population models for count data, part II: a new SAEM algorithm[END_REF][START_REF] Savic | Implementation and Evaluation of the SAEM Algorithm for Longitudinal Ordered Categorical Data with an Illustration in Pharmacokinetics-Pharmacodynamics[END_REF] showed estimates were improved with the use of AGQ_SAS or SAEM_MLX_tun;

these approaches gave good results here too. Another investigation [START_REF] Karlsson | Performance of the LAPLACE method in repeated time-to-event modeling[END_REF] highlighted that for cases with low information content LAP_NM had problems that disappeared when SAEM_NM was used. Again, this was only retrieved for variances of random effects, but was accordingly the case for the sparse design scenarios S3A and S3P. The impact of initial conditions had not been explored before, and this study showed the lack of robustness of some otherwise accurate estimation methods. Notwithstanding, it is important to realize that none of the NONMEM nor MONOLIX methods has been tested before, as the sofware have been updated since previous publications (from versions NONMEM VI and MONOLIX 2.4, respectively). Another comparison [START_REF] Chan | The use of the SAEM algorithm in MONOLIX software for estimation of population pharmacokineticpharmacodynamic-viral dynamics parameters of maraviroc in asymptomatic HIV subjects[END_REF] presenting EM methods as alternatives to gradientbased methods in terms of computation rates and runtimes was recently published (based on real data).

CONCLUSIONS

For standard dose-response models analyzed through mixed-effects models, differences could be identified in the performance of estimation methods available in current software.

Along with the exploration of different settings, designs and initial conditions, the strength of the present investigation resides in the inclusion of a high number of estimation methods and software.

TABLES:

FIGURESFigure 1 : 2 and 3 :

 123 FIGURES

Figure 2 :

 2 Figure 2: Percentage of completion and number of instructions (in billions) obtained with the 9 investigated approaches for the true initial conditions. The barchart represents the median, and the arrows link the minimum to the maximum value of the range.

Figure 3 :

 3 Figure 3:Relative estimation error (RER) for the parameter ED 50 (3a) and its variance (3b), for the 8 scenarios R1A, R2A, R3A, R1P, R2P, R3P, S3A, and S3P referring to 2 simulation designs (R for rich and S for sparse), 3 Hill factor values (1, 2, 3), and 2 residual error models (A for additive and P for proportional), with the estimation from true initial conditions and altered initial conditions. The boxplot represents the median (middle bar) and the interquartile range (box limits), with points for the mean (black) and the outliers (grey).

Figure 4 :

 4 Figure 4: Standardized RRMSE of the 9 population parameters for the 8 scenarios R1A, R2A, R3A, R1P, R2P, R3P, S3A, and S3P, and 2 initial conditions: true and altered. The following colour code was used: FOCE_NM = red, FOCE_R = pink, LAP_NM = orange, LAP_SAS = light green, AGQ_SAS = dark green, SAEM_NM = dark blue, SAEM_NM_tun = light blue, SAEM_MLX = dark violet, SAEM_MLX_tun = light violet.

Figure 5 :

 5 Figure 5:Strip chart of the mean standardized RRMSE obtained with each approach for the 8 scenarios R1A, R2A, R3A, R1P, R2P, R3P, S3A, and S3P, and 2 initial conditions: true and altered, on a semi-log scale. The colour code used is described in the Figure4legend. The star symbol (*) represents the S3A estimate from SAEM_NM_tun that is above 45 units. The dashed line is drawn at the value 1.5 used for description purposes in the results section.

Table 2

 2

	and reported below.
	FOCE and LAPLACE in NONMEM 7.1.0 (FOCE_NM and LAP_NM) had the maximum
	number of iterations set to the highest possible value as done in common practice, and the
	option INTERACTION was added for the scenarios with a proportional error. FOCE in
	R 2.9.1 (FOCE_R) was using the nlme routine. LAPLACE and AGQ in SAS 9.2 (LAP_SAS
	and AGQ_SAS) were adaptive Gaussian quadrature respectively corresponding to a number
	of quadrature points (QPOINTS) of 1 and 9. Other settings listed in table 2 were adapted from
	the defaults (FTOL=1E-15.7 XTOL=0 TECH=QUANEW EBSTEPS=50 EBSUBSTEPS=20
	EBSSFRAC=0.8 EBTOL=2.2E-12 INSTEP=1) in SAS. These settings were used previously
	(22) to improve robustness in the conditional modes calculations (the EB options) or to reduce
	the very high default convergence criteria (for FTOL and XTOL).
	SAEM presents a number of settings the user is invited to modify, that can follow different
	terminologies depending on the software: NONMEM 7.1.0/MONOLIX 3.1. These include the
	numbers NBURN/K1 and NITER/K2 of iterations in the stochastic ( k = 1) and the cooling
	(decreasing  k ) phases, respectively, as well as the number ISAMPLE/nmc of chains in the
	MCMC procedure. Stopping rules can also be defined for the two software for the stochastic
	phase, and also for the cooling phase in MONOLIX only. A simulated annealing version of
	SAEM during the first iterations can be set in NONMEM while it is automatically performed
	in MONOLIX. Moreover, i can be defined as the log-transform of a Gaussian random vector
	to meet with constraints of positivity, which corresponds to mu-referencing in NONMEM and
	the default in MONOLIX. In light of these possibilities, SAEM was run with each software

  initial conditions, respectively), never taking longer than 15 BI, very closely followed by FOCE_R and LAP_SAS. LAP_NM was displaying equivalently short runtimes for the additive error models (median NI = 10.2 BI and 11.2 BI, starting from true and altered initial conditions, respectively), which were doubled (median NI = 22.7 BI and 27.3 BI,

	starting from true and altered initial conditions, respectively) for the proportional error
	models, the design having no noticeable impact. SAEM approaches with default settings were
	systematically slower than FOCE and LAPLACE, but it was faster in MONOLIX (median NI
	= 43.2 BI and 52.6 BI, starting from true and altered initial conditions, respectively) than in
	NONMEM (median NI = 147.7 BI and 287.8 BI, starting from true and altered initial
	conditions, respectively), by around 3 folds when the initial conditions were true and 6 folds
	when they were altered. The tuned version of the approach, SAEM_MLX_tun, took around

expressed as number of instructions (NI) ranged from 4 to 1614 billion instructions (BI), and are displayed for estimations starting from true initial conditions in Figure 2. FOCE_NM was the fastest approach (median NI = 7.2 BI and 9.6 BI, starting from true and altered 2.5 times longer (median NI = 117.6 BI) than the non-tuned version, whereas SAEM_NM_tun (median NI =79.9 BI) was almost 3 times faster than SAEM_NM and 1.5 times faster than SAEM_MLX_tun; both had very similar runtimes between true and altered initial conditions. The NI reached with AGQ_SAS was high (median NI = 674.8 BI and 864.1 BI, starting from true and altered initial conditions, respectively); it was consistently the slowest.

Table 1 :

 1 True initial conditions are the parameter values used for the simulation of 8 scenarios constructed with 3 different Hill factor (γ) values and 2 different residual error models: additive (A) and proportional (P). True and altered initial conditions were used for the estimation of the simulated datasets.

	Parameters	E 0	E max	ED 50	γ	ω 2 (E 0 )	ω 2 (E max )	Cov(E max ,ED 50 )	ω 2 (ED 50 )		σ 2
	True initial conditions	5	30	500	1 2 3	0.090	0.490	0.245	0.490	A: 4	P: 0.010
	Altered initial conditions	10	60	1000	1 1 1	0.100	0.100	0.010	0.100	A: 1	P: 0.0625

Table 2 :

 2 Approach settings changed from defaults

	Approach	Algorithm	Software

ACKNOWLEDGMENTS E.L. Plan was supported by UCB Pharma, Belgium The authors would like to acknowledge Dr. R. Bauer (ICON Development Solutions, Ellicott City, USA), Dr. M. Lavielle (INRIA Saclay, University Paris-Sud, Orsay, France), and Dr. J. Pinheiro (Johnson & Johnson, Raritan, New Jersey, USA) for their valuable comments.


Julien Henry
email: julien.henry@imag.fr

AND David Monniaux
email: david.monniaux@imag.fr

Matthieu Moy
email: matthieu.moy@grenoble-inp.fr

SUCCINCT REPRESENTATIONS FOR ABSTRACT INTERPRETATION COMBINED ANALYSIS ALGORITHMS AND EXPERIMENTAL EVALUATION

interpretation techniques can be made more precise by distinguishing paths inside loops, at the expense of possibly exponential complexity. SMT-solving techniques and sparse representations of paths and sets of paths avoid this pitfall.

We improve previously proposed techniques for guided static analysis and the generation of disjunctive invariants by combining them with techniques for succinct representations of paths and symbolic representations for transitions based on static single assignment.

Because of the non-monotonicity of the results of abstract interpretation with widening operators, it is difficult to conclude that some abstraction is more precise than another based on theoretical local precision results. We thus conducted extensive comparisons between our new techniques and previous ones, on a variety of open-source packages.

Introduction

Static analysis by abstract interpretation is a fully automatic program analysis method. When applied to imperative programs, it computes an inductive invariant mapping each program location (or a subset thereof) to a set of states represented symbolically [START_REF] Cousot | Abstract interpretation frameworks[END_REF]. For instance, if we are only interested in scalar numerical program variables, such a set may be a convex polyhedron (the set of solutions of a system of linear inequalities) [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][START_REF] Bagnara | The Parma Polyhedra Library, version 0.9[END_REF][START_REF] Bagnara | The Parma Polyhedra Library: Toward a complete set of numerical abstractions for the analysis and verification of hardware and software systems[END_REF].

In such an analysis, information may flow forward or backward; forward program analysis computes super-sets of the states reachable from the initialization of the program, backward program analysis computes super-sets of the states coreachable from some property of interest (for instance, the violation of an assertion). In forward analysis, control-flow joins correspond to convex hulls if using convex polyhedra (more generally, they correspond to least upper bounds in a lattice); in backward analysis, it is control-flow splits that correspond to convex hulls.

It is a known limitation of program analysis by abstract interpretation that this convex hull, or more generally, least upper bound operation, may introduce states that cannot occur in the real program: for instance, the convex hull of the intervals [-2, -1] and [START_REF] Bagnara | The Parma Polyhedra Library, version 0.9[END_REF][START_REF] Bagnara | Precise widening operators for convex polyhedra[END_REF] is [-2, 2], strictly larger than the union of the two. Such introduction may prevent proving desired program properties, for instance = 0. The alternative is to keep the union symbolic (e.g. compute using [-2, -1] ∪ [START_REF] Bagnara | The Parma Polyhedra Library, version 0.9[END_REF][START_REF] Bagnara | Precise widening operators for convex polyhedra[END_REF]) and thus compute in the disjunctive completion of the lattice, but the number of terms in the union may grow exponentially with the number of successive tests in the program to analyze, not to mention difficulties for designing suitable widening operators for enforcing the convergence of fixpoint iterations [START_REF] Bagnara | The Parma Polyhedra Library, version 0.9[END_REF][START_REF] Bagnara | The Parma Polyhedra Library: Toward a complete set of numerical abstractions for the analysis and verification of hardware and software systems[END_REF][START_REF] Bagnara | Widening operators for powerset domains[END_REF]. The exponential growth of the number of terms in the union may be controlled by heuristics that judiciously apply least upper bound operations, as in the trace partitioning domain [START_REF] Rival | The trace partitioning abstract domain[END_REF] implemented in the Astrée analyzer [START_REF] Blanchet | A static analyzer for large safety-critical software[END_REF][START_REF] Cousot | The ASTR ÉE analyzer[END_REF].

Assuming we are interested in a loop-free program fragment, the above approach of keeping symbolic unions gives the same results as performing the analysis separately over every path in the fragment. A recent method for finding disjunctive loop invariants [START_REF] Gulwani | The reachability-bound problem[END_REF] is based on this idea: each path inside the loop body is considered separately. Two recent proposals use SMT-solving [START_REF] Kroening | Decision procedures[END_REF] as a decision procedure for the satisfiability of first-order arithmetic formulas in order to enumerate only paths that are needed for the progress of the analysis [START_REF] Gawlitza | Improving strategies via SMT solving[END_REF][START_REF] Monniaux | Using bounded model checking to focus fixpoint iterations[END_REF]. They can equivalently be seen as analyses over a multigraph of transitions between some distinguished control nodes. This multigraph has an exponential number of edges, but is never explicitly represented in memory; instead, this graph is implicitly or succinctly represented: its edges are enumerated as needed as solutions to SMT problems.

An additional claim in favor of the methods that distinguish paths inside the loop body [START_REF] Gulwani | The reachability-bound problem[END_REF][START_REF] Monniaux | Using bounded model checking to focus fixpoint iterations[END_REF] is that they tend to generate better invariants than methods that do not, by behaving better with respect to the widening operators [START_REF] Cousot | Abstract interpretation frameworks[END_REF] used for enforcing convergence when searching for loop invariants by Kleene iterations. A related technique, guided static analysis [START_REF] Gopan | Guided static analysis[END_REF], computes successive loop invariants for increasing subsets of the transitions taken into account, until all transitions are considered; again, the claim is that this approach avoids some gross over-approximation introduced by widenings.

All these methods improve the precision of the analysis by keeping the same abstract domain (say, convex polyhedra) but changing the operations applied and their ordering. An alternative is to change the abstract domain (e.g. octagons, convex polyhedra [START_REF] Miné | The octagon abstract domain[END_REF]), or the widening operator [START_REF] Bagnara | Precise widening operators for convex polyhedra[END_REF][START_REF] Halbwachs | Verification of real-time systems using linear relation analysis[END_REF]. This article makes the following contributions:

(1) We recast the guided static analysis technique from [START_REF] Gopan | Guided static analysis[END_REF] on the expanded multigraph from [START_REF] Monniaux | Using bounded model checking to focus fixpoint iterations[END_REF], considering entire paths instead of individual transitions, using SMT queries and binary decision diagrams (See §3).

(2) We improve the technique for obtaining disjunctive invariants from [START_REF] Gulwani | The reachability-bound problem[END_REF] by replacing the explicit exhaustive enumeration of paths by a sequence of SMT queries (See §4). (3) We implemented these techniques, in addition to "classical" iterations and the original guided static analysis, inside a prototype static analyzer. This tool uses the LLVM bitcode format [START_REF] Lattner | LLVM: A compilation framework for lifelong program analysis & transformation[END_REF][START_REF]LLVM Language Reference Manual[END_REF] as input, which can be produced by compilation from C, C++ and Fortran, enabling it to be run on many real-life programs. It uses the APRON library [START_REF] Jeannet | Apron: A library of numerical abstract domains for static analysis[END_REF], which supports a variety of abstract domains for numerical variables, from which we can choose with minimal changes to our analyzer. (4) We conducted extensive experiments with this tool, on real-life programs.

Bases

2.1. Static Analysis by Abstract Interpretation. Let X be the set of possible states of the program variables; for instance, if the program has 3 unbounded integer variables, then X = Z 3 . The set P(X) of subsets of X, partially ordered by inclusion, is the concrete domain. An abstract domain is a set X ♯ equipped with a partial order ⊑ (the associated strict order being ⊏); for instance, it can be the domain of convex polyhedra in Q 3 ordered by geometric inclusion. The concrete and abstract domains are connected by a monotone concretization function γ : X ♯ , ⊑ → (P(X), ⊆): an element x ♯ ∈ X ♯ represents a set γ(x ♯ ).

We also assume a join operator ⊔ : X ♯ × X ♯ → X ♯ , with infix notation; in practice, it is generally a least upper bound operation, but we only need it to satisfy γ(x ♯ ) ∪ γ(y ♯ ) ⊆ γ(x ♯ ⊔ y ♯ ) for all x ♯ , y ♯ .

Classically, one considers the control-flow graph of the program, with edges labeled with concrete transition relations (e.g. x ′ = x + 1 for an instruction x = x+1;), and attaches an abstract element to each control point. A concrete transition relation τ ⊆ X × X is replaced by an abstract forward abstract transformer τ ♯ :

X ♯ → X ♯ , such that ∀x ♯ ∈ X ♯ , x, x ′ ∈ X, x ∈ γ(x ♯ ) ∧ (x, x ′ ) ∈ τ =⇒ x ′ ∈ γ • τ ♯ (x ♯ ).
It is easy to see that if to any control point p ∈ P we attach an abstract element x ♯ p such that (i) for any p, γ(x ♯ p ) includes all initial states possible at control node p (ii) for any p, p ′ , τ ♯ p,p ′ (x ♯ p ) ⊑ x ♯ p ′ , noting τ p,p ′ the transition from p to p ′ , then (γ(x ♯ p )) p∈P form an inductive invariant : by induction, when the control point is p, the program state always lies in γ(x ♯ p ). Kleene iterations compute such an inductive invariant as the stationary limit, if it exists, of the following system: for each p, initialize x ♯ p such that γ(x ♯ p ) is a superset of the initial states at point p; then iterate the following: if

τ ♯ p,p ′ (x ♯ p ) ⊑ x ♯ p ′ , replace x ♯ p ′ by x ♯ p ′ ⊔ τ ♯ p,p ′ (x ♯ p )
. Such a stationary limit is bound to exist if X ♯ has no infinite ascending chain a 1 ⊏ a 2 ⊏ . . . ; this condition is however not met by domains such as intervals or convex polyhedra.

Widening-accelerated Kleene iterations proceed by replacing

x ♯ p ′ ⊔ τ ♯ p,p ′ (x ♯ p ) by x ♯ p ′ ▽(x ♯ p ′ ⊔τ ♯ p,p ′ (x ♯ p ))
where ▽ is a widening operator : for all x ♯ , y ♯ , γ(y ♯ ) ⊆ γ(x ♯ ▽ y ♯ ), and any sequence u

♯ 1 , u ♯ 2 , . . . of the form u ♯ n+1 = u ♯ n ▽ v ♯ n
, where v ♯ n is another sequence, become stationary. The stationary limit (x ♯ p ) p∈P , defines an inductive invariant (γ(x ♯ p )) p∈P . Note that this invariant is not, in general, the least one expressible in the abstract domain, and may depend on the iteration ordering (the successive choices p, p ′ ).

Once an inductive invariant γ((x ♯ p ) p∈P ) has been obtained, one can attempt decreasing or narrowing iterations to reduce it. In their simplest form, this just means running the following operation until a fixpoint or a maximal number of iterations are reached: for any p ′ , replace x ♯ p ′ by x ♯ p ′ ∩ p∈P τ ♯ p,p ′ (x ♯ p ) . The result also defines an inductive invariant. These decreasing iterations are indispensable to recover properties from guards (tests) in the program in most iteration settings; unfortunately, certain loops, particularly those involving identity (no-operation) transitions, may foil them: the iterations immediately reach a fixpoint and do not decrease further (see example in §2.3). Sections 2.4 and 2.5 describe techniques that work around this problem.

2.2. SMT-solving. Boolean satisfiability (SAT) is the canonical NP-complete problem: given a propositional formula (e.g. (a ∨ ¬b) ∧ (¬a ∨ b ∨ ¬c)), decide whether it is satisfiable -and, if so, output a satisfying assignment. Despite an exponential worst-case complexity, the DPLL algorithm [START_REF] Kroening | Decision procedures[END_REF][START_REF]Handbook of satisfiability[END_REF] solves many useful SAT problems in practice.

SAT was extended to satisfiability modulo theory (SMT): in addition to propositional literals, SMT formulas admit atoms from a theory. For instance, the theories of linear integer arithmetic (LIA) and linear real arithmetic (LRA) have atoms of the form a 1 x 1 +• • •+a n x n ⊲⊳ C where a 1 , . . . , a n , C are integer constants, x 1 , . . . , x n are variables (interpreted over Z for LIA and R or Q for LRA), and ⊲⊳ is a comparison operator =, =, <, ≤, >, ≥. Satisfiability for LIA and LRA is NP-complete, yet tools based on DPLL(T) approach [START_REF] Kroening | Decision procedures[END_REF][START_REF]Handbook of satisfiability[END_REF] solve many useful SMT problems in practice. All these tools provide a satisfying assignment if the problem is satisfiable.

2.3.

A Simple, Motivating Example. Consider the following program, adapted from [START_REF] Monniaux | Using bounded model checking to focus fixpoint iterations[END_REF], where input(a, b) stands for a nondeterministic input in [a, b] (the controlflow graph on the right depicts the loop body, s is the start node and e the end node):

1 void r a t e l i m i t e r ( ) { 2 i n t x o l d = 0 ; 3 while ( 1 ) { 4 i n t x = i n p u t ( -100000 , 100000) ;

5 i f ( x > x o l d +10) x = x o l d +10; 6 i f ( x < x o l d -10) x = x o l d -10; 7 x o l d = x ; 8 } } s e
This program implements a construct commonly found in control programs (in e.g. automotive or avionics): a rate or slope limiter.

The expected inductive invariant is x old ∈ [-100000, 100000], but classical abstract interpretation using intervals (or octagons or polyhedra) finds x old ∈ (-∞, +∞) [START_REF] Cousot | The ASTR ÉE analyzer[END_REF]. Let us briefly see why.

Widening iterations converge to x old ∈ (-∞, +∞); let us now see why decreasing iterations fail to recover the desired invariant. The x > x old+10 test at line 6, if taken, yields x old ∈ (-∞, 99990); followed by x = x old+10, we obtain x ∈ (-∞, 100000), and the same after union with the no-operation "else" branch. Line 7 yields x ∈ (-∞, +∞).

We could use "widening up to" or "widening with thresholds", propagating the "magic values" ±100000 associated to x into x old, but these syntactic approaches cannot directly cope with programs for which x ∈ [-100000, +100000] is itself obtained by analysis. The guided static analysis of [START_REF] Gopan | Guided static analysis[END_REF] does not perform better, and also obtains x old ∈ (-∞, +∞).

In contrast, let us distinguish all four possible execution paths through the tests at lines 6 and 7. The path through both "else" branches is infeasible; the program is thus equivalent to a program with 3 paths:

1 void r a t e l i m i t e r ( ) { 2 i n t x o l d = 0 ; 3 while ( 1 ) { 4 i n t x = i n p u t ( -100000 , 100000) ; 5 i f ( x > x o l d +10) x o l d = x o l d +10; 6 else i f ( x < x o l d -10) x o l d = x o l d -10; 7 else x o l d = x ; 8 } } s e
Classical interval analysis on this program yields x old ∈ [-100000, 100000]. We have transformed the program, manually pruning out infeasible paths; yet in general the resulting program could be exponentially larger than the first, even though not all feasible paths are needed to compute the invariant.

Following recent suggestions [START_REF] Gawlitza | Improving strategies via SMT solving[END_REF][START_REF] Monniaux | Using bounded model checking to focus fixpoint iterations[END_REF], we avoid this space explosion by keeping the second program implicit while simulating its analysis. This means we work on an implicitly represented transition multigraph ; it is succinctly represented by the transition graph of the first program. Our first contribution ( §3) is to recast the "guided analysis" from [START_REF] Gopan | Guided static analysis[END_REF] on such a succinct representation of the paths in lieu of the individual transitions. A similar explosion occurs in disjunctive

i n t x = 0 , y = 0 ; while ( 1 ) { i f ( x <= 50) y ++; else y--; i f ( y < 0 ) break ; x ++; } x y y ≤ x ∧ y ≤ 102 -x ∧ y ≥ 0. Figure 1.
Example program and its invariant: the piecewise linear, solid line is the strongest invariant, the grayed polyhedron is its convex hull. invariant generation, following [START_REF] Gulwani | The reachability-bound problem[END_REF]; our second contribution ( §4) applies our implicit representation to their method.

2.4. Guided Static Analysis. Guided static analysis was proposed by [START_REF] Gopan | Guided static analysis[END_REF] as an improvement over classical upward Kleene iterations with widening. Consider the program in Fig. 1, taken from [START_REF] Gopan | Guided static analysis[END_REF].

Classical iterations on the domain of convex polyhedra [START_REF] Cousot | Automatic discovery of linear restraints among variables of a program[END_REF][START_REF] Bagnara | Precise widening operators for convex polyhedra[END_REF] or octagons [START_REF] Miné | The octagon abstract domain[END_REF] start with x = 0 ∧ x = 0, then continue with x = y ∧ 0 ≤ x ≤ 1. The widening operator extrapolates from these two iterations and yields x = y ∧ x ≥ 0. From there, the "else" branch at line 4 may be taken; with further widening, 0 ≤ y ≤ x is obtained as a loop invariant, and thus the computed loop postcondition is

x ≥ 0 ∧ y = 0. Yet the strongest invariant is (0 ≤ x ≤ 51 ∧ y = x) ∨ (51 ≤ x ≤ 102 ∧ x + y = 102),
and its convex hull, a convex polyhedron (Fig. 1).

Intuitively, this disappointing result is obtained because widening extrapolates from the first iterations of the loop, but the loop has two different phases (x ≤ 50 and x > 50) with different behaviors, thus the extrapolation from the first phase is not valid for the second.

Gopan and Reps' idea is to analyze the first phase of the loop with a widening and narrowing sequence, and thus obtain 0 ≤ x ≤ 50 ∧ y = x, and then analyze the second phase, finally obtaining invariant (2.4); each phase is identified by the tests taken or not taken.

The analysis starts by identifying the tests taken and not taken during the first iteration of the loop, starting in the loop initialization. The branches not taken are pruned from the loop body, yielding:

while ( 1 ) { i f ( x <= 50) y ++; else break ; / * n o t taken i n phase 1 * / i f ( y < 0 ) break ;
x ++; } Analyzing this loop using widening and narrowing on convex polyhedra or octagons yields the loop invariant 0 ≤ x ≤ 51 ∧ y = x. Now, the transition at line 4 becomes feasible; and we analyze the full loop, starting iterations from 0 ≤ x ≤ 51 ∧ y = x, and obtain invariant (2.4) in Fig 1.

More generally, this analysis method considers an ascending sequence of subsets of the transitions in the loop body ; for each subset, an inductive invariant is computed for the program restricted to it. The starting subset consists in the transitions reachable in one step from the loop initialization. If for a given subset S in the sequence, no transitions outside S are reachable from the inductive invariant attached to S, then iterations stop; otherwise, add these transitions to S and iterate more. Termination ensues from the finiteness of the control-flow graph.

2.5. Path-focusing. Monniaux & Gonnord's path-focusing [START_REF] Monniaux | Using bounded model checking to focus fixpoint iterations[END_REF] technique distinguishes the different paths in the program in order to avoid loss of precision due to merge operations. Since the number of paths may be exponential, the technique keeps them implicit and computes them when needed using SMT-solving. The (accelerated) Kleene iterations ( §2.1) are computed over a reduced multigraph instead of the classical transition graph.

Let P be the set of control points in the transition graph, P W ⊆ P the set of widening points such that removing the points in P W gives an acyclic graph. One can choose a set P R such that P W ⊆ P R ⊆ P .

The set of paths is kept implicit by an SMT formula ρ expressing the semantics of the program, assuming that the transition semantics can be expressed within a decidable theory. For an easy construction of ρ, we also assume that the program is expressed in SSA form, meaning that each variable is only assigned once in the transition graph. This is not a restriction, since there exists standard algorithms that transform a program into an SSA format.

This formula contains Boolean reachability predicates b i for each control points

p i / ∈ P R , b s i and b d i for each p i ∈ P R , so that a path p i1 → p i2 → • • • → p in between two points p i1 , p in ∈ P R can easily be expressed as the conjunction b s i1 ∧ 2≤k<n b i k ∧ b d
in . The Boolean b s i is true when the path starts at point p i , whereas b d i is true when the path arrives at p i . In other words, we split the points in P R into a source point, with only outgoing transitions, and a destination point, with only incoming transitions, so that the resulting graph is acyclic and there are no paths going through control points in P R .

In order to find focus paths, we solve an SMT formula which is satisfiable when there exists a path starting at a point p i ∈ P R in a state included in the current invariant candidate X i , and arriving at a point p j ∈ P R in a state outside X j . In this case, we construct this path using the model and update X j . When p i = p j , meaning that the path is actually a self-loop, we can apply a widening/narrowing sequence, or even compute the transitive closure of the loop (or an approximation thereof, or its application to X i ) using abstract acceleration [START_REF] Gonnord | Combining widening and acceleration in linear relation analysis[END_REF].

We assume that we can encode the concrete semantics of the program into the SMT formula, or at least an abstraction thereof at least as precise as the one applied by the abstract interpreter (in simple terms: we want to avoid the case where the SMT solver exhibits a possible path, but the static analyzer realizes that this path is infeasible; this would lead to nontermination, because the SMT solver would exhibit the same path on the next iteration). A workaround would be to apply satisfiability modulo path programs [START_REF] Harris | Program analysis via satisfiability modulo path programs[END_REF]: from each path ruled infeasible by abstract interpretation, extract a blocking clause for the SAT solver underlying the SMT-solver.

Guided Analysis over the Paths

Guided static analysis, as proposed by [START_REF] Gopan | Guided static analysis[END_REF], applies to the transition graph of the program. We now present a new technique applying this analysis on the implicit multigraph from [START_REF] Monniaux | Using bounded model checking to focus fixpoint iterations[END_REF], thus avoiding control flow merges with unfeasible paths. In this section, we use the same notations as §2.5.

The combination of these two techniques aims at first discovering a precise inductive invariant for a subset of paths between two points in P R , by the mean of ascending and narrowing iterations. When an inductive invariant has been found, we add new feasible paths to the subset and compute an inductive invariant for this new subset, starting with the results from the previous analysis. In other words, our technique considers an ascending sequence of subsets of the paths between two points in P R . We iterate the operations until the whole program (i.e all the feasible paths) has been considered. The result will then be an inductive invariant of the entire program.

The ascending iteration applies path-focusing [START_REF] Monniaux | Using bounded model checking to focus fixpoint iterations[END_REF] to a subset of the multigraph. As [START_REF] Gopan | Guided static analysis[END_REF], we do some narrowing, to recover precision lost by widening, before computing and taking into account new feasible paths. Thus, our technique combines the advantages of Guided Static Analysis and Path-focusing.

Algorithm 1 performs Guided static analysis on the implicitly represented multigraph. I p denotes a set of initial states at program point p (thus ∅ for most p). The current working subset of paths, noted P and initially empty, is stored using a compact representation, such as binary decision diagrams. We also maintain two sets of control points:

• A ′ : points in P R that may be the starting points of new feasible paths.

• A : points in P R on which we apply the ascending iterations. When the abstract value of a control point p is updated, p is added to both A and A ′ .

Algorithm 1 Guided static analysis on implicit multigraph Select p i ∈ A ′ 10:

1: A ′ ← {p|P R /I p = ∅} 2: A ← ∅ 3: P ← ∅ //
A ′ ← A ′ \ {p i } 11:
ComputeNewPaths(p i ) // Update A, A ′ and P We distinguish three phases in the main loop of the analysis:

(1) We start finding a new relevant subset P of the graph. Either the previous iteration or the initialization led us to a state where there are no more paths in the previous subset P , starting at p i , that make the abstract values of the successors grow (otherwise, the SMT solver would not have answered "unsat "). Narrowing iterations preserve this property. However, there may exist such paths in the entire multigraph, that are not in P . This phase computes these paths and adds them to the subset. This phase is described in 3.2 and corresponds to lines in 8 to 12 in Algorithm 1. (2) Given a new subset P , we search for paths starting at point p i ∈ P R , such that these paths are in P , i.e are included in the working subgraph. Each time we find a path, we update the abstract value of the destination point of the path. This is the phase explained in 3.1, and corresponds to lines 14 to 18 in Algorithm 1. (3) We perform narrowing iterations the usual way (line 19 in algorithm 1) and reiterate from step 1 unless there are no more points to explore, i.e. A ′ = ∅. The order of steps is important: narrowing has to be performed before adding new paths, or spurious new paths would be added to P . Starting with the addition of new paths avoids doing the ascending iterations on an empty graph.

3.1.

Ascending Iterations by Path-focusing. For computing an inductive invariant over a subgraph, we use the Path-focusing algorithm from [START_REF] Monniaux | Using bounded model checking to focus fixpoint iterations[END_REF] with special treatment for self loops (line 17 in algorithm 1).

In order to find which path to focus on, we construct an SMT formula f (p i ), whose model when satisfiable is a path that starts in p i , goes to a successor p j ∈ P R of p i , such that the image of X i by the path transformation is not included in the current X j . Intuitively, such a path makes the abstract value X j grow, and thus is an interesting path to focus on. We loop until the formula becomes unsatisfiable, meaning that the analysis of p i is finished.

If we note Succ(i) the set of indices j such that p j ∈ P R is a successor of p i in the expanded multigraph, and X i the abstract value associated to p i :

f (p i ) = ρ ∧ b s i ∧ j∈P R j =i ¬b s j ∧ X i ∧ j∈Succ(i) (b d j ∧ ¬X j )
The difference with [START_REF] Monniaux | Using bounded model checking to focus fixpoint iterations[END_REF] is that we do not work on the entire transition graph but on a subset of it. Therefore we conjoin the formula f (p i ) with the actual set of working paths, noted P , expressed as a Boolean formula, where the Boolean variables are the reachability predicates of the control points. We can easily construct this formula from the binary decision diagram using dynamic programming, and avoiding an exponentially sized formula. In other words, we force the SMT solver to give us a path included in P . Each time the invariant candidate of a point p j has been updated, p j is inserted into A ′ since it may be the start of a new feasible paths.

Adding New Paths.

Our technique computes the fixpoint iterations on an ascending sequence of subgraphs, until the complete graph is reached. When the analysis of a subgraph is finished, meaning that the abstract values for each control point has converged to an inductive invariant for this subgraph, the next subgraph to work on has to be computed. This new subgraph contains all the paths from the previous one, and also new paths that become feasible regarding the current abstract values. The new paths in P are computed one after another, until no more path can make the invariant grow. This is line 11 in Algorithm 1, which corresponds to Algorithm 2. We also use SMT solving to discover these new paths, but we subtly change the SMT formula given to the SMT solver: we now try to find a path that is not yet in P , but is feasible and makes the invariant candidate of its destination grow. We thus check the satisfiability of the formula f ′ (p i ), where:

f ′ (p i ) = f (p i ) ∧ ¬P
X j is updated using an abstract union when the point p j is the target of a new path. This way, further SMT queries do not compute other paths with the same source and destination if it is not needed (because these new paths would not make X j grow, hence would not be returned by the SMT solver).

When a new path has been found, it is immediately added into P . We then have to add p i and p j into A (since we do not apply widening in this section) and p j into A ′ , since p j may be the starting point of a new feasible path. Compute the path e from the model 7:

X j ← X j ⊔ τ e (X i ) 8:

P ← P ∪ {e} 9: A ← A ∪ {p i } 10:
A ′ ← A ′ ∪ {p i } 11: end while 3.3. Termination. Termination of this algorithm is guaranteed, because: (1) the subset of paths P strictly increases at each loop iteration, and is bounded by the finite set of paths in the entire graph. ( 2) when computing new paths, we cunjunct our formula with ¬P , meaning that we obtain each possible path only once. The number of path is finite, so this computation always terminates. (3) the Pathfocusing iterations terminate because of the properties of widening.

3.4.

Example. We revise the rate limiter described in 2.3. In this example, Pathfocusing works well because all the paths starting at the loop header are actually self loops. In such a case, the technique performs a widening/narrowing sequence or accelerates the loop, thus leading to a precise invariant. However, in some cases, there also exists paths that are not self loops, in which case Path-focusing applies widening. This widening may induce unrecoverable loss of precision.

Suppose the main loop of the rate limiter contains a nested loop like:

void r a t e l i m i t e r ( ) { i n t x o l d = 0 ; while ( 1 ) { i n t x = i n p u t ( -100000 , 100000) ;

i f ( x > x o l d +10) x = x o l d +10; i f ( x < x o l d -10) x = x o l d -10; x o l d = x ; while ( w a i t ( ) ) {} } }
We choose P R as the set of loop headers of the function, plus the initial state. In this case, we have three elements in P R .

The main loop in the expanded multigraph has then 4 distinct paths going to the header of the nested loop.

Guided static analysis from [START_REF] Gopan | Guided static analysis[END_REF] yields, at line 3, x old ∈ (-∞, +∞). Pathfocusing [START_REF] Monniaux | Using bounded model checking to focus fixpoint iterations[END_REF] also finds x old ∈ (-∞, +∞). Now, let us see how our technique performs on this example.

Figure 2 shows the sequence of subset of paths during the analysis. The points in P R are noted p i , where i is the corresponding line in the code: for instance, p 3 corresponds to the header of the main loop.

(1) The starting subgraph is depicted on Figure 2 Step 1. At the beginning, this graph has no transitions. (2) We compute the new feasible paths that have to be added into the subgraph.

We first find the path from p 1 to p 3 and obtain at p 3 x old = 0.

p 1 p 3 p 8
Step 1

p 1 p 3 p 8
Step 2

x old ← 0

-10000 ≤ x ≤ 10000 x old -10 ≤ x x ≤ x old + 10/ x old ← x Figure 2.
Ascending sequence of subgraphs

The image of x old = 0 by the path that goes from p 3 to p 8 , and that goes through the else branch of each if-then-else, is -10 ≤ x old ≤ 10. This path is then added to our subgraph.

Moreover, there is no other path starting at p 3 whose image is not in -10 ≤ x old ≤ 10.

Finally, since the abstract value associated to p 8 is -10 ≤ x old ≤ 10, the path from p 8 to p 3 is feasible and is added into P . The final subgraph is depicted on Figure 2 Step 2.

(3) We then compute the ascending iterations by path-focusing. At the end of these iterations, we obtain -∞ ≤ x old ≤ +∞ for both p 3 and p 8 . ( 4) We now can apply narrowing iterations, and recover the precision lost by widening: we obtain -10000 ≤ x old ≤ 10000 at points p 3 and p 8 . ( 5) Finally, we compute the next subgraph. The SMT-solver does not find any new path that makes the abstract values grow, and the algorithm terminates. Our technique gives us the expected invariant x old ∈ [-10000, 10000]. Here, only 3 paths out of the 6 have been computed during the analysis. In practice, depending on the order the SMT-solver returns the paths, other feasible paths could have been added during the analysis.

In this example, we see that our technique actually combines best of Guided Static Analysis and Path Focusing.

Disjunctive Invariants

While many (most?) useful program invariants on numerical variables can be expressed as conjunctions of inequalities and congruences, it is sometimes necessary to introduce disjunctions. For instance, the loop for ( int i =0; i <n; i++) {...} has head invariant 0 ≤ i ≤ n ∨ (i = 0 ∧ n < 0). For this very simple example, a simple syntactic transformation of the control structure (into i =0; if ( i <n) do {...} while ( i <n)) is sufficient, but in more complex cases more advanced analyses are necessary [START_REF] Balakrishnan | Refining the control structure of loops using static analysis[END_REF][START_REF] Jeannet | Dynamic partitioning in linear relation analysis: Application to the verification of reactive systems[END_REF][START_REF] Sharma | Simplifying loop invariant generation using splitter predicates[END_REF][START_REF] Monniaux | Modular abstractions of reactive nodes using disjunctive invariants[END_REF]; in intuitive terms, they discover phases or modes in loops.

Gulwani & Zuleger [START_REF] Gulwani | The reachability-bound problem[END_REF] proposed a technique for computing disjunctive invariants, by distinguishing all the paths inside a loop. In this section, we propose to improve this technique by using SMT queries to find interesting paths, the objective being to avoid an explicit exhaustive enumeration of an exponential number of paths.

For each control point p i , we compute a disjunctive invariant 1≤j≤mi X i,j . We denote by n i the number of distinct paths starting at p i . To perform the analysis, one chooses an integer δ i ∈ [1, m i ], and a mapping function

σ i : [1, m i ] × [1, n i ] → [1, m i ].
The k-th path starting fom p i is denoted τ i,k . The image of the j-th disjunct X i,j by the path τ i,k is then joined with X i,σ(j,k) . Initially, the δ i -th abstract value contains the initial states of p i , and all other abstract values contain ∅.

For each control point p i ∈ P R , m i , δ i and σ i can be defined heuristically. For instance, one could define σ i so that σ i (j, k) only depends on the last transition of the path, or else construct it dynamically during the analysis.

Our method improves this technique in two ways :

• Instead of enumerating the whole set of paths, we keep them implicit and compute them only when needed. • At each loop iteration of the original algorithm [START_REF] Gulwani | The reachability-bound problem[END_REF], an image by each path inside the loop is computed for each disjunct of the invariant candidate. Yet, many of these images may be redundant: for instance, if our invariant candidate is (0 ≤ x ≤ 10 ∧ 0 ≤ y ≤ 1000) ∨ (x < -10 ∧ y < -10), then there is no point enumerating paths whose image is included in this invariant candidate. In our approach, we compute such an image only if it makes the resulting abstract value grow. Our improvement consists in a modification of the SMT formula we solve in 3. We introduce in this formula Boolean variables {d j , 1 ≤ j ≤ m}, so that we can easily find in the model which abstract value of the disjunction of the source point has to be chosen to make the invariant of the destination grow. The resulting formula that is given to the SMT solver is defined by g(p i ). When the formula is satisfiable, we know that the index j of the starting disjunct that has to be chosen is the one for which the associate Boolean value d j is true in the model. Then, we can easily compute the value of σ i (j, k), thus know the index of the disjunct to join with.

g(p i ) = ρ ∧ b s i ∧ j∈P R j =i ¬b s j ∧ 1≤k≤mi (d k ∧ X i,k ∧ l =k ¬d l ) ∧ j∈Succ(i) (b d j ∧ 1≤k≤mi (¬X j,k ))
In our algorithm, the initialization of the abstract values slightly differs from algorithm 1 line 5, since we now have to initialize each disjunct. Instead of Line 5, we initialize X i,k with ⊥ for all k ∈ {1, .., m i } \ {δ i }, and X i,δi with ← I pi .

Furthermore, the Path-focused algorithm (line 17 from algorithm 1) is enhanced to deal with disjunctive invariants, and is detailed in algorithm 3.

The Update function can classically assign X i,σi(j,k) ▽(X i,σi(j,k) ⊔ τ i,k (X i,j )) to X i,σi(j,k) , or can integrate the special treatment for self loops proposed by [START_REF] Monniaux | Using bounded model checking to focus fixpoint iterations[END_REF], with widening/narrowing sequence or acceleration. Compute the path τ i,k from res Update(X i,σi (j,k) ) 9: end while We experimented with a heuristic of dynamic construction of the σ i functions, adapted from [START_REF] Gulwani | The reachability-bound problem[END_REF]. For each control point p i ∈ P R , we start with one single disjunct (m i = 1) and define δ i = 1. M denotes an upper bound on the number of disjuncts per control point.

The σ i functions take as parameters the index of the starting abstract value, and the path we focus on. Since we dynamically construct these functions during the analysis, we store their already computed image into a compact representation, such as Algebraic Decision Diagrams. σ i (j, k) is then constructed on the fly only when needed, and computed only once. When the value of σ i (j, k) is required but undefined, we first compute the image of the abstract value X i,j by the path indexed by k, and try to find an existing disjunct of index j ′ so that the least upper bound of the two abstract values is exactly their union (using SMT-solving). If such an index exists, then we set σ i (j, k) = j ′ . Otherwise:

• if m i < M , we increase m i by 1 and define σ

i (j, k) = m i • if m i = M , we define σ i (j, k) = M
The main difference with the original algorithm [START_REF] Gulwani | The reachability-bound problem[END_REF] is that we construct σ i (j, k) using SMT queries instead of enumerating a possibly exponential number of paths to find a solution.

Implementation and Experimental Comparisons

We have implemented our proposed solutions inside a prototype of intraprocedural static analyzer called PAGAI, as well as the classical abstract interpretation algorithm, and the state-of-the-art techniques Path Focusing [START_REF] Monniaux | Using bounded model checking to focus fixpoint iterations[END_REF] and Guided Static Analysis [START_REF] Gopan | Guided static analysis[END_REF]. It is available online at https://forge.imag.fr/projects/pagai/. The implementation is documented in [START_REF] Henry | Static analysis by path focusing[END_REF].

PAGAI operates over LLVM bitcode [START_REF]LLVM Language Reference Manual[END_REF][START_REF] Lattner | LLVM: A compilation framework for lifelong program analysis & transformation[END_REF], which is a target for several compilers, most notably Clang (supporting C and C++) and llvm-gcc (supporting C, C++, Fortran and Ada). Abstract domains are provided by the APRON library [START_REF] Jeannet | Apron: A library of numerical abstract domains for static analysis[END_REF], and include convex polyhedra (from the builtin Polka "PK" library), octagons, intervals, and linear congruences. For SMT-solving, our analyzer uses Yices [START_REF] Dutertre | A fast linear-arithmetic solver for DPLL(T)[END_REF] or Microsoft Z3 [START_REF] Mendonça De Moura | Z3: An efficient SMT solver[END_REF].

PAGAI currently neither models the memory heap nor performs interprocedural analysis. Instead, LLVM optimization phases are applied prior to analysis, in order to inline non-recursive function calls and lift certain memory accesses to operations on explicit numerical variables (e.g. y=t [0] * t [0]; preceded by t [0]=x; without any aliased write in between is replaced by y=x * x;). The remaining memory reads are considered as indeterminates, and memory writes are ignored; this is a sound abstraction.

We conducted extensive experiments on real-life programs in order to compare the different techniques, mostly on open-source projects (Fig. 3) written in C, C++ and Fortran. These results confirm that our combined technique improve the analysis in comparison with the two techniques taken individually, at a reasonable cost. The extension with disjunctive invariants increases precision in many cases, but with higher cost in terms of execution time.

Conclusion and Future Prospects

Roughly, an analysis by abstract interpretation is defined by the choice of an iteration strategy and an abstract domain. In this article, we demonstrated that changes in the iteration algorithm can significantly improve precision, sometimes while improving analysis times.

A common criticism of analysis techniques based on SMT-solving is that they do not scale up. Yet, our experiments show that, for numerical properties, they scale up to the size of typical functions and loops. It is however quite certain that, naively applied, they cannot scale to the kind of programs targeted by e.g. the Astrée tool, that is, a dozens or hundreds of thousands of lines of code in a single loop operating over similar numbers of remanent variables. Actually, for such applications, only (quasi-)linear algorithms scale up, and "cheap" abstract domains such as octagons (O(n 3 ) where n is the number of variables) are not applied to the full variable set, but to restricted subsets thereof. It thus seems reasonable that techniques such as considering "packs" of related variables, slicing, etc. may similarly help SMT We compared the precision of different techniques and abstract domains by comparing the invariants for the inclusion ordering. A better metric is perhaps to take a client analysis -such as the detection of overflows and array bound violations -and compare the rates of alarms.

We focused on numerical properties, because they are supported by easily available abstract libraries. Yet, in most programs, properties of data structures are important for proving interesting properties. Further investigations are needed not only on good abstractions for pointers (many are already known) but also on their conversion to SMT problems.

Paths in the current subset 4 : 5 :

 45 for all p i ∈ P R do X i ← I pi 6: end for 7: while A ′ = ∅ do 8:while A ′ = ∅ do 9:

15 :

 15 Select p i ∈ A 16: A ← A \ {p i } 17:PathFocusing(p i ) // Update A and A ′ end while 21: return {X i , i ∈ P R }

Algorithm 2 1 :

 21 ComputeNewPaths while true do 2:res ← SmtSolve [f ′ (p i )]

Algorithm 3

 3 Disjunctive invariant computation with implicit paths 1: while true do 2: res ← SmtSolve [g(p i )]

7 :

 7 Take j ∈ {l|d l = true} 8:

Figure 3 .

 3 Figure 3. Comparison of the abstract values obtained on several open-source projects. The table shows their respective number of lines of code, number of control points in P R , and execution time on various techniques. Techniques are classical abstract interpretation (S), Guided Static Analysis (G), Path-focused technique (PF), our combined technique (G+PF), and its version with disjunctive invariants (DIS). The bars (resp. ) gives the percentage of invariants stronger (more precise; smaller with respect to inclusion) with the left-side (resp. right-side) technique, and "uncomparable" gives the percentage of invariants that are uncomparable, i.e neither greater nor smaller; the code points where both invariants are equal make up the remaining percentage.

Table 1 .

 1 -based techniques to scale to global analyses. Execution times for various techniques

		Size			Execution time (seconds)
	Name	kLOC |P R |	S	G	PF G+PF DIS
	a2ps-4.14	55	2012	23	74	34	115	162
	gawk-4.0.0	59	902	15	46	12	40	50
	gnuchess-6.0.0	38	1222	50	220	81	312	351
	gnugo-3.8	83	2801	77	159	92	766	1493
	grep-2.9	35	820	41	85	22	65	122
	gzip-1.4	27	494	22	268	91	303	230
	lapack-3.3.1	954 16422 294 3740 3773	8159 10351
	make-3.82	34	993	67	108	53	109	257
	tar-1.26	73	1712	37	218	115	253	396

This work was partially funded by ANR project "ASOPT".