Julien Worms

Going beyond mean and median programs performances

Keywords: programs performances modelling, gaussians mixture, programs performances variability, parametric statistics

.

I. INTRODUCTION AND MOTIVATION

When you read books or articles on computer architectures, performance analysis, operating systems (OS) or compilation, some people still think that the execution time of a program P on a fixed machine M with fixed data input I is stable around a value, that can be noted as a single number ExecutionTime(P,I,M). This situation was true in old computers, but nowadays nobody really observes constant execution times, except in rare situations: ideal execution environment, special processor architectures designed for performance stability, bare mode execution, etc. Even SPEC organisation resumes the execution times of a program as a single number (computed as the sample median of a set of observed values).

The consequence of performance variability is that the reported constant numbers of program performances in the literature are not reproducible, and it becomes more and more difficult to decide about the most effective program version or OS configuration. As in car industry, a random car driver can hardly ever reproduce the gas consumption reported by the vendor.

There are many factors that explain why the observed programs performances are instable: technological factors (electronics, variable clock frequency, input/output peripherals, etc), micro-architectures (hardware prefetching, out of order execution, speculation, memory hierarchy), thread and process competition on shared ressources or data, operating systems services and actions, algorithmic factors (some application are designed to not behave exactly from one run to another), etc. If the user has complete control of its execution environment, he may be able to stabilise performances by configuring the BIOS (to deactivate some hardware features), by configuring the OS, by reducing the workload of the shared machine, etc. Unfortunately regular users cannot get such privileged rights, they are just allowed to use some computing power. Such users observe that their programs have unstable performances. These users may rely on rigorous statistical methods to decide which code version they should use for their production computation.

The instability of programs performances makes difficult to check whether a code version is better than another. So people usually make standard statistics based on mean or median values. That is, people resume all the observed performance numbers by a single value, which is not adequate in practice. The reason is that the programs performances are generally not of simple gaussian type in practice, as demonstrated in [START_REF] Worms | Parametric and Non-Parametric Statistics for Program Performance Analysis and Comparison[END_REF]. So what is the consequence if the performances of a program do not follow a gaussian distribution ? In theory, it means that some statistical tests designed for gaussians (such as the Student t-test) would compute a wrong risk level. Likewise, the well known formula that computes the confidence interval of an average would correspond to a wrong confidence level. Of course, if the performances sample is large enough, the statistical error must be asymptotically bounded (but nobody is able to define how large it should be).

In practice, the observed performances are often multimodal data distributions (see Sect. II). This suggests us to use gaussian mixtures to model the density distributions. Based on this modeling, we are able to do rigorous and flexible parametric statistics. We define and compute new performance metrics that go beyond mean and median values. These metrics will provide additional decision criteria to the user to analyse and to compare programs performances.

Now, what does the word parametric mean for a statistical model ? This is fairly simple: if X = (x 1 , . . . , x n ) denotes some execution times data (n repetitions of program executions), we suppose that they are independent realisations of some probability distribution having a probability density function (p.d.f.) f X . The approach will be parametric if we suppose that this p.d.f. belongs to some family of distributions f θ where θ is a vector parameter composed of d real-valued sub-parameters, which exact values are unknown: d is the dimension of the model. Since the set of all possible distributions is infinite dimensional, we thus reduced the problem of estimating f X to the problem of estimating a finite-dimensional parameter θ: this is the main purpose of the parametric approach.

Let us also mention that, despite the appeal of the nonparametric approach, the parametric approach is still very heavily used in the scientific activities involving statistics thanks to numerous reasons:

1) The parameters of a parametric model can often be interpreted (which helps understanding and summarising the data); 2) Some computations cannot be performed, or some mathematical results cannot be proved, without parametrically specifying some parts of the model (for example, the perturbation part); 3) For really large datasets, it is often proved that the parametric techniques are indeed more accurate than their non-parametric alternatives; 4) For high dimensional data, non-parametric techniques may suffer from a lack of accuracy which is called the "curse of dimensionality", and which is out of the scope of our present research work. It is important to notice that our statistical methods treat any kind of continuous performance data (execution times, energy consumption, network traffic, memory bandwidth, etc). We make experiments and demonstrations with programs execution times because this sort of performance is the most common and easiest to collect. Any other kind of continuous performance can be analysed using the parametric statistics we present in this article.

Our article is organised as follows. Sect. II presents our method of modelling the performance data using gaussian mixtures. Sect. III formally defines and computes four new performance metrics, and explain why they are useful in practice. Sect. IV resumes our experimental study. Related work and discussion are given in Sect. V, then we conclude.

II. MODELING PROGRAMS PERFORMANCES WITH GAUSSIAN MIXTURES

The central idea of our work in this direction stems from the following remark: in our own experience in parallel and sequential application performance analysis, the variability of such data is not necessarily of the "deviation around a single mean value" kind, it often exhibits a clear clustering pattern: in other words, the observed performances often vary from each other by clustering around two or more central values. Examples of such observed performances are given in pages 16 and 17 in [START_REF] Worms | Parametric and Non-Parametric Statistics for Program Performance Analysis and Comparison[END_REF]. Therefore, we cannot choose for our modelling classical families of distributions such as the Gaussian, Exponential, Gamma or Weibull family, which are by nature monotonic or unimodal.

We propose to model the performances by mixtures of gaussian distributions: this family of distributions has demonstrated to be an essential tool in many areas of scientific activities for many years now (biology, engineering, astronomy, among many others), particularly in the image analysis and pattern recognition fields. Mixtures of gaussian distributions is a highly flexible family of distributions which can fit a great variety of data: it is not only naturally adequate for modelling data which exhibits clusters, but it can also handle the problem of possible skewness in the data, despite the symmetry of the gaussian components of the mixtures.

A. Definition of the gaussian mixtures family

Remind that we consider data sample X = (x 1 , . . . , x n ) which are independent realisations of a probability density function (p.d.f.) f X . We say that the data are issued from a finite mixture of gaussian distributions (or simply a gaussian mixture, which we will abbreviate by GM from now on) if f X is equal to some p.d.f. g θ,K (parametrized by θ and K described below) of the form:

g θ,K (x) = π 1 ϕ(x; µ 1 ; σ 1 ) + . . . + π K ϕ(x; µ K ; σ K ) = K k=1 π k ϕ(x; µ k ; σ k ) (1) where:
-π 1 , . . . , π K ∈ R + are the mixture weights, which are positive and sum to 1; -µ 1 , . . . , µ K and σ 1 , . . . , σ K are the mean values and standard deviations of the mixture individual components; -ϕ( • ; µ k ; σ k ) denotes the p.d.f. of the gaussian/normal distribution N (µ k , σ k ). Also, for the sequel of this article, we note Φ the cumulative distribution function (c.d.f.) of ϕ( • ; µ k ; σ k ); -K ∈ N is the number of components in this mixture; -θ is a vector gathering all the parameters (except K) in a single notation,

θ = (π 1 , . . . , π K ; µ 1 , . . . , µ K ; σ 1 , . . . , σ K )
We will denote by F GM the set of all mixtures of gaussian distributions, and say that X is GM-distributed if its cumulative distribution function (c.d.f.) F X belongs to F GM , which means there exists some parameters K and θ

= (π 1 , . . . , π K ; µ 1 , . . . , µ K ; σ 1 , . . . , σ K ) such that ∀x ∈ R, F X (x) equals F θ,K (x) = K k=1 π k Φ(x; µ k ; σ k ) (2)
which is (of course) itself equivalent to f X being equal to the density g θ,K defined in Equ. 1. Naturally, the more components the mixture has, the more flexible the shape of the distribution can be (but this has a cost: the model has more parameters to be estimated).

The next section presents the method for building a GM model based on a data sample. This is called clustering in the literature.

B. Clustering method

Our aim is then, for the moment, to estimate the parameters of the distribution described in Equ. 1, from which we assume our data are issued: the parameters are the clusters weights (π k ) k=1..K , the clusters means (µ k ) k=1..K and the clusters standard deviations (σ k ) k=1..K , and they are gathered in one single notation, θ (which is 3K-dimensional). The formal details of our clustering method are given in [START_REF] Worms | Parametric and Non-Parametric Statistics for Program Performance Analysis and Comparison[END_REF], this section is a synthesis.

1) Estimation of the parameters of the mixture for fixed K: How then is θ estimated? Remind that the dimension of θ depends on the value of K, which is fixed for the moment. We naturally adopt a parametric approach, and intend to compute the maximum likelihood estimator θ of θ, which is defined as the (global) maximiser of the log-likelihood function. This function, given the observations x 1 , . . . , x n and Equ. 1, is defined by:

L(θ) = n i=1 log g θ,K (x i ) = K k=1 π k ϕ(x i ; µ k ; σ k )
Maximising it consists in calculating the various partial derivatives of L with respect to the different parameters π k , µ k , σ k (k = 1, . . . , K), and equalising them to 0 to obtain the so-called score equations. It is rather clear that any attempt to directly solve these score equations will turn out to be an unsolvable problem, due to the presence of a log of a sum.

This major obstacle was overcome thanks to the approach synthesised in the celebrated paper [START_REF] Dempster | Maximum likelihood from incomplete data via the em algorithm[END_REF]: the EM algorithm. The acronym EM means a succession of E-steps (E for Expectation) and M-steps (M for Maximisation), which eventually lead to obtaining (numerically) the value of the maximum likelihood estimator θ. We do not detail this algorithm here, we sketch/vulgarise/explain it in [START_REF] Worms | Parametric and Non-Parametric Statistics for Program Performance Analysis and Comparison[END_REF].

2) Determination of the number K of components of the mixture: Let us explain how the number K of components of the gaussian mixture model can be chosen. In practice, several candidate values are considered for K, and one of them, noted K, is chosen so that the corresponding GM model best fits the data at hand. The determination of K is nearly the most important issue in clustering analysis, and in this work we adopt a simple and widespread strategy: using the BIC criterion (BIC stands for Bayesian Information Criterion).

The principle of the BIC criterion for determining K is the following. If, for a given K ≥ 1, we note L K the maximum value of the log likelihood for the model with K components, then it should be easy to conceive that the greater K is, the greater the value L K will be: indeed, for instance, if you consider the model with K + 1 components which best fits your data, then it will certainly better fit your data than the best model having K components (maybe not far better, but better all the same). Therefore, choosing K which maximises L K would not work. The likelihood value needs to be penalised by a value which grows with K, in order to counterbalance the fitting gain that more complex models yield. That is the idea of the so-called information criterions, for instance the BIC criterion: to choose the value of K that minimises the value BIC(K) = -2L K + K log(n) (this value of the penalisation term K log(n) has theoretical justifications, which will not be detailed here). Therefore, if BIC( K) = max K≥1 BIC(K), then the model with K components will be a tradeoff between good data fitting and reasonable complexity. Note that, in practice, the maximum is chosen among a finite number of candidate values, for instance 1 ≤ K ≤ K max (where K max does not exceed 10 in general).

This section presented the skeleton of our clustering method based on EM algorithm. Now, we are able to build a GM model for any sample X. The next section provides two simple examples.

3) Simple examples: Let us consider two well known SPEC benckmarks which are named galgel and apsi. Running these benchmarks multiple times on the same low overhead machine results in variable execution times, as illustrated by histograms in Fig. 1. The histograms show the frequency of observed execution times (sample of 35 runs). After clustering, the approximate theoretical models based on GM are illustrated with the curves in Fig. 1. As can be seen, the execution times of these benchmarks exhibit multi-modal behaviour.

However, how can we check if this computed model fit well the data or not? This is the purpose of a goodnessof-fitting test. We developed such test in details in [START_REF] Worms | Parametric and Non-Parametric Statistics for Program Performance Analysis and Comparison[END_REF], and demonstrated its robustness. We will discuss it later in Sect. V. We draw in Fig . 2 the empirical CDF (step function) versus theoretical CDF (cuve). For both examples, the fitting is very good.

The next section presents new performance metrics for analysing and comparing programs performances, based on parametric statistics.

III. NEW PROGRAM PERFORMANCE METRICS

In the literature, people are mainly focused on the average or median program performance. However in practice, the average or the median value may not be the most interesting summary measure that reflects the program performance, or may not be the best performance metric to make a decision about selecting the most suitable program version.

For instance, consider the situation of a very long running application that a user executes very few times. The user has the choice between many code versions, which one should he select? If he bases his choice only on the expected average or median execution time, he may be disappointed if he executes his application very few times. If an application is rarely executed, the mean or median performances are not felt, especially if the performance distribution is multimodal. So, additional performance metrics can help him making better selection.

Clustering of apsi performances

Also, consider the situation where a user wants to know if the performances of his application are stable or not. Which metric can he use? The well known variance is a metric that measures how data is spread out around the average only: knowing how to interpret this metric is not so widespread among the practicians, and can be misleading when the data distribution is multi modal, because the variance is a measure of dispersion around a single value, the average. Therefore it cannot be the unique metric used for evaluating performance stability, additional metrics can be introduced and used, as we will see later.

This section provides new performance metrics that help the user to select a good program version based on performance analysis. Let X and Y two random variables, representing the performances of two code versions. Let X be a sample of X and Y be a sample of Y , meaning that X = (x 1 , . . . , x n ) and Y = (y 1 , . . . , y m ), with n and m denoting the respective sample sizes.

A. The metric I 1 : the mean difference

We may be interested in quantifying the average difference between the performances of two code versions. That is, we may be interested in computing the expected value E [ |X -Y | ]. This defines our first performance metric as

I 1 = E [ |X -Y | ].
Our parametric estimation, noted I 1 , assumes that both X and Y are modelled with gaussian mixture distributions. This means that the underlying p.d.f. f X and f Y of X and Y equal weighted combinations of gaussian p.d.f.

f X (x) = K i=1 π i ϕ(x; µ i ; σ i ) and f Y (x) = K j=1 π j ϕ(x; µ j ; σ j )
where we recall here that ϕ(x; µ; σ) denotes the p.d.f. of the gaussian distribution N (µ, σ), and K and K are the respective number of clusters of these gaussian mixtures. Under this model, we readily have

I 1 = |x -y| f X (x)f Y (y) dx dy = K i=1 K j=1 π i π j |x -y|ϕ(x; µ i ; σ i ) ϕ(y; µ j ; σ j ) dx dy = K i=1 K j=1 π i π j E |Z i -Z j |
where Z i and Z j denote independent gaussian variables with distributions N (µ i , σ i ) and N (µ j , σ j ). By classical properties of the gaussian family, Z i -Z j has distribution

N µ i -µ j , σ 2 i + (σ j ) 2
. Therefore, we use the following formula (proved in [START_REF] Worms | Parametric and Non-Parametric Statistics for Program Performance Analysis and Comparison[END_REF]): if Z has distribution N (µ, σ) then E [ |Z| ] = (2Φ(µ/σ)-1)µ+2σϕ(µ/σ). This entails the following formula for our theoretical performance metric:

I 1 = K i=1 K j=1 π i π j (µ i -µ j )   2Φ   µ i -µ j σ 2 i + (σ j ) 2   -1   + 2(σ 2 i + (σ j ) 2 ) π e -(µ i -µ j ) 2 /(2(σ 2 i + (σ j ) 2 ))
Consequently, using the estimations θ, K, θ , K of the parameters θ, K, θ , K (i.e. the parameters of the estimated gaussian mixtures distributions F GM X and F GM Y ), the parametric estimation I 1 of our first performance metric I 1 comes:

I 1 = K i=1 K j=1 πi π j (μ i -μ j )   2Φ   μi -μ j σ2 i + (σ j ) 2 )   -1   + 2(σ 2 i + (σ j ) 2 ) π e -(μ i -μ j ) 2 /(2(σ 2 i + (σ j ) 2 ))
B. The metric I 2 : the probability that a single program run is better than another

When the user needs to select which code version to execute, he may base his selection criteria on the expected average speedup for instance. But if his application is rarely executed, the average performance gain may not be interesting for him. He may be interested in executing a single time his application, and he wishes that his single run has the best chances of being the fastest between the two code versions. Formally, to help him to decide, we can compute P [X < Y ], the probability that a single run of X would be better than a single run of Y . This defines our second metric of program performances

I 2 = P [X < Y ].
Our parametric estimation assumes that both X and Y are modeled with gaussian mixture distributions.

I 2 = 1 x<y f X (x)f Y (y) dx dy = K i=1 K j=1 π i π j 1 x<y ϕ(x; µ i ; σ i ) ϕ(y; µ j ; σ j ) dx dy = K i=1 K j=1 π i π j P Z i -Z j < 0
where Z i and Z j denote independent gaussian variables such that Z i -Z j is gaussian distributed with expectation µ iµ j and variance σ 2 i + (σ j ) 2 . Since P [Z < 0] = Φ(-µ/σ) whenever Z has distribution N (µ, σ), we thus have

I 2 = K i=1 K j=1 π i π j Φ µ j -µi √ σ 2 i +(σ j ) 2 .
Consequently, plugging in the estimators of the parameters of the gaussian mixture distributions leads to the following parametric estimator of I 2 :

I 2 = K i=1 K j=1 πi π j Φ   μ j -μi σ2 i + (σ j ) 2  
Alternatively, we can also generalise this metric to consider a constant real shift ∆ ∈ R to check between X and Y , and thus consider I 2 = P [X < Y + ∆] with its parametric estimation:

I 2 = K i=1 K j=1 πi π j Φ   ∆ + μ j -μi ) σ2 i + (σ j ) 2  
C. The metric I 3 : the probability that a single run is better than all the others

In practice, a user may have more than only two code versions. How can he decide about the best code version among many others, for a single run only ? Comparing code versions two by two is misleading. All code versions must be compared together. Let X 1 , X 2 , • • • , X r denote r random variables corresponding to r distinct code versions. We propose to compute the probability that one code version, say the first one, executes faster than all the others for a single run only (not in average or in median), i.e. that X 1 < min(X 2 , . . . , X r ). This defines the following program performance metric:

I 3 = P [X 1 < min(X 2 , • • • , X r )] = E 1 X1<min(X2,••• ,Xr)
The parametric estimation of I 3 is noted I 3 . Here we assume that, for every given j ∈ {1, . . . , r}, the random variable X j is distributed as a gaussian mixture with parameters K = K j and θ = θ j = (π 1,j , . . . , π Kj ,j ; µ 1,j , . . . , µ Kj ,j ; σ 1,j , . . . , σ Kj ,j ).

As we did for I 1 and I 2 , we need to obtain a formula for the metric I 3 in terms of the parameters. Let us note Y = min(X 2 , • • • , X r ), and let G be the c.d.f. of Y . By the mutual independence of X 1 , X 2 , . . . , X r , the variables X 1 and Y are independent and therefore a classical probability property yields, since

P [x < Y ] = (1 -G)(x), I 3 = P [ X 1 < Y ] = E [1 X1<Y ] = E [(1 -G)(X 1 )] = ∞ -∞ (1 -G(x))f 1 (x)dx.
By independence of the variables X 2 , . . . , X r , and Equ. (2), we have:

(1 -G)(x) = P [ X 2 > x, . . . , X r > x ] = r j=2 P [X j > x] = r j=2 Kj i=1 π j (1 -Φ(x; µ i,j ; σ i,j )).
Our parametric estimator I 3 of the metric I 3 is then equal to the following integral (which we compute numerically, by using the R software for instance)

I 3 = ∞ -∞ (1 -Ĝ(x)) f1 (x)dx
where Ĝ(x) = r j=2 Kj i=1 πj (1 -Φ(x; μi,j ; σi,j )) and f1 (x) = K1 i=1 πi,1 ϕ(x; μi,1 ; σi,1 ) D. The metric I 4 : the variability level People do not always know how to quantify the variability of programs performances. By default, they use the variance, but they may not know how to interpret it. The variance measures how the data spread out around the average: but if the data are multi-modal or present clusters, then the average is not necessarily a good measure of the variability, especially when the modes or clusters are particularly distant from each other, and therefore the variance loses its attractiveness.

We propose to consider an alternative or complementary measure of the variability of programs performances: the number of modes of the underlying p.d.f. of the data, which we will note I 4 . It is simply equal to the number of local maxima of the p.d.f., which is supposed to represent the different values around which the data are spreading or clustering. A local maxima is called a mode in statistics.

Thanks to the gaussian mixture modelling, which yields an explicit formula for the estimated p.d.f., we can compute a parametric estimation I 4 of I 4 , which is simply equal to the number of local maxima of the gaussian mixture p.d.f. f θ, K estimated from the data. Often, this estimation I 4 turns out to be equal to K, the estimated number of clusters of the fitted gaussian mixture distribution. But it is not necessarily always the case, since sometimes the gaussian mixture fitting algorithm proposes a higher value of K than the actual number of groups in the data, in order to flexibly account for assymetry. Formally, the variability level is computed as: I 4 =number of local maxima of the estimated gaussian mixture p.d.f., which is often equal to K, the number of components in this model.

IV. IMPLEMENTATION AND EXPERIMENTS

We implemented all our statistical methods using the R software, and we delivered our software with a free academic licence. We did a huge and extensive set of experiments and simulations, all are detailed in [START_REF] Worms | Parametric and Non-Parametric Statistics for Program Performance Analysis and Comparison[END_REF]. We collected during more than 10 years a great amount of performance data, resulted from many empirical studies: SPEC CPU applications (2001,2006), all SPEC OMP applications, NAS Parallel Benchmark, own micro-benchmarks, other parallel applications, various compilers versions and options, Linux versions, different HPC machines architectures and generations, etc. The number of samples is 2438, each one contains between 30 and 1000 execution times. In this section we present the list of the most important conclusions.

The gaussian mixture model fits well the data:: The first empirical study was devoted to validating the robustness of our goodness-of-it test, based on Kolmogorov-Smirnov distances. We used extensive statistical simulations to check the robustness and the error ratio of the test and bootstrap calibration, the results are very satisfying compared to what is used in practice for other tests [START_REF] Worms | Parametric and Non-Parametric Statistics for Program Performance Analysis and Comparison[END_REF]. Based on this test, we are able to ensure that 83% of the samples are very well fitted with gaussian mixtures. The remaining 17% samples were rejected from gaussian mixture modelling. After investigation, we found that the reason is mostly because they contain ties (identical values), that result from rounding errors when collecting performance data. It remains however a marginal set of samples that cannot be modelled by gaussian mixtures, other distributions should be considered.

The performance metrics are accurate enough:: we also did an extensive set of experiments to check the accuracy of our new performance metrics I 1 , I 2 , I 3 and I 4 . We computed by simulations their mean square errors and mean absolute percentage errors, compared to the theoretical exact values I 1 , I 2 , I 3 and I 4 . The accuracy of our metrics are empirically satisfactory [START_REF] Worms | Parametric and Non-Parametric Statistics for Program Performance Analysis and Comparison[END_REF].

Empirical study of variability levels of programs execution times:: Based on our metric for the estimation of variability level ( I 4 ), we computed the value of this metric for all the samples. We found:

• ≈ 37% of the samples have a variability level equal to one, which means that the execution times are spread around a single value.

V. RELATED WORK AND DISCUSSION

Program performance evaluation:: In the field of code optimisation and high performance computing, most of the published articles declare observed speedups or other performance metrics. Unfortunately, few studies based on rigorous statistics are really conducted to check whether the observations of the code performance improvements are statistically significant or not.

Program performance analysis and optimisation may rely on two well known books that explain digest statistics to our community [START_REF] Jain | The Art of Computer Systems Performance Analysis : Techniques for Experimental Design, Measurement, Simulation, and Modelling[END_REF], [START_REF] Lilja | Measuring Computer Performance: A Practitioner's Guide[END_REF] in an accessible way. These two books are good introductions for doing fair statistics for analysing performance data. Based on these two books, previous work on statistical program performance evaluation have been published [START_REF] Georges | Statistically rigorous Java performance evaluation[END_REF]. In the later article, the authors rely on the Student's t-test to compare between two average execution times (the two sided version of the student t-test) in order to test if two theoretical means are unequal. We improved the previous work in [START_REF] Mazouz | Study of Variations of Native Program Execution Times on Multi-Core Architectures[END_REF]: first, we showed how to conduct a one-sided Student's t-test to validate that µ X > µ Y . Second, we show how to check the normality if small samples and the equivalence of their variances (using the Fisher's F-test) in order to use the classical Student's t-test instead of the Welch's variant.

In [START_REF] Touati | The Speedup-Test: A Statistical Methodology for Program Speedup Analysis and Computation[END_REF], we presented a rigorous statistical methodology regarding program performance analysis. We rely on well known statistical tests (Shapiro-wilk's test, Fisher's Ftest, Student's t-test, Kolmogorov-Smirnov's test, Wilcoxon-Mann-Whitney's test) to study if the observed speedups are statistically significant or not. By fixing 0 < α < 1 a desired risk level, we are able to analyse the statistical significance of the average execution time as well as the median.

In the current article, we go beyond the classical mean and median performance factors. We provide additional interesting performance metrics, based on parametric statistics (gaussian mixture modelling).

References on gaussian mixtures:: References about mixture models and particularly gaussian mixture models are numerous, in the statistics literature, machine-learning literature, as well as in many statistics-using fields (in particular image analysis, bioinformatics, biology, medicine, etc); we therefore only cite the reference book on mixture models [START_REF] Maclachlan | Finite Mixture Models[END_REF], [START_REF] Mengersen | Mixture Estimation and Applications[END_REF], [START_REF] Titterington | Statistical Analysis of Finite Mixture Distributions[END_REF].

The use of the EM algorithm as a solution to finite mixtures fitting is a classical subject in the statistics and pattern recognition literature, since the release of the breakthrough paper [START_REF] Dempster | Maximum likelihood from incomplete data via the em algorithm[END_REF] (note that the strength of the EM algorithm is that it is not restricted to the estimation of gaussian mixtures, but extend to other mixtures, often with a computational cost though). In the present work, we decided to address the important issue of choosing the appropriate number K of components by relying on the BIC criterion: there are however several popular alternatives, for which some references are [START_REF] Fujita | A nonparametric method to estimate the number of clusters[END_REF], [START_REF] Tibshirani | Estimating the number of clusters in a data set via the gap statistic[END_REF].

Discussion:: Concerning the performance of the goodness-of-fit test we introduced in [1], we found out that it was very satisfying. We nonetheless observed that, in the presence of too much a proportion of equal data in the dataset, our fitting test tends to artificially reject the gaussian mixture model too often: a simple solution to this problem would be to increase the precision of the measurement method so that the risk of observing exaequo is reduced.

Gaussian mixtures seem to be good model for most of the performances data that we observed in practice. Also, such data distributions provide interesting perspective regarding multi-dimensional performance data. Indeed, an interesting performance analysis must consider the performance as multi-dimensional data, each dimension corresponds to a specific performance nature. For instance we can consider the triplet {Execution time, energy consumption, network traffic}. Fortunately gaussian mixtures are a good solution for modelling multi-dimensional data in order to analyse the relationship between multiple dimensions. Some people are interested in extreme values statistics (Worst Case Execution Times, Best Case Execution Times). In that case, gaussian mixture modelling is not necessarily the adequate way of addressing this issue: either mixtures of other data distributions must be considered, or alternative methods should be considered (extreme value analysis techniques in particular, although they require quite a large number of data values to be truly reliable).

VI. CONCLUSION

In the presence of program performance variability, we must rely on formal statistics to decide about the best code version. Currently, people analyse mainly mean and median performances only. In the current research work, we extend these metrics as follows:

1) We build a statistical modelling based on gaussian mixtures to fit multi-modal data distributions; 2) We build a statistical test to quantify the quality of data-model fitting, based on Kolmogorov-Smirnov distances; 3) We define new code performance metrics that go beyond mean and median performances.; 4) We implemented all our statistical methods using R software, and we demonstrate its practical efficiency. 5) Our software, called VARCORE, is publicly distributed as a free open source code. Our gaussian mixture modelling provides a new and interesting metric to evaluate the variability of performances. Indeed, instead of only considering means and variances (the variances being, by the way, difficult to interpret in practice), we propose to consider as well the modes of the data distributions. Thus, the variability level of performances can be measured by the number of these modes, which we can compute with a parametric method based on gaussian mixtures. The number of modes is a natural way of giving an idea of the performance variability: a data distribution with a single mode means that the performances are quite stable around a single value; with two modes, it means that the performances are varying around two values, etc. In addition, if the number of modes is greater than one, this can be a good indication of being careful with the interpretation (and comparison) of the variance (since the usual interpretation generally assumes that the data are issued from a unimodal distribution). Moreover, the existence of these modes can be further investigated by trying to explain them with auxiliary measurements made during the execution of the program, which is certainly the most fruitful perspective of this research work.

  Figure 1. Gaussians mixture modelling

Figure 2 .

 2 Figure 2. Quality of data fitting (empirical vs. theoretical CDF)


Jean-Philippe Tavella

Mathieu Caujolle

Stephane Vialle

Cherifa Dad

Charles Tan

Gilles Plessis

Mathieu Schumann

Arnaud Cuccuru

Sebastien Revol

Toward an Accurate and Fast Hybrid Multi-Simulation with the FMI-CS Standard

Abstract-Multi-simulation in the context of future smart electrical grids consists in associating components modeling different physical domains, but also their local or global control. Our DACCOSIM multi-simulation environment is based on the version 2.0 of the FMI-CS (Functional Mock-up Interface for Co-Simulation) standard maintained by the Modelica Association. It has been specifically designed to run large-scale and complex systems on a single PC or a cluster of multicore nodes. But it is quite challenging to accurately simulate FMUs-composed systems involving predictable and unpredictable events while preserving the system overall performance. This paper presents some additions to the FMI-CS standard aiming to improve the accuracy and the performance of distributed multi-simulations involving a mix of both time steps and various kinds of events. The proposed FMI-CS primitives are explained, as well as the Master Algorithm strategies to exploit them efficiently.

I. MOTIVATIONS AND RELATED WORKS

A smart power grid is an example of complex Cyber-Physical Systems. For EDF (the major French utility company), it will extensively rely on new control functions (i) to increase the grid efficiency, reliability, and safety, (ii) to enable better integration of new assets (e.g. distributed generation and alternative energy sources), (iii) to support market dynamics and manage new interactions between established and new energy players. We consider vital to develop agile modeling and simulation in order to design and validate novel distributed operating functions a long time before performing tests on experimental sites, and to reuse models and tools often limited to specific areas of expertise.

Heterogeneity of smart power grids led EDF to investigate coupling standards such as the Functional Mock-up Interface (FMI) initiated by Daimler AG within the ITEA2 MODELISAR project and now maintained by the Modelica Association 1 . More precisely, EDF chose the FMI-CS (FMI for Co-Simulation) standard because this operation mode allows to export models as active components called FMUs (Functional Mock-up Units), each FMU being a self-contained archive file including a model and a numerical solver with the additional benefit of protecting the model intellectual property (IP). Simulation of a smart power grid requires to associate eventbased and time stepped components to achieve hybrid multi-1 https://www.fmi-standard.org/downloads simulations, but managing events with the current FMI-CS standard is not straightforward.

EPOCHS [START_REF] Hopkinson | EPOCHS: a platform for agent-based electric power and communication simulation built from commercial off-the-shelf components[END_REF] was an example of distributed multi-simulation environment based on the well known HLA middleware standard2 . It mixed events and time steps, but delayed event processing at the end of each computation time step (in simulated time) and lacked accuracy. At the opposite, GECO [START_REF] Lin | GECO: Global Event-Driven Co-Simulation Framework for Interconnected Power System and Communication Network[END_REF] accurately processed all events in a global event queue without introducing delay in the simulated time. But it did not support distribution on several PCs, and so can not scale largely. In 2013, some researchers started to mix distributed HLA standard and FMI 1.x standards [START_REF] Awais | The high level architecture RTI as a master to the functional mock-up interface components[END_REF] but could not manage events before the end of the time steps. In 2014 the C2WT environment mixed again HLA and FMI standards [START_REF] Neema | Model-Based Integration Platform for FMI Co-Simulation and Heterogeneous Simulations of Cyber-Physical Systems[END_REF], and interconnected the telecom network simulator OMNET++ with FMUs modeling continuous time systems, but it used only constant time steps and its execution algorithm did not address time step interruption for unexpected event processing.

In 2016, some members of our consortium proposed a generic computation strategy and algorithm [START_REF] Camus | Hybrid Co-simulation of FMUs using DEV&DESS in MECSYCO[END_REF] to achieve accurate event location while respecting the FMI-CS 2.0 standard. But despite optimization, the current FMI-CS features can lead to numerous rollbacks and reruns of FMU computations to process unexpected events. So, such a solution could limit the performance and scalability achieved on multi-core PC clusters by our pure FMI-based distributed multi-simulations [START_REF] Dad | Scaling of Distributed Multi-Simulations on Multi-Core Clusters[END_REF]. Moreover, discrete signals often used by cyber-physical components in electrical grids cannot be well addressed by this strategy.

This paper proposes new FMI-CS functionalities, in order to improve both the accuracy of cyber-physical co-simulations involving a mix of time steps and various kinds of event signals, and their performance and scalability on distributed architectures.

II. FMI-CS STANDARD AND DACCOSIM ENVIRONMENT

A. Functional Mockup Interface features and limits

In FMI-CS a component is a self-contained object that includes the model description and a numerical solver provided by the design environment where it was conceived. Cosimulation can encompass a large number of active components (FMUs) potentially developed from different tools. The data exchange between these components is controlled by a program called Master Algorithm (MA) whose behavior is not specified by the FMI standard. The MA controls the time progress of a system co-simulation with a sequence of step by step calculations. During a computation step, each FMU independently simulates a part of the system and, at the end of the step, the outputs from some FMUs provide new input values to some other FMUs.

However, our practice shows that all kind of signals are not correctly handled by the current version 2.0 of the FMI-CS standard, especially in an hybrid co-simulation context mixing discrete and continuous components (see section III-B).

B. DACCOSIM, an FMI-based co-simulation software

DACCOSIM has been designed to achieve multisimulations of continuous time systems discretized with time steps, and run solvers using constant or variable time steps, the latter maintaining the requested accuracy with the minimal amount of computations whatever the system dynamics. DACCOSIM parallel and distributed runtime architecture is performance-and scalability-oriented, and takes maximal advantage of any cluster of multi-core nodes. Each simulation executes a series of time steps composed of three stages: the time step computations of the FMUs (independent computations), the communication of the computed FMUs outputs values to the connected inputs (many small communications), and the simulation control by the hierarchical control masters (information gathering, next operations decision and order broadcasting). DACCOSIM architecture distributes FMUs on cluster nodes, encapsulates each FMU in a multithreaded wrapper, and implements a hierarchical (and distributed) control master to manage all threads and FMU operations (Fig. 1). A unique global master controls local ones, and each local master controls a set of FMUs on its node. All master tasks run concurrently. Each multithreaded wrapper runs the computation function of its FMU to achieve time step progress, manages inter-FMUs and FMU-Master communications, and stores simulation results on disk. It parallelizes and optimizes the orchestration of these tasks to minimize the global execution time. The Java version of DACCOSIM relies on JavaFMI3 developed by SIANI institute (Las Palmas University), and is available for both Windows and Linux operating systems, either 32-bit or 64-bit.

III. FMI EVOLUTION TO CONSIDER VARIOUS SIGNALS

A. Almost everywhere differentiable signals Fig. 2 illustrates the different kinds of almost everywhere differentiable signals we need to consider in our cyber-physical systems. Signal s1 is a continuous & piecewise differentiable signal. It is present at each time t i ∈ R + , continuous on R + , but not differentiable at all the points t j no having a null ordinate. Signal s2 is a piecewise continuous & differentiable signal. It is present at each time t i ∈ R + , but is neither continuous nor differentiable at time t1 and t2. Signal s3 is a piecewise constant signal. It is present at each time t i ∈ R + , and is constant on disjoint and continuous time slots I j so that j I j = R + , with a discontinuity appearing at each time slot switch t j . Signal s4 is a discrete event signal. It is present for a set of definition D, being a discrete time set t i ∈ D, with D ⊂ R + . This signal can be confused with the events it generates.

Later in this article, we only focus on signals that are always present all the time of the simulation (signal s1, s2 and s3 on figure 2) as these cases correspond to the ones encountered in cyber-physical systems like smart grids, smart buildings and smart factories.

B. Proposal of FMI-CS standard evolution

The FMI-CS standard defines an application programming interface (API) that components conforming to the standard must implement. The environment interacts with the FMU only by means of its API, which consists primarily of the fmi2DoStep() primitive which attempts to perform a simulation step on the FMU, given a specific time step. The problem with the current version of the FMI-CS standard is that event handling is not very accurate: when an event occurs during a step calculation, the event is only reported at the end of the current step. If the simulation step is too large, events could even be missed. Different strategies can be considered to approximately locate an event (progress with minimum step size [START_REF] Galtier | FMI-Based Distributed Multi-Simulation with DAC-COSIM[END_REF], bisectional method [START_REF] Camus | Hybrid Co-simulation of FMUs using DEV&DESS in MECSYCO[END_REF]). All of them require rollbacks and decrease the overall computational performance (Fig. 3).

Specifically, events can be either state events or time events. It is assumed here that state events are unpredictable breakpoints (they cannot be predicted without exploring the future) while time events are predictable break-points.

For unpredictable break-points, an evolution of the fmi2DoStep(...) is proposed in the form of a new primitive fmi21DoStep(stepSize, nextEventTime). It allows the FMU solver to stop at the first unpredictable event with a new return code fmi21Event and the time event is given by the returned value nextEventTime (Fig. 4). If no event occurs, nextEventTime is the computation step end time. So, the MA can go on computing state variables without performing any rollback on this component.

But when one FMU prematurely stops, the other FMUs involved in the co-simulation will proceed with their calculation causing a computation time waste. The fmi21BreakPendingStep(earlierTime, stopTime) primitive we propose is meant to reduce this computational burden by interrupting an FMU pending calculation before its normal end. A pending step execution will be stopped by an asynchronous call of the primitive at earlierTime if it is not yet reached, or upon receipt in the opposite case returning the actual stop time value stopTime. It will also ignore any earlierTime earlier than its current start time.

For predictable break-points, a primitive named fmi21GetNextEventTime(currentTime, stopTime, eventTime)

is proposed to report in advance the time event precise date (Fig. 5). This time information could be exploited to optimize the value of the next step size in view of the exact date of the next predictable event.

To be complete, we propose four additional primitives fmi21GetXXX() to get the value of variables according to their type XXX (Boolean, Integer, Real or String) at a discontinuity point. Prototyping of these primitives makes sense in particular in the context of super-dense time system [START_REF] Liu | Modeling timed concurrent systems[END_REF], but not only. After a fmi21DoStep(), the FMU stops at t + i , allowing the new fmi21GetXXX() to get the values of variables (at t + i ) after all the events occurring at the same Newtonian time t, exactly as the current fmi2GetXXX() do. But in addition and without changing the reached Newtonian time all previous signal values are also returned by the new primitive whatever the events that happened between t - i and t + i . This mechanism respects the constraint of a positive non zero step size for step computation as required by the current FMI standard but contradicts a recent proposal from UC Berkeley [START_REF] Cremona | FIDE -An FMI Integrated Development Environment[END_REF].

All these new primitives have been introduced by EDF in June 2016 in the Working Group Clock and Hybrid Cosimulation of the FMI project, now hosted by the Modelica Association4 , and discussions are ongoing.

IV. EXPECTED IMPACT ON MULTI-SIMULATION

A. Illustrative use-case presentation

The academic use case we intend to use in order to implement and test the proposed FMI primitives is a barrelfiller factory (Fig. 6). Involved models, detailed in [START_REF] Camus | Hybrid Co-simulation of FMUs using DEV&DESS in MECSYCO[END_REF], consist of:

• A queue of barrels waiting on a conveyor to be filled.

• A tank storing the water to fill the barrels.

• A controller c1 managing the opening of the valve between the tank and the barrel. • A controller c2 regulating the whole filling process. It mixes continuous dynamics (water flow from the tank) and several types of non-continuous signals discussed in III-A, such as piecewise continuous and differentiable signals (water level in the tanks), piecewise constant and potentially discrete events (command signals), see Fig. 7.

B. Expected impact on result accuracy

Results showed in Fig. 8 come from a DACCOSIM multisimulation. They highlight the current accuracy insufficiencies of FMI-CS 2.0. Even when applying adaptive time step strategies relying on rollbacks to ensure accurate simulation and correct detection of events, the location of events can only be approximated. Accuracy is limited by the minimum time step size authorized by the masters. But, when using the proposed fmi21DoStep(), the accuracy of the event time will no longer be limited by the master algorithm control strategies, it will only depend on the ability of the FMU internal solver to manage these events.

The proposed fmi21GetXXXEvent() primitives complement the fmi21DoStep() function by allowing to fetch the FMU output values computed by the internal solver just before and after a time event. Exact characterization of the FMU dynamics before and after an event can thus be achieved without having to approximate them by performing small step computations. Signals such as the barrel water volume or the abort signal (Fig. 7-a and-d) will thus be accurately characterized.

Thanks to the proposed add-ons to the FMI-CS standard, similar accuracy performance should be achieved within a domain-specific simulation tool and within a FMI 2.0 based generic multi-simulation tool such as DACCOSIM.

C. Expected impact on performances and parallelization

One of the major challenges when managing events with the FMI standard in hybrid multi-simulations is avoiding CPU time waste. The new FMI functions we propose (see section III-B) will help to avoid some rollbacks and re-computations:

• With function fmi21GetNextEventTime() a control FMU will be able to communicate the date of its next known time event. Master algorithms, such as the one implemented in DACCOSIM, should be able to exploit this information when defining the size of the new computation step: its end will be aligned with the exact event time if initially included in the considered time step. So, at the end of the time step, the event will be routed to other FMUs without any delay nor rollback. • With function fmi21DoStep(), an FMU will be able to stop its computation just after it has generated a state event, and with function fmi21BreakPendingStep() it will be able to request all other FMUs they stop their computations as soon as they reach this event time. Thus the new DACCOSIM master algorithms will take advantage of this behavior and ask only the fastest FMUs, i.e. those whose computation process had already exceeded the earliest event time before getting the corresponding break signal, to roll back and rerun a computation up to the event time. Fig. 9 illustrates this last mechanism, considering 3 FMUs running their doStep() functions concurrently on different CPU cores. FMU1 generates a state event at the simulated time t i+1 and signals to all others FMU to stop their computation at t i+1 (calling fmi21BreakPendingStep()). But FMU2 achieves very light computations and has already reached the end of the simulated time step (t i+1 ). Then, FMU2 will roll back to t i and re-run its computation up to t i+1 . At the opposite, FMU3 is more computation intensive and benefits from the FMU2 break signal to stop its computation when reaching t i+1 , avoiding unnecessary computations (up to t i+1 ), rollback and re-computation (up to t i+1 ). Finally, complete execution time on a multi-core processor will be shorter when resizing the ongoing FMU computations. The fmi21BreakPendingStep() function is also under discussion in the Working Group Clock and Hybrid Co-simulation of the FMI project.

Better support of events in FMI-CS will allow to introduce more FMU components devoted to control. These components consume very few CPU time (fast control algorithms) but increase communications (reading many sensor outputs and sending command signals). So, this evolution should impact FMU distribution algorithms on PC clusters, which optimize both load balancing and inter-node communications [START_REF] Dad | Scaling of Distributed Multi-Simulations on Multi-Core Clusters[END_REF].

V. FUTURE WORKS ON SMART GRID SIMULATIONS

Until now the FMI-CS standard allowed us to compute different kinds of Smart Grid simulation from building thermics [START_REF] Dad | Scaling of Distributed Multi-Simulations on Multi-Core Clusters[END_REF] to MV network voltage management based on stateestimation. These simple use cases permitted us to identify the limits of the current standard and design new ways to overcome them. In fact, integrating the proposed primitives into the standard will be required in order to enable efficient and accurate large, multi-physics and control-centric simulations.

In a near future, complex Smart Grid use cases involving (i) distribution networks to control Distribution System Operators (DSO) assets with new smart functions to ensure the absence of constraints on the power grid, (ii) buildings and their energy systems with control over elements such as heating, ventilation and air conditioning, and (iii) energy actors to control the flexibility of buildings and electrical vehicles, will be performed (Fig. 10). So, important work is going to be A preliminary version of this work has been presented in May 2016 in the Annex 60 project, an international project conducted under the umbrella of the International Energy Agency (IEA) within the Energy in Buildings and Communities (EBC) Programme. Annex 60 will develop and demonstrate new generation computational tools for building and community energy systems based on Modelica, FMI and BIM standards, contributing to smart power grid design. Then, a more complete version has been introduced in June 2016 in the Working Group Clock and Hybrid Co-simulation of the FMI project, and discussions are ongoing to attempt to improve the FMI-CS standard.

Fig. 1 .

 1 Fig. 1. DACCOSIM distributed architecture

Fig. 2 .

 2 Fig. 2. Different kinds of almost everywhere differentiable signals

Fig. 4 .Fig. 5 .

 45 Fig. 4. Behavior of the new primitive fmi21DoStep(...)

Fig. 6 .Fig. 7 .

 67 Fig. 6. Use case mixing discrete and continuous components

Fig. 9 .

 9 Fig. 9. Parallel execution of FMUs with different CPU loads, on a multi-core processor

http://standards.ieee.org/findstds/standard/1516-2000.html

https://bitbucket.org/siani/javafmi/wiki/Home

https://www.modelica.org/

ACKNOWLEDGMENT

Authors thank Region Grand Est and RISEGrid institute for their support to this research.