Abir Ben

Khaled-El Feki

Laurent Duval
email: laurent.duval@ifpen.fr

Cyril Faure
email: cyril.faure@cea.fr

Daniel Simon
email: daniel.simon@inria.fr

Mongi Ben
email: mongi.ben-gaid@ifpen.fr

Ben Mongi

Gaid

CHOPtrey: contextual online polynomial extrapolation for enhanced multi-core co-simulation of complex systems

Keywords: parallel simulation, Functional Mockup Interface, smoothed online prediction, causal polynomial extrapolation, context-based decision, internal combustion engine

  L'archive ouverte pluridisciplinaire

Introduction

Intricacy in engineered systems increases simulator complexity. However, most existing simulation software products are currently unable to exploit multi-core platforms, as they rely on sequential Ordinary Differential Equations (ODE) and Differential Algebraic Equations (DAE) solvers. Co-simulation approaches can provide significant improvements by allowing to jointly simulate models coming from different areas, and to validate both individual behaviors and their interactions [START_REF] Valášek | Modeling, simulation and control of mechatronical systems[END_REF]. Different simulators may be exported from original authoring tools, for instance as Functional Mock-up Units (FMUs), and then imported in a co-simulation environment. Hence, they cooperate at run-time, thanks to Functional Mockup Interface (FMI) definitions [2] for their interfaces, and to the master algorithms of these environments.

Co-simulation has shown an important potential for parallelization (see [START_REF] Ben | Distributed real-time simulation of numerical models: application to power-train[END_REF] for a review). Meanwhile, synchronization between the different sub-models is required due to their mutual dependencies. It avoids the propagation of numerical errors in simulation results and guarantees their correctness. Unfortunately, synchronization constraints also lead some processors into waiting periods and idle time. This decreases the potential efficiency of the threaded parallelism existing in multi-core platforms.

To overcome this limitation and exploit the available parallelism more efficiently, the dependency constraints between parallel sub-models should be relaxed as far as possible while preserving an acceptable accuracy of the simulation results. This can be performed by a wellgrounded system decomposition, tailored to data dependency reduction between the sub-models. For instance, the method proposed in [START_REF] Sjölund | Towards efficient distributed simulation in Modelica using Transmission Line Modeling[END_REF] for distributed simulation uses transmission line modeling (TLM) based on bilateral delay lines: decoupling points are chosen when variables change slowly and the time-step of the solver is relatively small.

Unfortunately, often perfect decoupling cannot be reached and data dependencies still exist between the different blocks. Thus, tight synchronization is required between blocks using small communication steps. This greatly limits the possibilities to accelerate the simulation and eventually to reach the real-time.

We propose in this work a Computationally Hasty Online Prediction framework (CHOPred) to stretch out synchronization steps with negligible precision changes in the simulation, at lowcomplexity. It is based on a Contextual & Hierarchical Ontology of Patterns (CHOPatt) that selects appropriate Causal Hopping Oblivious Polynomials (CHOPoly) for signal forecasting, allowing data exchange between sub-models beyond synchronization points.

This paper is organized as follows. We first review challenges as well as related work and summarize contributions in Section 2. Section 3 presents a formal model of a hybrid dynamical system and the motivations for twined parallel simulation and extrapolation. The background on prediction and the proposed Causal Hopping Oblivious Polynomials (CHOPoly) are developed in Section 4, with details in Appendix A. Then, the principles of the Contextual & Hierarchical Ontology of Patterns (CHOPatt) for the management of hybrid models are exposed in Section 5. Finally, the methodology's performance is assessed in Section 7 using an internal combustion engine model described in Section 6.

Related work and contributions

Continuous systems usually attain high integration speeds with variable-step solvers. The major challenge for hybrid systems resides in their numerous discontinuities that prevent similar performance. Since we are especially interested in modular co-simulation [START_REF] Ben | Distributed real-time simulation of numerical models: application to power-train[END_REF], it is shown in [5] that integrating each sub-system with its own solver allows avoiding interrupts coming from unrelated events. Moreover, event detection and location inside a sub-system can be processed faster because they involve a smaller set of variables. However, partitioning the original complex model into several lesser complex models may add virtual algebraic loops, therefore involving delayed outputs, even with an efficient execution order.

To take advantage of model splitting without adding useless delays, we propose in [6] a new co-simulation method based on a refined scheduling approach. This technique, denoted "RCosim", retains the speed-up advantage of modular co-simulation thanks to the parallel execution of the system's components. Besides, it improves the accuracy of the simulation results through an offline scheduling of operations that takes care of model input/output dynamics.

However, in practical applications, current co-simulation set-ups use a constant communication grid shared by all the models. In fact, the size of communication steps has a direct impact on simulation errors, and effective communication step control should rely on online error estimation induced by slackened exchange rates. Schierz et al. [START_REF] Schierz | Co-simulation with communication step size control in an FMI compatible master algorithm[END_REF] propose the use of adaptive communication step sizes to better handle the various changes in model dynamics. Meanwhile, the stability of multi-rate simulators with adaptive steps needs to be carefully assessed, for example, based on errors propagation inside modular co-simulations [8].

Data extrapolation over steps is also expected to enhance the simulation precision over large communication steps. Nevertheless, actual complex systems (CPS) usually present non-linearities and discontinuities. Therefore, their future behavior is hard to predict from past observations only. Moreover, the considered models are generated using the Simulink Coder target or the FMI for Model Exchange framework, which does not provide input derivatives (in contrast with the FMI for Co-Simulation architecture). Hence, polynomial prediction cannot always correctly extrapolate. For example, [9] uses a constant, linear or quadratic extrapolation and a linear interpolation to improve the accuracy of the modular time integration. This method is successful for non-stiff systems but fails in the stiff case. Data extrapolation over steps is also expected to enhance the simulation precision over large communication steps. Nevertheless, actual complex systems (CPS) usually present non-linearities and discontinuities. Therefore, their future behavior is hard to predict from past observations only. Moreover, the considered models are generated using the Simulink Coder target or the FMI for Model Exchange framework, which does not provide input derivatives (in contrast with the FMI for Co-Simulation architecture). Hence, polynomial prediction cannot always correctly extrapolate. For example, [9] uses a constant, linear or quadratic extrapolation and a linear interpolation to improve the accuracy of the modular time integration. This method is successful for non-stiff systems but fails in the stiff case.

In the purpose of defining a method for the parallel simulation of hybrid dynamical systems, our previous work on a single-level context-based extrapolation [START_REF] Khaled | Context-based polynomial extrapolation and slackened synchronization for fast multi-core simulation using FMI[END_REF] accounts for steps, stiffness, discontinuities or weird behaviors. It uses adapted extrapolation to limit excessively wrong prediction. It shows that properly-chosen context-based extrapolation, combined with model splitting and parallel integration, can potentially improve the speed versus precision trade-off needed to reach real-time simulation.

In [START_REF] Khaled | Context-based polynomial extrapolation and slackened synchronization for fast multi-core simulation using FMI[END_REF], context-based extrapolation is exclusively intended for FMU models and extrapolation is performed on integration steps only. This paper improves upon the preceding results with a better-grounded methodology, by:

• adding oblivion to past samples through a weighting factor in polynomial prediction;

• increasing the computational efficiency via a novel matrix formulation without inversion;

• adding an online error evaluation at each communication step to select the best context and weighting factor at the forthcoming step, further minimizing extrapolation errors;

• adding a hierarchy of decisional and functional contexts to improve prediction strategy.

Problem formalization and motivation

Model definition

Complex physical systems are generally modeled by hybrid non-linear ODEs or DAEs. The hybrid behavior is due to the discontinuities, raised by events triggered off by the crossing of a given threshold (zero-crossing). It plays a key role in the simulation complexity and speed. Indeed, more events slow down numerical integration.

Let us provide a formal model, considering a hybrid dynamic system Σ whose continuous state evolution is governed by a set of non-linear differential equations:

Ẋ = f (t, X, D, U ext ) for t n ≤ t < t n+1 , Y ext = g(t, X, D, U ext ) , where X ∈ R n X is the continuous state vector, D ∈ R n D is the discrete state vector, U ext ∈ R n U ext is the external input vector, Y ext ∈ R n Y
is the external output vector and t ∈ R + denotes the time. The sequence (t n ) n≥0 of strictly increasing time instants represents discontinuity points called state events, which are the roots of the equation:

h(t, X, D, U ext ) = 0 .
The function h is usually called zero-crossing function or event indicator, used for event detection and location [START_REF] Zhang | Zero-crossing location and detection algorithms for hybrid system simulation[END_REF]. At each time instant t n , a new continuous state vector can be computed as a result of the event handling:

X(t n ) = I(t n , X, D, U ext ) ,
and a new discrete state vector can be computed as a result of a discrete state update:

D(t n ) = J(t n-1 , X, D, U ext ) .
If no discontinuity affects a component of X(t n ), the right limit of this component will be equal to its value at t n . This hybrid system model is adopted by several modeling and simulation environments and is underlying the FMI specification [2]. We assume that Σ is well-posed, in the sense that a unique solution exists for each admissible initial conditions X(t 0 ) and D(t 0 ) and that consequently X, D, U ext , and Y ext are piece-wise continuous functions, i.e. continuous on each sub-interval ]t n , t n+1 [.

Model parallelization with modular co-simulation

To execute the system in parallel, the model must be split into several sub-models. For simplicity, assume that the system is decomposed into two separate blocks denoted Model 1 and Model 2, in Figure 1 with X = [X [1] X [2] ] T and D = [D [1] D [2] ] T , where T denotes the matrix transpose. Therefore, the sub-systems can be written as: [1] (t, X [1] , D [1] , U [1] , U ext ) , Y [1] = g [1] (t, X [1] , D [1] , U [1] , U ext ) ,

Ẋ[1] = f
Ẋ[2] = f [2] (t, X [2] , D [2] , U [2] , U ext ) , Y [2] = g [2] (t, X [2] , D [2] , U [2] , U ext ) . Y [1] U [2] Y [2]
U [1] Model 1

Model 2

Global model

Global model Σ

U ext Y ext U ext Y ext Model split Σ 1 Σ Σ 2
Figure 1: System splitting into sub-models for parallelization.

Here, U [1] are the inputs needed for Model 1 (Σ 1 ), directly provided by the outputs Y [2] produced by Model 2 (Σ 2 ). Similarly, U [2] are the inputs needed for Model 2 directly provided by the outputs Y [1] produced by Model 1. Our approach generalizes to any decomposition of the system Σ into B blocks, where each block is indexed by b ∈ {1, . . . , B}.

Model of computation

We define here a simulator as a model equipped with a solver. Different simulators run in parallel to perform the numerical integration of the whole multi-variable system. Figure 2 depicts an instance based on the two sub-models illustrated in Figure 1. Each of these simulators needs to exchange, at communication (or synchronization) points t s b (b = 1 or b = 2), the data needed by the other. To speed up integration, the parallel branches must be as independent as possible, so that they are synchronized at a rate H [b] = t s b +1t s b , by far slower than their internal integration step h [b] n b (H [b] h [b] n b ). Therefore, between communication points, each simulator integrates at its own rate (assuming a variable-step solver).

Σ 1 Σ 2

Initialization Exchange 1 Exchange 2

Integration step h [1] n 1

Communication step H [1] Communication step H [2] Integration step h Special case: H [1] =H [2] =H t s 1 =t s 2 =t s H = t s+1 -t s Two issues arise whether in the mono-rate or the multi-rate case. In the first case (Figure 2), each simulator considers that data incoming from other simulators is held constant within the common exchange rate. Thus, they all suffer from similar splitting-induced update errors.

In the multi-rate case illustrated in Figure 3, Model 1 needs to communicate with an external model two times faster than Model 2 (H

[b 1 ] = 1 2 H [b 2 ]
). Data incoming from Model 2 is thus held constant during 2H [1] , potentially causing aliasing. Faster simulators are additionally impacted by slower ones.

Σ 1 Σ 2

Initialization Exchange 1 Exchange 2

Integration step h [1] n 1

Communication step H [1] Communication step H [2] Integration step h H [2] =2.H [1] t s 2 +1 =t s 1 +2

t s 1 +1 t s 1 +3
, t s 2 +2 =t s 1 +4 , etc. While larger and multi-rate communication intervals are likely to speed up simulation, they may result in increased integration errors and poorer confidence in the final result.

External model

Motivation for extrapolation

Modeling the errors induced by slackened synchronization is a first direction in improving the trade-off between integration speed and accuracy.

Integration errors and parallelism

Error evaluation and convergence analysis were performed in [START_REF] Ben | Distributed real-time simulation of numerical models: application to power-train[END_REF]Chapter 9] for different schemes (sequential and parallel modular co-simulation, models with real and artificial loops, mono-rate and multi-rate co-simulation, etc.). Assuming that the b th sub-model is connected to B c -1 other sub-models, its global error on states (1) and outputs (2) are bounded in Landau-Bachmann O() notation. Errors are functions of the communication step (or synchronization interval) H = max H [b] , b ∈ {1, . . . , B c }, the integration time-step h and the order of accuracy p of the used numerical solver [3, Chapter 4]:

X [b] (t n+1 ) -X [b] n+1 ≤ O(h p ) + O(H) , (1)
Y [b] (t n ) -Y [b] n ≤ O(h p ) + O(H) . (2)
Both global errors are clearly bounded by two terms. The first term is related to the applied numerical solver, more specifically the time-step and the order. The worst case scenario on the bounding would be the maximum used integration step h and order p. The second term is related to the size of the communication step H. However, the weight of each term on the error is deeply related to the size of the communication step relatively to the integration step. Based on the same approach as in [9], it is clear that the communication step H dominates the error when H h. Therefore, considering a split model and a parallel execution, a trade-off must be found between acceptable simulation errors, thanks to tight enough synchronization, and simulation speed-up thanks to decoupling between sub-models.

Contribution of the extrapolation

To add a degree of freedom to this trade-off achievement, we propose to extrapolate model inputs to compensate the communication steps stretching between sub-models. In fact, it was proven previously that numerical solutions are first order accurate, O(H), when choosing larger communication steps H. Holding inputs constant between two synchronization intervals plays the role of a zeroth-order hold (ZOH, constant extrapolation). To generalize the error bound, the O(H) term can then be replaced with O H δ+1 , where δ is the extrapolation degree. Using for example linear (δ = 1) or quadratic (δ = 2, Appendix A.1) extrapolation instead of constant update (δ = 0) has the potential to reduce the bound of simulation errors, as proven in [9, p. 238 sq.]. However, the difficulty with extrapolation is that it could be sensitive for different reasons:

• there exists no universal prediction scheme, efficient for every signal;

• prediction should be efficient: causal, sufficiently fast and reliable;

• standard polynomial prediction may fail in stiff cases [12] (cf. Section 5 for details).

We choose to base our extrapolation on polynomial prediction, which allows fast and causal calculations. In this situation, the rationale is that the computing cost of a low-order polynomial predictor would be by far smaller than the extra model computations needed by shorter communication steps. Since such predictions would be accurate neither for any signal (for instance, blocky versus smooth signals) nor for any signal behavior (slow variations versus steep onsets), we borrow the context-based approach from lossless image encoders [START_REF] Weinberger | The LOCO-I lossless image compression algorithm: principles and standardization into JPEG-LS[END_REF], such as GIF (Graphics interchange format) or PNG (Portable Network Graphics) compression formats. The general aim of these image coders is to predict a pixel value based on a pattern of causal neighboring pixels. Compression is achieved when the prediction residues possess smaller intensity values, and more generally a sparser distribution (concentrated around close-to-zero values) than that of pixels in the original image. They may, therefore, be coded on smaller "bytes", using entropy coding techniques. In images, one distinguishes basic "objects" (smooth-intensity varying regions, edges with different orientations). Based on simple calculations over prediction patterns, different behaviors are inferred (e.g. flat, smooth, +45 o or -45 o edges, etc.). Look-up table predictors are then used, depending on the context.

In the proposed approach, we build a heuristic table of contexts (Section 5) based on a short frame of past samples, and affect pre-computed polynomial predictors to obtain contextdependent extrapolated values. We now review the principles of extrapolation.

CHOPoly: Causal Hopping Oblivious Polynomials for low-complexity extrapolation

Background on prediction, real-time constraints and extrapolation strategies

The CHOPtrey framework requires a dedicated instance of forecasting for discrete time series. The neighboring topics of prediction, interpolation or extrapolation represent a large body of knowledge in signal processing [START_REF] Wiener | Extrapolation, interpolation, and smoothing of stationary time series[END_REF]15], econometrics [START_REF] Brown | Smoothing, Forecasting and Prediction of Discrete Time Series[END_REF], control [START_REF] Box | Time Series Analysis: Forecasting and Control, Probability and Statistics[END_REF][START_REF] Väliviita | Polynomial predictive filtering in control instrumentation: a review[END_REF] and numerical analysis [START_REF] Richardson | The approximate arithmetical solution by finite differences of physical problems involving differential equations, with an application to the stresses in a masonry dam[END_REF]. They are paramount in complex systems that live, sample and communicate at different rates. Their use in simulation might be milder. There exists a natural and common belief that signal extrapolation may introduce additional error terms and cause numerical instability [START_REF] Arnold | Error analysis and error estimates for co-simulation in FMI for model exchange and co-simulation V2.0[END_REF][START_REF] Stettinger | A model-based approach for prediction-based interconnection of dynamic systems[END_REF]. Building upon our previous work [START_REF] Khaled | Context-based polynomial extrapolation and slackened synchronization for fast multi-core simulation using FMI[END_REF], two additions diminish the importance of extrapolation caveats: past sample weighting (oblivious polynomials) and a hierarchy of contexts defined by a pattern ontology.

As we operate under real-time conditions, implying strong causality, only a tiny fraction of time series extrapolation methods are practically applicable. For instance, the theory of multirate filter banks bears formal similarities with distributed simulation or co-simulation, and has long been deployed for interpolation, extrapolation, and reduction of computational needs [22]. Such systems allow optimized noise handling at variable sampling rates [START_REF] Gauthier | Optimization of synthesis oversampled complex filter banks[END_REF]. The overall propagation delay is harmless in traditional signal and image processing (for compression or signal restoration). However, it hinders their use in distributed simulation. Additionally, the decomposition of signals into frequency bands is not adapted to simulation data, that sometimes exhibit stiff behavior, preventing the use of band-limited extrapolation [START_REF] Liu | Finite length band-limited extrapolation of discrete signals[END_REF]. We alternatively propose to feed a bank of weighted smoothing polynomials, whose outputs are selected upon both simulation and real-time behaviors, organized in a Contextual & Hierarchical Ontology of Patterns (CHOPatt, Section 5).

ZOH (zeroth-order hold) or nearest-neighbor extrapolation is probably the most natural, the less hypothetical, and the less computationally expensive forecasting method. It consists in using the latest known sample as the predicted value. It possesses small (cumulative) errors when the time series is relatively flat or when its sampling rate is sufficiently high, with respect to the signal dynamics. In other words, it is efficient when the time series is sampled fast enough to ensure small variations between two consecutive sampling times. However, it indirectly leads to under-sampling or aliasing related disturbances [25]. They affect the signal information content, including its derivatives, and consequently its processing, for instance through differential equation systems. They appear as quantization-like noise [26, p. 609 sq.], delay induction, offset or bump flattening. As a remedy, higher order extrapolation schemes can be used. They consist in fitting a δ-degree polynomial through δ + 1 supporting points [27, p. 15 sq.]. The most common are the predictive first-order [START_REF] Arnold | Multi-rate time integration for large scale multibody system models[END_REF][START_REF] Hoher | A contribution to the real-time simulation of coupled finite element models of machine toolsa numerical comparison[END_REF] or the second-order hold (FOH or SOH) [25]. They result in linear or parabolic extrapolation. Polynomial parameters are obtained in a standard way with Waring-Lagrange interpolation. Hermite polynomials can be used when one is interested in extrapolating derivatives as well, and their relative stability has been studied [30, p. 61 sq.]. Both, although being exact at synchronization points, tend to oscillate amidst them. To avoid introducing discontinuities, extrapolation can be smoothed with splines [START_REF] Oh | A co-simulation framework for power system analysis[END_REF] or additional interpolation [START_REF] Busch | Continuous approximation techniques for co-simulation methods: Analysis of numerical stability and local error[END_REF].

In our co-simulation framework, communication intervals are not chosen arbitrarily small for computational efficiency. Thus, the slow variation of inputs and outputs cannot be ensured in practice. Provided a cost versus error trade-off is met, borrowing additional samples from past data and using higher-order extrapolation schemes could be beneficial. Different forecasting methods of various accuracy and complexity may be efficiently evaluated. We focus here on online extrapolation with causal polynomials, for simplicity and ease of implementation, following initial work in [START_REF] Faure | Real-time simulation of physical models toward hardware-in-the-loop validation[END_REF]Chap. 16]. Their main feature is that they are obtained in a least-squares fashion: we do not impose that they pass exactly through supporting points. However, they do so when the degree δ and the number of supporting points δ + 1 is set as above. We improve on [START_REF] Khaled | Context-based polynomial extrapolation and slackened synchronization for fast multi-core simulation using FMI[END_REF] by adding a time-dependent oblivion faculty to polynomial extrapolation with an independent power weighting on past samples. It accounts for memory depth changes required to adapt to sudden variations. Computations are performed on frames of samples hopping at synchronization steps, allowing a low-complexity matrix formulation. The resulting Causal Hopping Oblivious Polynomials (CHOPoly) are described next, evaluated on the case study described in Section 6 and tested in Section 7.

Notations

The first convention abstracts a novel sample indexing framework. It becomes independent of the actual sampling period H and the running time index. This amounts to considering a unit communication step for any signal u(t), i.e. H = 1, and to using a zero-reindexing convention: the last available sample is indexed by 0 (u 0 ), and the previous samples are thus indexed backwards: u -1 , u -2 , u -3 ,. . . This simplification possesses two main traits:

1. local indices hop at each communication step and become independent of the actual timing and communication rate;

2. some intermediate computations may thus be performed only once, reducing the risk of cumulative numerical errors [START_REF] Beckermann | The sensitivity of least squares polynomial approximation[END_REF].

We can either use a recursive formulation with infinite memory or a finite prediction frame. The first option is used, for instance, in adaptive filtering or Kalman estimation. It generally includes oblivion, implemented for instance using a scalar forgetting factor w l , l ∈] -∞, . . . , -1, 0], assigning a decreasing weight, often exponential [START_REF] Gardner | Exponential smoothing: The state of the art -part II[END_REF], to older error samples. We instead consider the second option: a finite frame of the λ last consecutive past samples {u 1-λ , u 2-λ , . . . , u 0 }. Limited-sized buffers are more natural for polynomial forecasting, and they also appear beneficial to reset computations in the case of sudden changes in the data, as discussed in Section 5.2. To emulate a variable memory depth weight, we choose a piece-wise power weighting with order ω ≥ 0. It can be expressed as follows, without λ and ω indices to lighten notations:

w l =            0 if l < 1 -λ , λ + l λ ω 1 -λ ≤ l ≤ 0 . (3)
Their behavior for different powers ω is illustrated in Figure 4. For ω = 0, all samples in the frame are assigned the same importance. The higher the power, the smaller the influence of older samples.

CHOPoly: Type I, symbolic and Type II implementations

Although the forthcoming derivations are elementary, they are rarely exposed with their complexity evaluated. We consider polynomial predictors, performed in the least-squares sense [26, p. 227 sq.]. We denote by P δ,λ,ω the least-squares polynomial predictor of degree δ ∈ N, frame length λ ∈ N * and weighting factor ω. Its stable estimation requires that λ > δ. The polynomial P δ,λ,ω is defined by the vector of length δ+1, with coefficients a δ ; hence u(t) = a δ +a δ-1 t+• • •+a 0 t δ . A causal prediction consists in the estimation of unknown data at a future and relative time τ ≥ 0. Generally τ ∈]0, 1[, i.e., inside the time interval between the last known sample u 0 and the forthcoming communication step. The predicted value at time τ is loosely denoted by u P δ,λ,ω (τ), or simply u(τ) when the context is self-explanatory. We refer to Appendix A.1 for an introductory example on constant weighting, degree-two or parabolic prediction (P 2,λ,0 ). 9 0.8 0.9 In the most generic context, we minimize the least-squares prediction error:

1 λ = 2, ω = 1/8 λ = 2, ω = 1/2 λ = 2, ω = 1 λ = 2, ω = 2 λ = 3, ω = 1/8 λ = 3, ω = 1/2 λ = 3, ω = 1 λ = 3, ω = 2 λ = 5, ω = 0 λ = 5, ω = 1/8 λ = 5, ω = 1/2 λ = 5, ω = 1 λ = 5, ω = 2 λ = 7, ω = 5
e(a δ ) = 0 l=1-λ λ + l λ ω ×        u l -        δ d=0 a d l δ-d               2 . ( 4
)
If one chooses ω = 0, and λ = δ + 1, one exactly recovers Waring-Lagrange polynomials, and the error e(a δ ) is exactly zero. Consequently, CHOPoly encompasses standard extrapolation used in co-simulation (e.g. ZOH, FOH, SOH or higher-order). Choosing λ > δ + 1 induces a form of smoothing, that reduces oscillations. Choosing ω > 0 generally promotes a better extrapolation near the last know sample, making smoothing lag-dependent. Both parameters take into account the possibility that the signal values at communication steps could be imperfect, due to jitter in the actual sampling, noise or round-off errors for instance. As λ ω is a non-null constant, the derivation of (4) with respect to a i yields δ + 1 equations, for each i ∈ {0, . . . , δ}:

0 l=1-λ (λ + l) ω l δ-i u l = δ d=0         0 l=1-λ (λ + l) ω l 2δ-d-i         a d .
We define the sums of powers z d,λ = λ-1 l=0 l d (see Appendix A.1) and the weighted sums of powers z d,λ,ω = λ-1 l=0 (λl) ω l d . They can be rewritten, if ω ∈ N, as:

z d,λ,ω = ω o=0 (-1) o o ω λ ω-o z o+d,λ ,
where o ω denotes the binomial coefficient. The associated weighted moments write, accordingly, m d,λ,ω = λ-1 l=0 (λl) ω l d u -l . We refer to Appendix A.1 for additional information on power sums z d,λ = λ-1 l=0 l d and moments. The generic extrapolation pattern takes the first following matrix form (Type I):

u(τ) = 1 τ • • • τ δ                     z 0,λ,ω -z 1,λ,ω • • • (-1) δ z δ,λ,ω -z 1,λ,ω . . . . . . . . . . . . . . . . . . . . . (-1) δ z δ,λ,ω • • • • • • z 2δ,λ,ω                     -1                  m 0,λ,ω -m 1,λ,ω . . . (-1) δ m δ,λ,ω                  . (5)
We note

τ δ = [1, τ, • • • , τ δ ]
T the vector of τ powers, and

m δ,λ,ω = m 0,λ,ω , -m 1,λ,ω , • • • , (-1) δ m δ,λ,ω T .
The inverse of the (δ + 1) × (δ + 1) Hankel matrix Z δ,λ,ω in (5):

                    z 0,λ,ω -z 1,λ,ω • • • (-1) δ z δ,λ,ω -z 1,λ,ω . . . . . . . . . . . . . . . . . . . . . (-1) δ z δ,λ,ω • • • • • • z 2δ,λ,ω                    
is denoted by Z -δ,λ,ω . Hence, the Type-I formula for CHOPoly writes:

u(τ) = τ T δ Z -δ,λ,ω m δ,λ,ω . (6)
A polynomial of degree δ involves terms of degree 2δ. For long runs of simulation, this could lead to huge values computed from large sample indices. This potentially induces roundoff errors and subsequent issues in matrix inversion. Such matrices, small and theoretically of full rank, become numerically ill-conditioned and rank deficient [START_REF] Kaltofen | Fast estimates of Hankel matrix condition numbers and numeric sparse interpolation[END_REF]. This situation is avoided by the frame-hopping zero-reindexing convention. Moreover, the matrix Z -δ,λ,ω may be computed beforehand, as a look-up table. The compact formulation in (6) thus removes numerical instabilities.

The moment-based formulation conceals more direct symbolic formulae. They are linear in past samples and polynomial in τ. Some examples are given in Appendix A.2. They can be implemented with Ruffini-Horner's rule [START_REF] Peña | On the multivariate Horner scheme[END_REF] for higher extrapolation orders. However, symbolic polynomial evaluations are not always handy to implement.

We note u λ = [u 0 , u -1 , . . . , u 2-λ , u 1-λ ] T . Then (6) can be rewritten in a Type II form. For instance, since τ T 0 = 1, extrapolation with polynomial P 0,λ,1 (A.3) rewrites as a matrix product with Π 0,λ,1 :

u(τ) = 2 λ(λ + 1) λ λ -1 • • • 1 u λ = τ T 0 Π 0,λ,1 u λ , (7)
The general Type-II formula for predictor matrices Π δ,λ,ω of size (δ + 1) × λ implements (6) as:

u(τ) = τ T δ Π δ,λ,ω u λ . (8)
Such a formulation allows an efficient storage of matrices Π λ,δ,ω with different weighting factors. Examples are provided in Table A.5 in Appendix A.4. Estimates of the generic number of elementary operations required for each extrapolated τ are compared in Table A.4 in Appendix A.3. Type I and Type II are roughly quadratic in (δ, λ). Although the rectangular matrix Π δ,λ,ω is larger than Z -δ,λ,ω , since λ > δ, the economy in lazy matrix addressing combined with the number of elementary operations make Type II implementation more practical and efficient.

CHOPatt: Contextual & Hierarchical Ontology of Patterns

Pattern representation and two-level context hierarchy definitions

We now introduce the pattern-based approach, borrowed from common lossless image encoders (Section 3.4.2), by providing a contextual and hierarchical framework for CHOPoly extrapolation. To this aim, the role of contexts is to cover all possible scenarios (an ontology) of signal's evolution for a hybrid dynamical system. For instance, slow and steep variations must be included, the same goes for blocky and smooth signals.

On the first level, we define a functional context that differentiates a set of typical patterns, represented in Figure 5. We are interested in the dynamics of the most recent samples, depicted with links in dotted lines with respect to the past behavior (links in solid lines). It defines six mutually exclusive entities, illustrated by a short name: "flat", "calm", "move", "rest", "take" and "jump". 1.

"Flat" addresses steady signals. "Calm" represents a sufficiently sampled situation, where value increments over time remain below fixed thresholds. "Move" defines a formerly "calm" signal whose novel value changes rapidly, above a pre-defined threshold. "Rest" handles signals previously varying above a threshold and becoming "calm". The "take" context addresses signals with constantly high variations. The "jump" pattern is a "take" with a sign change as for bumps and dips. They are detected in practice by comparing consecutive differences on past samples (d i , cf. (9)), indexed by i, to thresholds (γ i , cf. (10)). Their formal definition is provided in Table 1 and their adaptive selection in Section 5.2.

On the second level, we define a meta-or decisional context that determines if extrapolation would be beneficial or detrimental. Indeed, a major difficulty for hybrid complex systems resides in sharp and fast variations in signal patterns induced by stiffness and discontinuities. We thus detect if the signal can be characterized by a seventh "cliff" pattern, for which functional CHOPoly prediction would fail, by taking into account the importance of the signal's amplitude. This decision context is detected by comparing a ratio ρ of differences based on past samples (12) to a threshold (denoted by Γ, cf. ( 13)). The above seven patterns define a hierarchical context selection, composed of the decisional context embedding the functional context on a second level, summarized in Figure 6.

Functional context selection

As a concrete example, used in the remaining of our work, we propose two simple measures of variation based on the last three samples only: u 0 , u -1 and u -2 (with the zero-reindexing convention from Section 4.2), the last and previous differences:

d 0 = u 0 -u -1 and d -1 = u -1 -u -2 . ( 9
)
Their absolute values are compared with thresholds γ 0 and γ -1 , respectively, defined in [START_REF] Khaled | Context-based polynomial extrapolation and slackened synchronization for fast multi-core simulation using FMI[END_REF]. To build the different contexts, three complementary conditions are defined:

• O if |d i | = 0; • C i if 0 < |d i | ≤ γ i ; • C i if |d i | > γ i .
Table 1 formally defines the six entities from the functional context and presents examples of "default" ω-parametrized CHOPoly families. Since the "flat" context addresses steady signals, a n(ame

) # |d -1 | |d 0 | d -1 .d 0 (δ, λ, ω) f(lat) 0 O O O (0, 1, .) c(alm) 1 C 1 C 2 any (2, 5, .) m(ove) 2 C 1 C 2 any (0, 1, .) r(est) 3 C 1 C 2 any (0, 2, .) t(ake) 4 C 1 C 2 > 0 (1, 3, .) j(ump) 5 C 1 C 2 < 0 (0, 1, .)
Table 1: Functional context: entities definition and associated prediction.

mere ZOH suffices, hence P 0,1,ω . The "calm" context represents smooth signals, in this case, it could be approximated by a quadratic polynomial, for instance, P 2,5,ω . For the "flat" and "jump" contexts, an additional procedure consists in resetting the extrapolation to prevent inaccurate prediction. For example, when context 1 is chosen just after context 4, the quadratic extrapolation P 2,5,ω requires 5 valid samples, whereas the last 3 only are relevant, justifying our finite-length frame option. The choice of thresholds γ 0 and γ -1 is potentially crucial. For instance, fixed values may reveal inefficient under important amplitude or scale variation in signals. Hence, we have chosen to update them adaptively, based on the statistics of a past frame {u 1-Λ , . . . , u -3 } (Λ denotes the maximum frame size).

CHOPatt

With excessively low thresholds, high-order extrapolations would rarely be chosen, losing the benefits of prediction. Overly high thresholds would, in contrast, suffer from any unexpected jump or noise. As contexts are based on backward derivatives, we have used in the simulations presented here the mid-range statistical estimator of their absolute values. This amounts to set:

γ 0 = γ -1 = 1 2 max i∈[1-Λ,...,-3] (|u i -u i+1 |) , (10)
which appeared to be sufficiently robust in our test-cases.

Next communication period: weighting factor and decisional context selection

To further improve our algorithm in [START_REF] Khaled | Context-based polynomial extrapolation and slackened synchronization for fast multi-core simulation using FMI[END_REF] and to decrease even more prediction induced integration errors, the oblivious weighting factor from (3) makes the algorithm more subtly aware of data freshness. Decisions for the next communication period are taken with respect to the new updated input value u 0 , compared to past predictions. It is gauged with two error measures, illustrated in Figure 7. They assess, a posteriori, what would have been the best prediction. First, ∆ worst denotes the worst case scenario without extrapolation, the predicted value being held constant and equal to u -1 :

∆ worst = |u 0 -u -1 |.
Second, we estimate the best prediction pattern in retrospect, obtained by optimizing the free parameter ω in the prediction polynomial defined by the last selected context. We choose here a subset of weighting factors in Ω = 0, 1

A.4 for details):

∆ best = min ω∈Ω |u 0 -ûω -1 |.
The best weighting factor ω best defined in ( 11) is then selected during the next communication interval:

∆ best = |u 0 -ûω best -1 |. (11)
Δ worst Δ best u -1 u 0 û ω1 -1 û ω3 -1 û ω2 -1
Figure 7: Illustration of two error measures: ∆ worst when there is no extrapolation and ∆ best for the best prediction pattern.

Regarding the selection of the decisional context called "cliff", we first define a ratio:

ρ = ∆ best ∆ worst . (12)
It specifies whether there is a sharp fast variation or not by comparing it to a pre-defined threshold and then decides to extrapolate or not. The threshold denoted Γ is defined as follows:

0.7 ≤ Γ < 1. (13)
In fact, when ρ > Γ (e.g. with Γ = 90 %), it means that the input value is only enhanced by 1 -Γ (e.g. 10 %) regarding the "true" value. This is the case when there is a sharp and fast variation or a weird behavior. The decisional context "cliff" is then selected and activated with its associated heuristic polynomial predictor P δ,λ,ω = P 0,1,ω . On the other hand, when ρ ≤ Γ, the conventional functional context table introduced in Table 1 of Section 5.2 is used.

Case study

Engine simulator

In this study, a Spark Ignition (SI) RENAULT F4RT engine has been modeled with 3 gases (air, fuel and burned gas). It is a four-cylinder, in-line Port Fuel Injector (PFI), engine in which the engine displacement is 2000 cm 3 . The combustion is considered as homogeneous. The air path (AP) consists in a turbocharger with a mono-scroll turbine controlled by a waste-gate, an intake throttle and a downstream compressor heat exchanger. This engine is equipped with two Variable Valve Timing (VVT) devices, for intake and exhaust valves, to improve the engine efficiency (performance, fuel and emissions). The maximum power is about 136 kW at 5000 rpm.

The F4RT engine model was developed using the ModEngine library [38]. ModEngine is a Modelica [39] library that allows the modeling of a complete engine with diesel and gasoline combustion models.

Requirements for the ModEngine library were defined and based on the already existing IFP-Engine library. The development of the IFP-Engine library was performed several years ago at "IFP Energies nouvelles" and it is currently used in the AMESim1 tool. ModEngine contains more than 250 sub-models. It has been developed to allow the simulation of a complete virtual engine using a characteristic time-scale based on the crankshaft angle. A variety of elements is available to build representative models for engine components, such as turbocharger, wastegate, gasoline or diesel injectors, valve, air path, Exhaust Gas Recirculation (EGR) loop etc. ModEngine is currently functional in Dymola2 .

The engine model and the split parts were imported into xMOD model integration and virtual experimentation tool [START_REF] Gaïd | Heterogeneous model integration and virtual experimentation using xMOD: Application to hybrid powertrain design and validation[END_REF], using the FMI export features of Dymola. This cyber-physical system has 118 state variables and 312 event indicators (of discontinuities).

Decomposition approach

The partitioning of the engine model is performed by separating the four-cylinder from the air path, then by isolating the cylinders (C i , for i ∈ {1, . . . , 4}) from each other. This kind of splitting allows for the reduction of the number of events acting on each sub-system. In fact, the combustion phase raises most of the events, which are located in the firing cylinder. The solver can process them locally during the combustion cycle of the isolated cylinder, and then enlarge its integration time-step until the next cycle.

From a thermodynamic point of view, the cylinders are weakly coupled, but a mutual data exchange does still exist between them and the air path.

The dynamics of the air path is slow (it produces slowly varying outputs to the cylinders, e.g. temperature) compared to those of the cylinders (they produce fast outputs to the air path, e.g. the global torque). Besides, unlike cylinder outputs, most air path outputs are not a direct function of the air path inputs (they are called Non-Direct Feedthrough (NDF) outputs). This results in choosing the execution order of the split model from the air path to the cylinders (in accordance with the analysis of the behavior of NDF to Direct Feedthrough (DF) in [START_REF] Ben | Distributed real-time simulation of numerical models: application to power-train[END_REF]Chapter 9]).

The model is split into 5 components and governed by a basic controller denoted CTRL. It gathers 91 inputs and 98 outputs.

Tests and results

Tests are performed on a platform with 16 GB RAM and an "Intel Core i7" 64-bit processor, running 4 cores (8 threads) at 2.70 GHz.

Reference simulation

The model validation is based on the observation of some quantities of interest as the pressure, the gas mass fraction, the enthalpy flow rate, the torque, etc. These outputs are computed using LSODAR 3 , a variable time-step solver with a root-finding capability that detects the events occurring during the simulation. It also has the ability to adapt the integration method to the observed system stiffness.

The simulation reference Y ref is built from the integration of the entire engine model, the absolute solver tolerance (tol) being decreased until stable results. The latter are reached for tol = 10 -7 , at the cost of an unacceptable slow simulation speed.

Then, to explore the trade-offs between simulation speed and precision, simulations are run with increasing values of the solver tolerance until reaching a desired relative integration error Er, defined by ( 14)

Er(%) = 100 N . N-1 i=0 Y ref (i) -Y(i) Y ref (i) , (14)
with N the number of saved points during 1 s of simulation. Iterative runs showed that the relative error converge to a desired error (Er ≤ 1 %) for tol = 10 -4 . The single-thread simulation of the whole engine with LSODAR and tol = 10 -4 provides the simulation execution time reference, to which the parallel version is compared. When using the split model, each of its 5 components is assigned to a dedicated core and integrated by LSODAR with tol = 10 -4 .

Automatic detection of fast and sharp variations

Adding a hierarchy of decisional and functional contexts overcomes previous limitations in [START_REF] Khaled | Context-based polynomial extrapolation and slackened synchronization for fast multi-core simulation using FMI[END_REF]. We illustrate this with two signals denoted "Out1" and "Out2". They are built in Matlab/Simulink as shown in Figure 8a. They exhibit different variations over time, as illustrated in Figure 8b. From Figure 9, extrapolation of "Out1" fails around t = 8 s at the sharp variation. Here extrapolation is detrimental since it increases errors instead of minimizing them. To fix this limitation on high jumps, we apply hierarchical context selection to detect the "cliff" context. Using the ratio ρ (12) and comparing it to a threshold Γ (13) equal to 90 %, the improvement of extrapolation at this step is found lesser than 10 %. The "cliff" context is then activated to avoid extrapolation. As a result, the decisional context prevents additional prediction errors and the result of the simulation is improved as it is shown in Figure 10.

Automatic selection of the weighting factor

To determine the best weighting factor regarding error minimization, we first test them separately on "Out1" and "Out2". This first test is quite simple since there is no interaction between blocks. It means that extrapolated signals have no effect on each other. The purpose here is to show the difference between all weighting factors.

Figure 11 shows the absolute error, which is the absolute value of the difference between the reference in Figure 8b built with a small communication step H = 10 µs and signals simulated with a larger communication step H = 100 µs, extrapolated or not. It can be inferred that the higher the weighting factor, the smaller the error.

Table 2 shows the cumulative absolute error on a long simulation run (during 10 s). It confirms that the weighting factor ω = 2 is the best regarding error reduction. In fact, the extrapolation is enhanced from ω = 0 to ω = 2 by reducing the error of prediction by 20.10 % for "Out1" and by 11.39 % for "Out2".

The same experience is now applied on the F4RT engine model and extrapolation with different weighting factors is applied separately on all engine inputs. Figure 12 represents one of the cylinder 1 outputs, "the enthalpy flow rate", for the different extrapolation modes.

We notice more clearly from Figure 13 that for each communication step, there is a different best weighting factor that minimizes the absolute error. Besides, the computation of the cumulative integration error, during a long simulation run, shows that there is no unique best weighting factor.

The weighting factor ω is then chosen dynamically during the simulation as described in Section 5.3. At each communication step, the weighting factor that minimizes the last error is selected and used for the current step.

Thanks to dynamic error evaluation and weighting factor selection, the cumulative integration error is almost always the lowest one for the different outputs as shown in Figures 14a and14b. However, the worst enhancement of the error depends on the output, for instance, it is obtained with ω = 1 2 for the air mass fraction of the cylinder 1 and with ω = 2 for the fuel mass fraction of the cylinder 1. This confirms that for complex coupled systems, there is no unique best weighting factor, hence the necessity and the usefulness of combining different ones.

Besides, for the "burned gas mass fraction" of cylinder 1 (see Figure 15), the dynamic weighting factor selection decreases the error of prediction by 32 % compared to the previous work (with ω = 0, in [START_REF] Khaled | Context-based polynomial extrapolation and slackened synchronization for fast multi-core simulation using FMI[END_REF]) as well as the simulation error by 40 % compared to the non-extrapolated signal.

Regarding now the achievement on the simulation speed-up, Table 3 shows the acceleration compared with the single-threaded reference. Firstly, the speed-up is supra-linear with respect to the number of cores when the model is split into 5 threads integrated in parallel on 5 cores. Indeed, the containment of events detection and handling inside small sub-systems allows for solvers accelerations, enough to over-compensate the multi-threading costs. relative error increase (e.g. +341 % for the fuel density). Thirdly, the combination of model splitting with expanded communication steps (use of H = 250 µs) as well as CHOPtrey allows to keep the same extra speed-up while decreasing the relative error to values close to, or below, those measured with H = 100 µs and ZOH. We can conclude that the enhancement brought with CHOPtrey allows for improved performance on both sides: simulation time and results' accuracy.

Summary and perspectives

The main objective of CHOPtrey is to provide a framework for hybrid dynamical systems cosimulation speed-up [START_REF] Broman | Requirements for hybrid cosimulation[END_REF] overheads. In addition, sustained or even improved simulation precision is obtained without noticeable instability. This work can be extended in different directions. Keeping with data extrapolation, simulated signals can be cleaned from long range trends [START_REF] Ning | Chromatogram baseline estimation and denoising using sparsity (BEADS)[END_REF], to better detect subtle behavioral modifications, and subsequently, adapt detection thresholds. They can be processed on different time grids [START_REF] González | On the effect of multirate co-simulation techniques in the efficiency and accuracy of multibody system dynamics[END_REF] with multi-scale techniques [START_REF] Chaux | Noise covariance properties in dual-tree wavelet decompositions[END_REF]. Acting as local pseudo-derivatives, the latter can decompose signals into morphological components such as polynomial trends, singularities and oscillations. This would allow improvements in context assignment by measuring sharp variations and spurious events with data-relative sparsity metrics [START_REF] Repetti | Euclid in a taxicab: Sparse blind deconvolution with smoothed 1 / 2 regularization[END_REF].

Moreover, the discrimination of cliff behaviors could be further improved by using knowledge of the plant model, allowing to discard out-of-bound values, e.g., non-negative variables.

Finally, simulation results suggest that widening the communication steps is an important source for integration acceleration. Beyond equidistant communication grids, adaptive, context and/or error-based closed-loop control of communication steps [START_REF] Schierz | Co-simulation with communication step size control in an FMI compatible master algorithm[END_REF][START_REF] Sadjina | Energy conservation and power bonds in co-simulations: Non-iterative adaptive step size control and error estimation[END_REF] is a promising research objective.

The system in (A.1) rewrites:

          m 0,λ -m 1,λ m 2,λ           =           z 0,λ -z 1,λ z 2,λ -z 1,λ z 2,λ -z 3,λ z 2,λ -z 3,λ z 4,λ           ×           a 2 a 1 a 0           or m = Z 2,λ × a .
We now want to find the value predicted by P 2,λ,0 at time τ. Let τ 2 = [1, τ, τ 2 ] T be the vector of τ powers. Then u(τ) is equal to a 2 + a 1 τ + a 0 τ 2 = τ T 2 × a 2 . This system might be solved with standard pseudo-inverse techniques [START_REF] Greville | Some applications of the pseudoinverse of a matrix[END_REF], by premultiplying by the transpose of Z 2,λ . This is not required, as Z 2,λ is always invertible, provided that λ > δ. Its inverse is denoted Z -2,λ . It thus does not need to be updated in real-time. It may be computed offline, numerically or even symbolically. Hence: u(τ) = τ T 2 × Z -2,λ × m 2,λ . The vector τ 2 and Z -2,λ are fixed, and the product τ T 2 × Z -2,λ may be stored at once. Thus, for each prediction, the only computations required are the update of vector m 2,λ and its product with the aforementioned stored matrix.

ω Π 0,2,ω Π 1,3,ω Π 2,5,ω 0 1 2 1 1 1 6        5 2 -1 3 0 -3        1 70             62
                  
0.90 0.23 -0.09 -0.12 0.07 0.80 -0.24 -0.58 -0.32 0.34 0.15 -0.09 -0.14 -0.05 0.13 A.5: Type II matrices Π 0,2,ω , Π 1,3,ω and Π 2,5,ω , with integer and rational weighting powers ω ∈ 0, 1 8 , 1 4 , 1 2 , 1, 2 .

           

Figure 2 :

 2 Figure 2: Σ split into Σ 1 and Σ 2 for parallel mono-rate simulation.

Figure 3 :

 3 Figure 3: Parallel multi-rate simulation with an external model.

Figure 4 :

 4 Figure 4: Effect of the weighting power ω with different choices of λ on the memory depth.

Figure 5 :

 5 Figure 5: Pattern representation for the functional context in Table1.

Figure 6 :

 6 Figure 6: Hierarchical and functional context selection flowchart.

  (a) Construction of Out1 and Out2 in Matlab/Simulink.

  Illustration of Out1 and Out2.

Figure 8 :

 8 Figure 8: A test sample.

Figure 9 :

 9 Figure 9: Failure of the old extrapolation.

Figure 10 :

 10 Figure 10: Success of the new extrapolation.

Figure 12 :Figure 13 :

 1213 Figure 12: Enthalpy flow rate output of cylinder 1.

  Secondly, model splitting combined with enlarged communication steps, from H = 100 µs to H = 250 µs, allows around +12.50 % extra speed-up. Unfortunately, this extra speed-up is obtained at the cost of

Figure 14 :

 14 Figure 14: Absolute error in cylinder 1.

Figure 15 :

 15 Figure 15: Absolute error of the burned gas mass fraction of cylinder 1.

[2] n 2 Σ t s 1 t s 1 +1 t n 1 t n 1 +1 t s 1 +2 t n 2 t n 2 +1 t s 2 t s 2 +1 t s 2 +2

[2] n 2 Σ t s 1 t s 1 +2 t n 1 t n 1 +1 t s 1 +4 t n 2 t n 2 +1 t s 2 t s 2 +1 t s 2 +2

Communication step H [ex] t s ex t s ex +2 t s ex +4 t s ex +1 t s ex +3 Exchange 1 Exchange 2 Exchange 3 Init

Table 2 :

 2 Figure 11: Behavior of the absolute error at each communication step. Cumulative absolute error during 10 s and relative improvement.

	Absolute error for Out1	0 0.05 0.1 0.15	2.14	2.16	2.18	2.2	2.22		2.24	2.26	2.28	2.3 ZOH ω = 0 ω = 1/8 ω = 1/4 ω = 1/2 ω = 1 ω = 2	2.32
	Absolute error for Out2	0 0.05 0.1	2.14	2.16	2.18	2.2	time (s) 2.22	2.24	2.26	2.28	2.3 ZOH ω = 0 ω = 1/8 ω = 1/4 ω = 1/2 ω = 1 ω = 2	2.32
							time (s)
				Type		Error			Improvement
					Out1 Out2		Out1		Out2
				ZOH	0.572 0.399		-		-
			ω = 0	0.204 0.158	64.34 % 60.40 %
			ω = 1/8 0.201 0.157	64.86 % 60.65 %
			ω = 1/4 0.198 0.155	65.38 % 61.15 %
			ω = 1/2 0.193 0.153	66.27 % 61.65 %
			ω = 1	0.182 0.148	68.18 % 62.91 %
			ω = 2	0.163 0.140	71.50 % 64.91 %

Table 3 :

 3 based on extrapolation. Cheaper slackened synchronization intervals are allowed by a combination of prediction and multi-level context selection. It aims at reaching real-time simulation while preserving result accuracy.It is implemented in combination with model splitting and parallel simulation on a hybrid dynamical engine model. It results in effective simulation speed-up with negligible computational 22 CHOPtrey performance: speed-up versus accuracy. The speed-up factor is compared with single-threaded reference. The relative error variation is compared with ZOH at 100 µs.

	Communication Prediction Speed-up	Relative error variation
	time		factor	Burned gas density Fuel density
	100 µs	ZOH	8.9	-	-
	250 µs	ZOH CHOPtrey	10.01 10.07	+7 % -26 %	+341 % +21 %

www.lmsintl.com/imagine-amesim-1-d-multi-domain-system-simulation

www.3ds.com/products/catia/portfolio/dymola

Short for Livermore Solver for Ordinary Differential equations with Automatic method switching for stiff and nonstiff problems, and with Root-finding[START_REF] Hindmarsh | Algorithms and software for Ordinary Differential Equations and Differential-Algebraic Equations, part II: Higher-order methods and software packages[END_REF].

Model Driven Physical Systems Operation.

Acknowledgment

This work was supported by the ITEA2 project MODRIO 4 , and funded in part by the "Direction Générale des Entreprises" of the French Ministry of Industry. The authors thank Bruno Lety, Nicolas Pernet and Yann Creff (IFP Energies nouvelles) for proofreading and discussions.

Author biographies

Abir Ben Khaled-El Feki is a research engineer at IFP Energies nouvelles, Rueil-Malmaison, France. Her current research interests include real-time and multi-core simulation of complex systems using approaches as model splitting, refined co-simulation scheduling and data extrapolation.

Laurent Duval is a research engineer in signal processing and image analysis at IFP Energies nouvelles, Rueil-Malmaison, France. His research interests emphasize on multi-rate filtering, filter bank and wavelet techniques, statistical processing, optimization and their applications to signal detection, restoration and data compression. He holds currently the position of project manager at the Mechatronics, Computer Science and Applied Mathematics Division of IFPEN, the French national research center in energy and transportation. Prior to this, he spent one year as a postdoctoral researcher at the Networked Control Systems project-team of INRIA Rhône-Alpes. His research interests include multi-domain simulation, real-time computing and control systems.

Appendix A.1. Toy parabolic CHOPoly extrapolation formulae

We estimate the best fitting parabola (i.e. δ = 2) with a uniform weighting (ω = 0): u(t) = a δ + a δ-1 t + a δ-2 t 2 to approximate the set of discrete samples {u 1-λ , u 2-λ , . . . , u 0 }. We consider here "uniform" weighting w l = 1 for 1 -λ ≤ l ≤ 0 (i.e. with weighting factor ω = 0). The prediction polynomial P 2,λ,0 is defined by the vector of polynomial coefficients: a 2 = [a 2 , a 1 , a 0 ] T .These coefficients are determined, in the least-squares sense [START_REF] Stigler | Gauss and the invention of least squares[END_REF], by minimizing the squared or quadratic prediction error (4):

Here indices l are non-positive, i.e. between 1 -λ and 0. The minimum error is obtained by solving the following system of equations (zeroing the derivatives with respect to each of the free variables a i ):

System (A.1) may be rewritten as:

Closed-form expressions exist for the sum of powers z d,λ , involving Bernoulli sequences [START_REF] De Bruyn | Formulas for 1 + 2 p + 3 p + . . . + n p[END_REF].

For instance, up to the 4 th power:

• z 0,λ = λ;

• z 1,λ = (λ -1)λ/2;

• z 2,λ = (λ -1)λ(2λ -1)/6;

Let m d,λ = m d,λ,0 = λ-1 l=0 l d u -l (here indices l are positive) denote the d th moment 5 of the samples u i , and m 2,λ the vector of moments [m 0,λ , -m 1,λ , m 2,λ ] T . We now form the Hankel matrix Z 2,λ of sums of powers (depending on δ = 2 and λ):

.

5 Definition: the d th moment of a real function f (t) about a constant c is usually defined as:

. The m d s may be interpreted as discrete versions of one-sided moments about t = 0 of the discrete function u(t); alternatively cf. definitions for z d,λ -the moments are sort of weighted (by u l ) sum-of-powers.

Appendix A.2. CHOPoly symbolic formulation

When only one polynomial predictor is required, actual computations do not require genuine matrix calculus, especially for small degrees δ. With δ = 0 and ω = 0 (or P 0,λ,0 ), one easily sees that:

that is, the running average of past frame values. It reduces to standard ZOH, u(τ) = u 0 , when λ = 1. For δ = 0 and ω = 1, one gets a weighted average giving more importance to the most recent samples:

With λ = 2, δ = 1 and ω = 0, P 1,2,1 (τ) = u 0 + (u 0u -1 )τ yields the simplest 2-point linear prediction or standard FOH. P 1,3,0 yields the simple estimator form:

For P 2,5,1 , we tediously get:

or with the Ruffini-Horner's method for polynomial evaluation, to slightly reduce the number of operations:

One easily remarks that, when the weighting exponent ω is an integer, prediction polynomials have rational coefficients, which limits floating-point round-off errors, especially when prediction times and variables (for instance quantized ones) are integer or rational as well.

Appendix A.3. CHOPoly Type I and II computational complexity The computational complexity of a single extrapolation is given in Table A. [START_REF] Sjölund | Towards efficient distributed simulation in Modelica using Transmission Line Modeling[END_REF]. It is evaluated by a number of elementary operations. These estimates are only meant to provide rough intuitions and guidelines on the actual implemented complexity. Their expressions for Type I (6) and Type II (8) Causal Hopping Oblivious Polynomials implementations are given terms of δ, λ (and ω). Type I direct implementation is not efficient: unoptimized moment computations yield a cubic complexity in (δ, λ, ω). Recurring results can be evaluated using call-by-need or lazy evaluation. For m d,λ,ω , the factors (λl) ω yield λ -2 powers (the λ -1 adds are unnecessary), since 0 ω and 1 ω are direct. For 2 ≤ l ≤ λ, the l d values are gathered in a (δ + 1) × λ array, with

2 -1 products. The (λl) ω l d terms can be stored in a (δ + 1) × λ matrix, involving (λ -2)(δ -1) products. These estimates are stored in the top-half of Table A.4, above the dashed line.

Upon hopping frame update and τ determination, using the simplification τ d+1 = ττ d , τ δ requires δ -1 products. The evaluation of the weighted moments entails only δ(λ -1) adds, and (λ-2)δ+1 products, since the lazy matrix storing (λ-l) ω l d contains some zeroes and ones (Table A.4, bottom-half). Both Type I and II are thus roughly quadratic in (δ, λ). Both are competitive with respect to Waring-Lagrange extrapolation, which requires O δ 2 or O(δ) with δ = λ + 1 depending on the implementation.

If we precisely compute the operation excess from Type I to II, we obtain Ξ(δ, λ) = 2δ 2 + δ + 3-2λ. For the parameters given in Appendix A.4, we have for instance Ξ(0, 2) = -1, Ξ(1, 3) = 0, Ξ(2, 5) = 3.


Question: What kind of symmetric encryption is the most appropriate?

Prior HE-friendly ciphers

Aim:

Minimize the multiplicative depth of the decryption function.

Concrete proposals:

•

Z G x 1 x t F F F • • • F z 1 z 2 z 3 • • • z t keystream = offline online m ⊕ m ⊕ keystream HE pk (k) IV G x 1 x t C F C F C F • • • C F HE pk (keystream) C ⊕ HE pk (m)
→ Minimize the multiplicative depth of F .

Instantiation with a counter

Expansion function G:

G(IV ) = (IV, IV 1, IV 2, . . . , IV (t -1))
Why not use for F a block cipher? security limited to 2 n/2 where n is the block size. → strong limitation for lightweight ciphers with n = 64 or 32. • No strong limitation of the size of the internal state.

Low

Trivium and Kreyvium two low-depth stream ciphers

recommended by the eSTREAM project

• transition function with degree 2

• key size = 80 bits

• IV size = 80 bits

• initialization = 1152 blank rounds

Multiplicative depth of Trivium

The keystream length which can be produced with a circuit of depth d, d ≥ 4, is

282 ×     d 3     +                      81 if d ≡ 0 mod 3 160 if d ≡ 1 mod 3 269 if d ≡ 2 mod 3

Multiplicative depth of Kreyvium

The keystream length which can be produced with a circuit of depth d, d ≥ 4, is

282 ×     d 3     +                      70 if d ≡

  Optimized implementations of AES [Gentry Halevi Smart 12][Cheon et al. 13] [Döroz Hu Sunar 14] • Lightweight block ciphers: SIMON [Lepoint Naehrig 14], PRINCE [Doröz et al.14] • Dedicated block cipher: Low-MC [Albrecht et al.

•

  We need a transition function Φ with a low multiplicative depth.

125 bits Experimental results Using HElib (on one core of a server with 4 x AMD Opteron 6172 processors)

	Using the Fan-Vercauteren scheme on a 48-core server
	0 mod 3 149 if d ≡ 1 mod 3 258 if d ≡ 2 mod 3 (sec.) (bits/min) 13 600 3650 1341 20 720 11380 516 13 682 3987 1272 20 480 12451 287 13 682 3369 3109 security level N used depth throughput (bits/min) Speed gain → 11 bits less than with Trivium • At depth 12, 46 bits Trivium-13 80 136 Kreyvium-13 128 125 LowMC-128 ? ≤ 118 256 1 core 48 cores Trivium-13 80 136 13 9.2 240 × 26.2 20 3.4 106 × 31.0 Kreyvium-13 128 125 13 5.7 151 × 26.5 20 2.2 77 × 34.0 LowMC-128 ? ≤ 118 256 14 10.0 90 × 9.0 • At depth 13, security level N used depth #slots latency throughput 20 480 9977 739 21 4.6 47 × 10.2

Conditional differential cryptanalysis [Knellwolf et al. 11] even in the weak-key setting, 64 bits can never be set to 0

Conclusions

• IV-based stream ciphers are the most appropriate ciphers.

• Good performances can be obtained with firmly-established symmetric ciphers.

Open question.

How many multiplicative levels are necessary to achieve a reasonable security level? → (impractical) approach based on discrete-log achieving a multiplicative depth of ( log κ + 1) for κ-bit security.


Hlib Mykhailenko
email: hlib.mykhailenko@inria.fr

Fabrice Huet
email: fabrice.huet@unice.fr

Giovanni Neglia
email: giovanni.neglia@inria.fr

Comparison of Edge Partitioners for Graph Processing

Keywords: vertex-cut, edge partitioning, GraphX, Spark

Deploying graph on a cluster requires its partitioning into a number of subgraphs, and assigning them to different machines. Two partitioning approaches have been proposed: vertex partitioning and edge partitioning. In the edge partitioning approach edges are allocated to partitions. Recent studies show that, for power-law graphs, edge partitioning is more effective than vertex partitioning. In this paper we provide an overview of existing edge partitioning algorithms. However, based only on published work, we cannot draw a clear conclusion about the relative performances of these partitioners. For this reason, we compare all the edge partitioners currently available for GraphX. Our preliminary results suggest that Hybrid-Cut partitioner provides the best performance.

Introduction

There are many frameworks for distributed graph processing, such as GraphX [START_REF] Xin | Graphx: A resilient distributed graph system on spark[END_REF], Pregel [START_REF] Malewicz | Pregel: a system for large-scale graph processing[END_REF], PowerGraph [START_REF] Gonzalez | Powergraph: Distributed graph-parallel computation on natural graphs[END_REF], GraphLab [START_REF] Low | Distributed GraphLab: a framework for machine learning and data mining in the cloud[END_REF], GraphBuilder [START_REF] Jain | Graphbuilder: scalable graph ETL framework[END_REF], Giraph [START_REF] Avery | Giraph: Large-scale graph processing infrastructure on hadoop[END_REF], etc. These frameworks require the graph to be divided in a number of different partitions, each assigned to a different machine of the computational cluster. Different partitioning algorithms are used to this purpose and in the rest of this paper we refer to them as partitioners.

A common approach to partition a graph is to assign a subset of vertices to each partition, and then as a consequences an edge may be cut when its vertices belong to different partitions. This approach is called vertex partitioning or also edge-cut partitioning.

Recently edge partitioning approach (also called vertex-cut partitioning) was proposed in [START_REF] Gonzalez | Powergraph: Distributed graph-parallel computation on natural graphs[END_REF]. In this case the partitioner assigns edges to partitions, and then a vertex may be cut into several pieces (from 2 till N pieces, where N is the number of partitions), where the number of vertex pieces equals the number of different partitions its edges are assigned. The presence of cut vertices obliges the different partitions to exchange periodically information during the computation to keep aligned the values for the vertices cut. The recent study [START_REF] Bourse | Balanced graph edge partition[END_REF] advocates that edge partitioning is better for power-law graphs (most of the real-world graphs are power-law graphs). Several graph processing frameworks rely on edge partitioning, such as GraphX [START_REF] Xin | Graphx: A resilient distributed graph system on spark[END_REF], PowerGraph [START_REF] Gonzalez | Powergraph: Distributed graph-parallel computation on natural graphs[END_REF], and Graph-Builder [START_REF] Jain | Graphbuilder: scalable graph ETL framework[END_REF]. However, there are very few studies dedicated to edge partitioners.

The first contribution of our paper is to provide an overview of all the edge partitioners we are aware of. Beside describing the specific algorithms proposed, we focus on existing comparisons of their relative performance. It is not easy to draw conclusions about which are the best partitioners, because research papers often consider a limited subset of partitioners, different performance metrics and different computational frameworks. For this reason a second contribution of this paper is a preliminary evaluation of all the edge partitioners currently implemented for GraphX, one of the most promising graph processing frameworks. Our current results suggest that Hybrid-Cut is probably the best choice, achieving significant reduction of the execution time with limited time required to partition the graph.

The paper is organized as follows. In Section 2 we present the notation used across the paper and we describe partitioning metrics. Section 3 presents all existing edge partitioning algorithms and is followed by a discussion in Section 4 about what can be concluded from the literature. Our experiments with GraphX partitioners are described in Section 5. Finally, Section 6 presents our conclusions.

Metrics

Let us denote the input directed graph as G = (V, E), where V is a set of vertices and E is a set of edges. An edge partitioner splits the set of edges into N disjoint subsets, E 1 , E 2 , . . . , E N . With some abuse of terminology, we can say that a vertex belongs to the partition E i if it is the source or the destination of an edge in E i . Using this convention, an edge can only belong to one partition, but a vertex is at least in one partition and at most in N partitions. Let V (E i ) denote the set of vertices that belong to partition E i . The vertices that appear in more than one partition are called frontier vertices. Each frontier vertex has thus been cut at least once. Let F (E i ) denote the set of frontier vertices that are inside partition E i and F (E i ) the set of vertices in partition E i that are not cut. Hence we have F (

E i ) = V (E i ) -F (E i ). N in (v)
is the set of source vertices for the edges incoming to vertex v.

There are different metrics which can help measure the quality of a partitioning. The most straightforward ones are execution metrics like the partitioning time, i.e. the time to partition the graph, the number of rounds performed by the partitioner to converge (for iterative partitioners), the execution time of a particular graph algorithm once the graph is partitioned, the network communication overhead, i.e. the amount of messages or bytes exchanged during the partitioning or graph algorithm execution. A limitation of execution metrics is that they are tightly coupled to specific applications and execution environments, making them unsuitable for general comparisons. Moreover, they can be costly to evaluate or measure.

For this reason, specific partitioning metrics have been proposed. These metrics try to evaluate the quality of the partition produced, assuming that certain characteristics of the partition will affect the final execution metrics. In the remaining of this section we list various partitioning metrics that are often used in the literature.

Partitioning metrics

Partitioning metrics can be split into two groups. In the first group there are the metrics that quantify how homogeneous the partitions' sizes are. The underlying idea is that if one partition is much larger that the others, the computational load on the corresponding machine is higher and then this machine can slow down the whole execution. The metrics in the second group quantify how much overlap there is among the different partitions, i.e. how many vertices appear in multiple partitions. This overlap is a reasonable proxy for the amount of inter-machine communication that will be required to merge the results of the local computations. The first two metrics below belong to the first group and we call them balance metrics, all the other ones are communication metrics.

Balance

It is measured as the ratio of the maximum number of edges in a partition to the average number of edges across all the partitions.

Balance = max i=1,...N |Ei| |E|/N .

Standard deviation of partition size (STD)

It is the normalized standard deviation of the number of edges in each partition.

STD = N i=1 |Ei| |E|/N -1 2 1 N 3 2.1.

Replication factor

It is the ratio of the number of vertices in all the partitions to the number of vertices in the original graph. It measures the overhead, in terms of vertices, induced by the partitioning.

Replication Factor = N i=1 |V (E i )| 1 |V| 2.1.

Communication cost

It is defined as the total number of frontier vertices.

Communication cost = N i=1 |F (E i )|

Vertex-cut

This metric measures how many times vertices were cut. For example, if a vertex is cut in 5 pieces, then its contribution to this metric is 4.

Vertex-cut = N i=1 F (E i ) -|V | + N i=1 F (E i ) 2.1.6 Normalized vertex-cut
It is a ratio of the vertex-cut metric of the partitioned graph to the expected vertex-cut of a randomly partitioned graph.

Expansion

Originally expansion was introduced in [START_REF] Kannan | On clusterings: Good, bad and spectral[END_REF][START_REF] Flake | Graph clustering and minimum cut trees[END_REF] in the context of vertex partitioning. In [START_REF] Kim | SBV-Cut: Vertex-cut based graph partitioning using structural balance vertices[END_REF] expansion was adapted to edge partitioning approach. It measures the relative contribution of a partition to the number of frontier vertices.

Expansion = max i=1,...N |F (Ei)| |V (Ei)|

Modularity

As expansion, modularity was proposed in [START_REF] Kannan | On clusterings: Good, bad and spectral[END_REF][START_REF] Flake | Graph clustering and minimum cut trees[END_REF] and later adapted to edge partitioning approach in [START_REF] Kim | SBV-Cut: Vertex-cut based graph partitioning using structural balance vertices[END_REF] as follows:

1 Modularity = N i=1   |V (Ei)| |V| - j =i |F (Ei)∩F (Ej )| |V| 2  
Intuitively, the higher the modularity is, the smaller is the fraction of vertices shared among different partitions.

Algorithms

In this section we will first provide a classification of the partitioners. We already mentioned the main distinction between edge partitioning and vertex partitioning, but the focus of this paper is on the former type.

Then, we can classify partitions by the amount of information they require; they can be online or offline [START_REF] Jain | Graphbuilder: scalable graph ETL framework[END_REF]. Online partitioners decide where to assign an edge on the basis of local information (e.g. about its vertices, and their degrees). For this reason they can be easily distributed. On the contrary, in offline partitioners (like METIS [START_REF] Karypis | METIS-unstructured graph partitioning and sparse matrix ordering system[END_REF]) the choice depends in general on the whole graph and multiple iterations may be needed to produce the final partitioning.

Finally, some partitioners can refine the partitioning during the execution of a graph processing algorithm (e.g. partitioners in Giraph [START_REF] Avery | Giraph: Large-scale graph processing infrastructure on hadoop[END_REF]). For example, after several iterations of the PageRank algorithm, the graph can be re-partitioned, based on the current execution time of the PageRank algorithm itself.

Most of the non-iterative algorithms have O(|E|) time complexity, and iterative algorithms usually have O(k|E|) time complexity, where k is the number of iterations of the partitioner.

In the rest of the section we introduce 17 partitioners. For them we propose the following classification based on their approach to partition. Table 2 provides a list of references where the following algorithms were first described or implemented.

Random assignment

Random assignment approach includes partitioners that randomly assign edges using a hash function based on some value of the edge or of its vertices. Sometimes, the input value of the hash function is not specified (e.g. [START_REF] Rahimian | Distributed Vertex-Cut Partitioning[END_REF]). Because of the law of large numbers, the random partitions will have similar sizes and then these partitioners achieve good balance.

RandomVertexCut

It assigns each edge using a hash function based on the pair of values source vertex id and destination vertex id.

CanonicalRandomVertexCut

It works is similar to previous one but first it orders the two ids.

EdgePartition1D

It was implemented in GraphX [START_REF] Xin | Graphx: A resilient distributed graph system on spark[END_REF] and it uses only the source vertex id in the hash function.

Randomized Bipartite-cut (BiCut)

This partitioner relies on the idea that random assigning the vertices of one of the two independent sets of a bipartite graph may not introduce any replicas. In particular BiCut selects the largest set of independent vertices, and then splits it randomly into N subsets. Then a partition is created from all the edges connected to the corresponding subset.

Segmenting the hash space

This approach complements random assignment with a segmentation of the hash space using some geometrical forms as grids, torus, etc. It maintains the good balance of the previous partitioners, while limiting at the same time the maximum replication factor.

Grid-based Constrained Random Vertex-cuts (Grid)

It partitions edges in a shard-grid G by using a simple hash function. For example the EdgePartition2D partitioner in [START_REF] Xin | Graphx: A resilient distributed graph system on spark[END_REF] maps the source vertex to a column of G and the destination vertex to a row of G. The edge is then placed in the shard at the column-row intersection and the shards are finally assigned to partitions in a round robin fashion. This partitioner guarantees that the replication factor is upper-bounded by 2 √ n -1. A variant of this algorithm is proposed in [START_REF] Jain | Graphbuilder: scalable graph ETL framework[END_REF], where two different column-row pairs are selected for the source and the destination and then the edge is randomly assigned to one of the shard belonging to both the pairs.

Torus-based Constrained Random Vertex-cuts (Torus)

It is similar to the Grid partitioner considered in [START_REF] Jain | Graphbuilder: scalable graph ETL framework[END_REF] but it relies on a 2D torus. Each vertex is mapped to one column and to 1 2 R + 1 shards of a given row, where R is the number of shards in a row. Then, as in previous partitioner, it considers the shards at the intersection of the two sets identified for the two vertices, and it randomly selects one of them. In this way, the replication factor has an upper bound equal to 1.5

√ n + 1.

Greedy approach

These partitioners assign edges greedily among the partitions, minimizing the current replication factor at each step.

Greedy Vertex-Cuts

To place the i th edge, this partitioner considers where the previous i -1 edges have been assigned. Basically, it tries to place an edge in a partition which already contains the source and the destination of this edge. If it is not possible, then it tries to allocate the edge to a partition which contains at least one vertex of this edge. If no such a partition can be found, the edge is assigned to the partition that is currently the smallest one. The authors also described a distributed implementation called Coordinated Greedy Vertex-Cuts (Greedy-Coordinated in what follows). This implementation requires communication about the different instances of the partioner. The authors also introduced a relaxed version called Oblivious Greedy Vertex-Cuts (Greedy-Oblivious in what follows), where no communication is required because each instance of the partitioner greedily assign an edge as described above but only considering its previous choices.

Grid-based Constrained Greedy Vertex-cuts (Grid Greedy)

As the Grid partitioner proposed in [START_REF] Jain | Graphbuilder: scalable graph ETL framework[END_REF] it relies on a shard grid, but the final shard is selected (from intersection) using the same greedy criterium in 3.3.1.

Torus-based Constrained Greedy Vertex-cuts (Torus Greedy)

The algorithm modifies the Torus algorithm in the same spirit the Grid Greedy partitioner modifies the Grid one.

Distributed Funding-based Edge Partitioning (DFEP)

This partitioner was proposed in [START_REF] Guerrieri | Dfep: Distributed funding-based edge partitioning[END_REF] and implemented for both Hadoop [START_REF] Bialecki | Hadoop: a framework for running applications on large clusters built of commodity hardware[END_REF] and Spark [START_REF] Zaharia | Resilient distributed datasets: A fault-tolerant abstraction for in-memory cluster computing[END_REF]. In the DFEP algorithm each partition receives initially a randomly assigned vertex and then tries to progressively grow by including all the edges of the vertices currently in the partition that are not yet assigned. In order to avoid a dishomogeneous growth of the partitions, a virtual currency is introduced and each partition receives an initial funding that can be used to bid on the edges. Periodically, each partition receives additional funding inversely proportional to the number of edges it has.

Cut hubs

Most of the real-world graphs are power-law graphs, where a relatively small percentage of nodes (hubs) concentrate most of the edges. This approach moves then from the observation that in any case hubs need to be cut to maintain balance among the partitions. These partitioners prefer then to cut hubs, trying to spare the large share of nodes that have a small degree.

Hybrid-Cut

The Hybrid-cut [START_REF] Chen | Powerlyra: Differentiated graph computation and partitioning on skewed graphs[END_REF] partitioner considers that a vertex is not a hub if its indegree is below a given threshold. ISn such case all the incoming edges are partitioned based on the hash value of this vertex and then they are placed in the same partition. Otherwise, all incoming edges are partitioned based on their source vertices. The algorithm was implemented for GraphX [START_REF]Hybrid-Cut and HybridCutPlus partitioners[END_REF].

Ginger

As Hybrid-cut this partitioner [START_REF] Chen | Powerlyra: Differentiated graph computation and partitioning on skewed graphs[END_REF] allocates the edges incoming to a no-hub vertex (v) to the same partition. The difference it that the partition P i is selected on the basis of the following heuristic:

i = argmax i=1,...N |N in (v) ∩ V (E i )| -1 2 (|V (E i )| + |V | |E| |E i |) .
The underlying idea is that the partitioner selects the partitions where there are already neighbors of node v, as far as this partition is not too large.

HybridCutPlus

It works as Hybrid-Cut whenever possible (i.e. if one vertex of an edge is a hub and another vertex is not a hub), otherwise it works as a Grid partitioner. It was implemented for GraphX [START_REF]Hybrid-Cut and HybridCutPlus partitioners[END_REF].

Aweto

It is a greedy implementation of BiCut using a heuristic similar to the one in the Ginger partitioner.

Iterative approach

While the previous partitioners assign once and for all a link to a partition, other partitioners may iterate multiple times on the partitioning, reassigning an edge to a different partition.

DFEPC

It is an iterative variant of DFEP that compensates for the possible negative effect of starting from an initial vertex that is poorly connected to the rest of the graph. DFEP is changed as follows: a partition whose size is by a given factor smaller than the current average partition size can also bid for edges that were already bought by other partitions, leading to their reassignments.

JA-BE-JA-VC

It was proposed in [START_REF] Rahimian | Distributed Vertex-Cut Partitioning[END_REF], and it is a vertex-cut version of an existing edge-cut partitioner JA-BE-JA [START_REF] Rahimian | Ja-be-ja: A distributed algorithm for balanced graph partitioning[END_REF]. The JA-BE-JA-VC algorithm starts by randomly assigning edges to partitions. Then each iteration of the algorithm is executed on each vertex in parallel. Two vertices are allowed to exchange the partition their edges belong to, provided that the exchange leads to a better cut of the graph. The potential benefit of an exchange is evaluated through two different heuristics (referred to as Dominant Color or Edge Utility). Simulated annealing can also be used to escape local minima.

In [START_REF] Rahimian | Distributed Vertex-Cut Partitioning[END_REF] the authors also consider a simple way to adapt JA-BE-JA [START_REF] Rahimian | Ja-be-ja: A distributed algorithm for balanced graph partitioning[END_REF] to produce an edge partitioning. First, JA-BE-JA is executed to produce vertex partitions. Then the edges that are cut are randomly assigned to one of the two partitions their vertices belong to.

Discussion

Table 1 shows which partitioners were directly compared and which metrics were considered in the comparison ("I" denotes execution metrics, while "II" denotes partitioning metrics). Bold fonts indicate that the comparison was performed by us (see Sec. 5). The table shows how many direct comparisons are missing (e.g. Hybrid-Cut was not compared to Torus, JA-BE-JA-VC, etc.) and our experiments in this paper contribute to provide a more complete picture. Even when a direct comparison is available, results are not necessarily conclusive. For example, in [START_REF] Rahimian | Distributed Vertex-Cut Partitioning[END_REF] the authors show that JA-BE-JA-VC outperforms DEFP, DFEPC, and the modified version of JA-BE-JA in terms of STD. However in terms of normalized vertex-cut, the modified version of JA-BE-JA is the best.

Table 2 indicates which specific metrics were considered in each paper. We can observe that in many cases execution metrics have not been considered due to their cost in terms of computation time. In particular, the lack of information about the partitioning time is problematic because in our experiments we have observed that it can vary by more than one order of magnitude across partitioners and contribute the most to the total execution time. The table also shows that it is difficult to perform an indirect comparison among partitioners using results from different papers, because there is often no common set of metrics considered across them.

We showed that it is hard to reach any conclusion regarding the benefits of existing partitioners for multiple reasons. First, not all partitioners have been compared. Second, execution metrics are not always provided and are tied to a specific computing environment. Finally, there is no study which links partitioning metrics to execution metrics, but for [START_REF] Mykhailenko | Which Metrics for Vertex-Cut Partitioning?[END_REF]. Motivated by these considerations, we decided to conduct experiments using the GraphX framework. Our results partially complete the pairwise comparison in Table 1 and are presented in the next section.

Experiments

We conducted our experiments on a computational cluster of nodes with dual-Xeon E5-2680@2.80GHz, 192GB RAM and 20 cores. We used Spark version 1.4.0 in standalone cluster mode. Our cluster configuration had eleven machines, one for master and ten for executors. We configured 4 Spark properties: spark.executor.memory (4 GB), spark.driver.memory (10 GB), spark.cores.max [START_REF] Guerrieri | Dfep: Distributed funding-based edge partitioning[END_REF] and spark.local.dir (set to a local directory).

We measured replication factor, communication cost, balance, STD as partitioning metrics, and the partitioning time and the execution time of Connected Components and PageRank (10 iterations) algorithms as execution metrics. As input we used the undirected com-Youtube graph (1,134,890 vertices/2,987,624 edges) from the SNAP project [START_REF]Stanford Large Network Dataset Collection[END_REF].

As a set of partitioners, we used all the GraphX built-in partitioners, i.e. RandomVertexCut, CanonicalRandomVertexCut, EdgePartition1D, Grid, the partitioner DFEP (whose code was gently made available by the authors of [START_REF] Guerrieri | Dfep: Distributed funding-based edge partitioning[END_REF]),

  Execution time of algorithms.

Figure 1 :

 1 Figure 1: Results for com-Youtube graph (CRVC -CanonicalRandomVertexCut; EP1D -EdgePartition1D; Greedy -Greedy-Oblivious; HCP -HybridCutPlus; RVC -RandomVertexCut).

Both for expansion and modularity could also be defined normalized to the number of edges rather than to the number of vertices.

Greedy-Oblivious (implemented by us), Hybrid-Cut and HybridCutPlus (both implemented by Larry Xiao [START_REF]Hybrid-Cut and HybridCutPlus partitioners[END_REF]).

In each experiment we first randomly selected a partitioner and a graph processing algorithm. Then, we applied the selected partitioner to the input graph. Finally, we executed the selected graph processing algorithm. We repeated every combination at least 30 times to obtain 95% confidence intervals for the execution time whose relative amplitude is always less than 10%.

Our experiments confirm that, as expected, the best balance, and STD metrics are obtained by random partitioners like RandomVertexCut and CanonicalRandomVertexCut. In terms of communication metrics, Fig. 1a shows that DFEP outperforms the other partitioners. This improvement comes at the cost of a much larger partitioning time, indeed in our setting DFEP partitions the graph in a few minutes, larger than the other partitioners considered by a factor 20. Moreover, these partitioning metrics are not necessarily good predictors for the final execution time. Indeed, Fig. 1b shows that while PageRank achieves the shortest execution time when the graph is partitioned by DFEP, Hybrid-Cut provides the best partitioning for Connected Components and HybridCutPlus performs almost as DFEP. This result questions the extent to which partitioning metrics can be used to predict the actual quality of a partition. From the practical point of view, Hybrid-Cut appears to achieve the best trade-off between partitioning time and reduction of the execution time.

Conclusion

In this paper we provided an overview of the existing edge partitioners. An analysis of the related literature shows they have been evaluated considering a disparate set of metrics and, moreover, they have been limitedly compared against each other. We presented some preliminary experimental results comparing a large number of edge partitioners for GraphX framework. We observed that partitioning metrics are not always suited to predict which partitioner will provide the shortest execution time and that simple and fast partitioners like Hybrid-Cut can outperform more sophisticated ones.

Acknowledgment

This work was partly funded by the French Government (National Research Agency, ANR) through the "Investments for the Future" Program reference #ANR-11-LABX-0031-01. We thank the authors of [START_REF] Guerrieri | Dfep: Distributed funding-based edge partitioning[END_REF] to have provided us their implementation of DFEP for Spark.